Named regions

A study area has to be positioned somewhere. You can give a bounding box, or your own shapefile - or you can name the ground: "Scotland", "Madagascar", "GREAT BRITAIN". The names come from Natural Earth, which is in the public domain, at 1:10m.

Two things use them. EcoSISTEM.boundingbox gives a name's extent, read from a table shipped with the package at no download. NaturalEarthSpec gives its actual outline, fetching the polygons the first time it is built. Both resolve a name by the same rule, so the box that one reports is the box the other's shape has.

A name needs its level

A name on its own often means more than one thing, so every region is identified by a level as well as a name.

using EcoSISTEM, Unitful, Unitful.DefaultSymbols
length(EcoSISTEM.naturalearth_levels())

31

The levels divide into political ones (ADMIN a country, SUBUNIT its constituent countries), statistical groupings (REGION_UN, SUBREGION), a code (ISO_A3_EH), and Natural Earth's physical vocabulary, prefixed Physical - Physical Island, Physical Continent, Physical Desert and the rest.

Where the levels agree you need not say which you mean:

EcoSISTEM.boundingbox("Scotland")

Extent(Y = (54.632°, 60.848°), X = (-13.692°, -0.754°))

Where they disagree, naming one is required rather than guessed at, and the error tabulates the choice:

try
    EcoSISTEM.boundingbox("Africa")
catch e
    println(sprint(showerror, e))
end

"Africa" means different regions at different levels, and the extents they give for AllTerritories() disagree:
    CONTINENT           Africa            W   -25.361  S  -46.966  E    51.418  N   37.346
    REGION_UN           Africa            W   -25.361  S  -53.193  E   105.715  N   37.346
    Physical Continent  AFRICA            W   -17.498  S  -34.786  E    51.381  N   37.318
Name one with `level = "..."`.

EcoSISTEM.naturalearth_regions(level) lists what exists at one, with each region's box and area:

EcoSISTEM.naturalearth_regions("Physical Tundra")

3 named regions at this level, by name:
  name                            W       S       E       N    area/km2 parts largest
  BARREN GROUNDS              -125.99   59.67   -87.0    70.1      830486     2
  Bol’shezemel’skaya Tundra     47.08   64.63   66.63   69.85      291911     2
  CANADIAN SHIELD             -119.73    42.8   -55.7    69.9     4459950    17

Passing a region where a level belongs is a common slip, and the error says which level that region is actually defined at rather than only that no such level exists.

How much of a name to take

A name almost never denotes one connected piece of ground. coverage says how much to take, and the default is everything the name covers, because that is what the source means by it.

france = EcoSISTEM.boundingbox("France", level = "ADMIN")

Extent(Y = (-21.371°, 51.088°), X = (-61.798°, 55.855°))

That reaches from French Guiana in South America to Mayotte in the Indian Ocean. For the ground most people mean, ask for the principal landmass:

EcoSISTEM.boundingbox("France", level = "ADMIN", coverage = LargestLandmass())

Extent(Y = (42.325°, 51.088°), X = (-4.785°, 8.201°))

LargestLandmass(count = n) takes the largest n - New Zealand's two main islands, Japan's four - and LandmassesAbove(threshold) takes everything clearing a threshold, which is how the United Kingdom loses Rockall. The threshold is either an area or a share of the region's own total:

LandmassesAbove(1km^2)     # an absolute size
LandmassesAbove(5percent)  # ...or a fraction of the region, which travels better between regions

A share must be a percentage rather than a bare number. 0.05 and 5percent are the same quantity, but only one of them says which it means when read beside 1km^2, so the bare number is refused.

Finding a region

investigate_regions asks which named regions relate to something you already have - a coordinate, a raster, a layer, a study area:

EcoSISTEM.investigate_regions(LatLong(55.95°, -3.19°), kind = :political, limit = 5)

5 named regions enclose it, smallest first:
  level       name                 W       S       E       N    area/km2 parts largest
  GEOUNIT     Scotland          -13.69   54.63   -0.75   60.85       78303    47    88%
  SUBUNIT     Scotland          -13.69   54.63   -0.75   60.85       78303    47    88%
  ADMIN       United Kingdom    -13.69   49.91    1.77   60.85      243783    57    90%
  SOVEREIGNT  United Kingdom   -130.75  -59.47    72.5   60.85      261906   141    83%
  SOVEREIGNT  Norway             -9.12  -54.46   33.64   80.77      382074   120    79%

Compared by bounding box, which costs no download but is coarse - Chile's box 
spans 43 degrees because of Easter Island. Pass `exact = true` to check against
the real outlines instead.

Encloses is the default. Overlaps(x) asks which regions your data reaches into, ordered by how much ground they share; Within(x) asks which lie entirely inside it, which is the question "what can I simulate in full with the data I have?"

A row of the report can be turned straight into a spec. Use only when the answer should be unique - it refuses if it is not - and first to take the best by the report's own ordering:

match = only(EcoSISTEM.investigate_regions(LatLong(55.95°, -3.19°), level = "SUBUNIT"))
NaturalEarthSpec(match)

NaturalEarthSpec("Scotland", level = "SUBUNIT")

Boxes, and then outlines

By default the query compares bounding boxes, which is what makes it free - and it is loose. Norway's box encloses Edinburgh, running west to Jan Mayen and north to Svalbard, so a query at that coordinate lists Norway beside the United Kingdom.

exact = true checks the survivors against the regions' real outlines instead. It downloads the geometry, so it is not free - but it removes those false positives, and it reaches the 54 regions that cross the antimeridian, which have no box a query can compare at all:

investigate_regions(LatLong(55.95°, -3.19°), kind = :political, exact = true)
# Norway is gone; CONTINENT Europe, which wraps, appears

Refinement is lazy and in box order, stopping as soon as the answer cannot change - refining only ever removes a match or shrinks its overlap, so a confirmed limit cannot be displaced by anything later. A continental query that matches 361 regions by box typically fetches a few dozen, and the report says how many.

Using a region as a study area

NaturalEarthSpec is a within mask, so it both restricts which cells are simulated and sets the grid's extent:

area = StudyArea(regime = temperature, supply = rainfall,
                 within = NaturalEarthSpec("Scotland", coverage = LargestLandmass()),
                 crs = EPSG(27700), cellsize = 1.0km)

outline = false takes the region's bounding box instead of its coastline, which is the cheaper thing to want when the region is only saying where to work.

Combining regions

Regions compose as geometry, through ConstructedShapeSpec - the vector mirror of ConstructedRasterSpec. The result is exact and carries no resolution of its own, so the study grid is still free to be decided afterwards.

uk = NaturalEarthSpec("United Kingdom", level = "ADMIN")
ireland = NaturalEarthSpec("Ireland", level = "ADMIN")
man = NaturalEarthSpec("Isle of Man", level = "ADMIN")
ConstructedShapeSpec(ShapeUnion(), uk, ireland, man)

ConstructedShapeSpec(ShapeUnion(), NaturalEarthSpec("United Kingdom", level = "ADMIN"), NaturalEarthSpec("Ireland", level = "ADMIN"), NaturalEarthSpec("Isle of Man", level = "ADMIN"))

Members may be any shape spec, including a ShapeSpec of your own study area and another ConstructedShapeSpec. Operations are ShapeUnion, ShapeIntersection and ShapeDifference for combining two or more, and ShapeBuffer, ShapeSimplify and ShapeConvexHull for transforming exactly one:

scotland = NaturalEarthSpec("Scotland", coverage = LargestLandmass())
ConstructedShapeSpec(ShapeBuffer(50km), scotland)

ConstructedShapeSpec(ShapeBuffer(50 km), NaturalEarthSpec("Scotland", level = "GEOUNIT", coverage = LargestLandmass(count = 1)))

An arbitrary function is accepted too, exactly as ConstructedRasterSpec takes any combine, so anything the geometry library offers is reachable without a new operation.

If you want X, use...

These are cartographic outlines

Natural Earth draws its shapes to be printed at a scale, not to be authoritative boundaries, and this bites in practice rather than in principle. Two measured examples:

• Its own BRITISH ISLES polygon stops at 59.80°N, so it omits Shetland. The union of the United Kingdom, Ireland and the Isle of Man reaches 60.85°N and does not.
• Its physical continents are landmass outlines, so EUROPE does not contain the British Isles at all - Paris is inside it, central England is not. An intersection with it is empty, and is refused rather than returned as a grid with nothing active.

A third to know about: the physical Channel Islands are the Californian ones.

Use a named region as a convenience. If your study needs a particular boundary, supply it as a ShapeSpec of your own - and combine it with a named one if that helps.