Supersonic Routes provides a
quick end-to-end introduction to using the himach
package
and is the place to start. This vignette provides advice on more
advanced use, explaining details that the introduction skates over.
Much of this vignette is optional or for occasional use, but the advice on saving and reading the cache is likely to be essential for a speedy workflow.
#the libraries needed for the vignette are
library(himach)
library(dplyr, quietly = TRUE, warn.conflicts = FALSE)
library(ggplot2)
library(sf)
# and we'll load a full set of test data
<- hm_get_test("coast")
NZ_coast <- hm_get_test("buffer")
NZ_buffer30 <- hm_get_test("nofly")
NZ_Buller_buffer40 <- hm_get_test("grid")
NZ_grid <- hm_get_test("route") NZ_routes
himach
uses caching to speed things up. Legs are cached
in route_cache
and arrival-departure links to airports are
cached in star_cache
(STAR is short for standard arrival
route, its counterpart being SID the standard instrument departure).
You will want to save and load the cache (meaning the combination of
route_cache
and star_cache
) as part of your
standard workflow. Quite where you save it is up to you, but a set of
routes is sensitive to (a) the route grid on which it is calculated (b)
the list of aircraft used, and their performance. The saving function
hm_save_cache
forces you to refer to these two datasets,
and uses metadata from them in the file name for the cache.
If you change either of these, then you can use
hm_clean_cache()
to empty the cache. You will also note
that if you run find_routes
and the map has changed, or
findToCToD
and the map or aircraft have changed, then the
cache will be cleared automatically.
For the vignette, we save to a temporary directory. You really don’t want to do this in practice ;-)
hm_clean_cache() #start without cache
# need to load some of the built-in data for this example
<- make_aircraft(warn = FALSE)
aircraft <- make_airports(crs = crs_Pacific)
airports #> Using default airport data: airportr::airport.
options("quiet"= 2) # for a little reporting
# how long does it take with an empty cache?
system.time(
<- find_route(aircraft[1, ],
routes make_AP2("NZAA", "NZDN", airports),
fat_map = NZ_buffer30,
route_grid = NZ_grid,
ap_loc = airports)
)#> Route:-NZAA<>NZDN----
#> Leg: NZAA<>NZDN Aircraft: SST M2.2
#> Cut envelope from lattice: 0.4
#> Running bidirectional Dijkstra...
#> Calculated phase changes
#> Done recursion
#> Checking Shortcuts
#> user system elapsed
#> 1.482 0.041 1.633
# test saving of cache to a disposable directory
<- tempdir()
tmp_dir # for convenience, hm_save_cache gives the full name, including path
<- hm_save_cache("test_v", NZ_grid, aircraft, path = tmp_dir)
full_filename
#empty cache - just to demonstrate the re-loading
# this isn't part of your normal workflow!
hm_clean_cache()
# but normally a session will begin with loading a cache like this
hm_load_cache(full_filename)
# how long does it take with a cache?
system.time(
<- find_route(aircraft[1, ],
routes make_AP2("NZAA", "NZDN", airports),
fat_map = NZ_buffer30,
route_grid = NZ_grid,
ap_loc = airports)
)#> Route:-NZAA<>NZDN----
#> user system elapsed
#> 0.045 0.003 0.055
# if you want to see a map
# map_routes(NZ_coast, routes, crs_Pacific, fat_map = NZ_buffer30, simplify_km = 2)
The cache just works invisibly in the background - you will notice it
speeds up finding of routes no end: in that example, from 1.5s (user) to
0.04s (user) on my machine. In particular, it helps with refuelling,
because the route_cache
quickly remembers the routes from
major hub airports to the main refuelling points, so they don’t need to
be calculated again.
Incidentally, if you add a new refuelling point, then the cache
remains valid because only legs are cached, not
routes. With a new refuelling point, find_route
will check both old legs and new (to the new refuelling points), gaining
where the legs are cached, before selecting the best combination of legs
to make the route.
There may be times when you have built a large cache, of several
thousand routes for a handful of aircraft. If you want to make a change
to the performance of one aircraft, say, you either need to
hm_clean_cache()
and start everything again, or do a bit of
housekeeping. The same is true if an airport location needs correcting:
ideally we delete all references to just that airport.
If you’re not comfortable housekeeping, then just let
himach
run overnight to redo the whole thing.
Housekeeping involves loading the cache manually, with
load
. It contains two environments route_cache
and star_cache
(segments to and from airports). The items
in the cache are named with strings which capture a number of run
parameters, including the aircraft ID. So we identify items for
deletion, and rm
them. Then save the cache manually with
save
.
# for this example, add a second route to the cache
<- find_route(aircraft[2, ],
routes make_AP2("NZAA", "NZDN", airports),
fat_map = NZ_buffer30,
route_grid = NZ_grid,
ap_loc = airports)
#> Route:-NZAA<>NZDN----
#> Leg: NZAA<>NZDN Aircraft: SST M1.6 8.8
#> Cut envelope from lattice: 0.4
#> Running bidirectional Dijkstra...
#> Calculated phase changes
#> Done recursion
#> Checking Shortcuts
# save the cache, which has NZAA-NZDN for 2 aircraft now
hm_save_cache("test_v", NZ_grid, aircraft, path = tmp_dir)
#now do housekeeping
load(full_filename) # filename from the previous chunk
ls(route_cache) # show the contents, just for information
#> [1] "M16_9-NZAA-NZDN--TRUE-TRUE-NA-TRUE-100-12"
#> [2] "M22-NZAA-NZDN--TRUE-TRUE-NA-TRUE-100-12"
# we want to delete instances of aircraft with ID that includes 'M22'
<- ls(route_cache, pattern="M22") |> as.list()
z length(route_cache) # before deletion
#> [1] 2
do.call(rm, z, envir = route_cache) # delete the M22 items
length(route_cache) #after deletion, 1 less
#> [1] 1
# then repeat for star_cache
<- ls(star_cache, pattern="M22") |> as.list()
z length(star_cache)
#> [1] 4
do.call(rm, z, envir = star_cache)
length(star_cache)
#> [1] 2
# then save the result (you might want to change the filename, or backup the old cache beforehand)
save("route_cache", "star_cache", file = full_filename)
It is not unusual for parts of the airspace to be closed, or be
considered unsafe for flying. himach
allows regions to be
marked as ‘avoid’. They will not feature in the grid, so routes will
avoid them, with one exception: an arrival or departure airport can be
inside a no-fly zone, as long as at least one connection point to the
grid is outside. So they might more precisely be called ‘no-overfly’
zones.
A no-fly zone is prepared in the same way as a map of land. If specific countries are to be avoided, this is where having a country name in the geographic data comes in handy.
One essential item is the avoid
attribute of the no-fly
zone. This is used to distinguish sets of legs with different, or no,
no-fly zone. Set
attr(your_avoid_map, "avoid") <- "your summary of that avoid map"
which will (a) remind you what was used (b) tell himach
to
recalculate all legs that have not already been calculated with that
value of avoid
. If you were to add an avoid area for North
Korean airspace, say, then in reality North Atlantic routes are not
affected, but currently himach
plays safe and assumes that
they are.
In this example, no offence is intended to the citizens of Buller District of New Zealand; it is a convenient example for showing how routes are forced to change when airspace is unavailable.
# using your own shp file
# NZ_Buller <- sf::read_sf("...../territorial-authority-2020-clipped-generalised.shp") %>%
# filter(TA2020_V_1 == "Buller District")
# NZ_Buller_u <- sf::st_union(sf::st_simplify(NZ_Buller, dTolerance = 1000))
# NZ_Buller_buffer50 <- sf::st_union(sf::st_buffer(NZ_Buller_u, 50 * 1000))
# attr(NZ_Buller_buffer50, "avoid") <- "Buller+50km"
# the quicker version, using a built-in no fly zone
# this uses data as in the previous code chunk
<- make_aircraft(warn = FALSE)
aircraft <- make_airports(crs = crs_Pacific)
airports #> Using default airport data: airportr::airport.
# run the same route, but with the avoid region
options("quiet"= 2) #just the progress bar
<- aircraft[c(1, 4), ]$id
ac <- find_routes(ac,
routes data.frame(ADEP = "NZAA", ADES = "NZDN"),
aircraft, airports,fat_map = NZ_buffer30,
route_grid = NZ_grid,
cf_subsonic = aircraft[3,],
avoid = NZ_Buller_buffer40)
#> Route:-NZAA<>NZDN----
#> Leg: NZAA<>NZDN Aircraft: SST M2.2
#> Cut envelope from lattice: 0.7
#> Running bidirectional Dijkstra...
#> Calculated phase changes
#> Done recursion
#> Checking Shortcuts
#> Adding subsonic, without range bounds.
#> Leg: NZAA<>NZDN Aircraft: 777-300ER
#> Running bidirectional Dijkstra...
#> Calculated phase changes
#> Done recursion
#> Checking Shortcuts
#>
#> Route:-NZAA<>NZDN----
#> Too far for one leg.
#> Adding subsonic, without range bounds.
#>
#this shows versions of the legs with and without no-fly
# ls(route_cache, pattern = "NZCH", envir = .hm_cache)
# create route summary
<- summarise_routes(routes, airports)
rtes
# draw a basic map
map_routes(NZ_coast, routes, crs_Pacific, fat_map = NZ_buffer30,
avoid_map = NZ_Buller_buffer40,
simplify_km = 2)
map_routes(NZ_coast, routes, show_route = "aircraft",
crs = crs_Pacific, fat_map = NZ_buffer30,
avoid_map = NZ_Buller_buffer40,
simplify_km = 2)
NA
routes?After a call to find_routes
, the output can have
NA
entries in some columns for some routes. There are two
reasons for this:
refuel = xxx
, then you will find other entries for the same
routeID
(eg “EGLL<>KSFO”) but with different
fullRouteID
(eg “EGLL<>PANC<>KSFO”) showing a
good route including refuelling.find_routes
.So these appear when the specified route is not possible.
Above 60 or 70 (North or South), the approximations used by the
st_buffer
function begin to show signs of exceeding their
limits. In particular, if you’re adding a 50km coastal buffer, for
example, there are separations between Canadian islands which are just
under 100km. Borden and Ellef Ringnes are examples. A buffer generated
by st_buffer
shows the strait between them as open water,
where it should be closed.
This can lead to over-optimistic routings: supersonic where they should not be.
The solution is to use the links from sf
to the
s2
package which come in more recent versions of the
sf
package. This does require you to use quite a high value
for the max_cells
parameter of
s2::s2_buffer_cells
.
<- s2::s2_data_countries(c("Greenland", "Canada", "Iceland"))
gr <- s2::s2_buffer_cells(gr, distance = 50000, max_cells = 20000) %>%
gr_buffer_s2 st_as_sfc()
<- ggplot(st_transform(gr_buffer_s2, crs_Atlantic)) + geom_sf(fill = "grey40") +
m_s2 geom_sf(data = st_transform(st_as_sfc(gr), crs_Atlantic))
sf_use_s2(FALSE) # to be sure
<- gr %>%
gr_transf st_as_sfc() %>%
st_transform(crs_Atlantic)
<- gr_transf %>%
gr_t_buffer st_buffer(dist = 50000)
<- ggplot(gr_t_buffer) + geom_sf(fill = "grey40") + geom_sf(data = gr_transf)
m_old
::plot_grid(m_old, m_s2, labels = c("bad", "good"),
cowplotncol = 1)
In fact, the problem of finding too many apparently over-ocean routes is broader than this. The other main contributor to this is missing islands from the map. See the comments in the first vignette.
An example of this is in the same place. Some maps omit small islands (well, larger ones like Killniq down to tiny ones like Goodwin Island) at the mouth of the Hudson Strait. This affects the apparent width of the opening. Given the islands, and a 50km buffer, the Strait is not open as the next example shows.
This uses a non-CRAN, but public package of hi-resolution maps,
rnaturalearthhires
. If you don’t want to load this package,
just note the results shown in the figure.
::sf_use_s2(TRUE)
sf<- sf::st_as_sf(rnaturalearthhires::countries10) %>%
hires filter(NAME %in% c("Greenland", "Canada", "Iceland"))
<- s2::s2_buffer_cells(hires, distance = 50000, max_cells = 20000) %>%
hires_buffer_s2 st_as_sfc()
<- ggplot(st_transform(hires_buffer_s2, crs_Atlantic)) +
m_hires geom_sf(fill = "grey40") +
geom_sf(data = st_transform(hires, crs_Atlantic))
::plot_grid(m_s2, m_hires, labels = c("good", "better"),
cowplotncol = 1)
There are a number of place in the vignettes, eg making an airport
dataset, where we have shown the use of a parameter to specify a
coordinate reference system. himach
has recently
transitioned to using spherical geometry directly using the
s2
package, both directly and through the sf
package. Before s2
was available in sf
there
was a constant need to align the coordinate reference systems of objects
before combining them.
Now, in theory, all geometrical operations use spherical geometry, so
a coordinate reference system should only be needed when you plot a map.
At that point, the coordinate reference system is saying how to move
from spherical coordinates to a flat projection. Four basic projections
are supplied crs_Atlantic
, crs_Pacific
,
crs_North
and crs_South
which you can use in
map_routes
to get the right map for your particular set of
routes. You can create others as shown in the vignette.
We will remove remaining references to coordinate reference systems
during route creation in later versions of himach
.