Choosing an area of interest¶
Every backend needs to know where. There are two channels, and aoi= accepts
seven different shapes — five of which appear in no other example.
lat_lim / lon_lim are the declared parameters, both [min, max] in degrees.
aoi= is the flexible channel: whatever you pass is normalised down to that same
pair, and a real polygon is kept alongside it as a clip mask.
This notebook shows each form resolving to the bbox it produces. GEBCO is used throughout because it is anonymous and small — no account, no API key.
Setup¶
resolve_aoi is the function every backend uses to interpret aoi=. Calling it
directly is the clearest way to see what each form means, with no download involved.
from earthlens.base.spatial import resolve_aoi
def show(label, **kwargs):
lat_lim, lon_lim, geometry = resolve_aoi(**kwargs)
mask = 'polygon mask' if geometry is not None else 'bbox only'
print(f'{label:24} lat {lat_lim} lon {lon_lim} ({mask})')
1. lat_lim / lon_lim — the declared pair¶
The base form. Both are [min, max] in degrees: lat_lim=[south, north],
lon_lim=[west, east]. Every other form below is reduced to exactly this.
lat_lim = [36.2, 38.0]
lon_lim = [-29.5, -27.7]
print(f'lat {lat_lim} lon {lon_lim}')
2. aoi= as a bounding box¶
Four values in [W, S, E, N] order — note that is not the same order as
the lat_lim / lon_lim pair above.
show('bbox [W, S, E, N]', aoi=[-29.5, 36.2, -27.7, 38.0])
3. aoi= as a point, with buffer=¶
Two values are read as (lon, lat). A point has no area, so buffer= (a
half-width in degrees) is required — omitting it raises rather than
assuming a default. The result is clamped to the valid lon/lat range.
show('point + buffer', aoi=(-28.6, 37.1), buffer=0.9)
try:
resolve_aoi(aoi=(-28.6, 37.1))
except ValueError as exc:
print(f'without buffer -> ValueError: {exc}')
4. aoi= as WKT¶
A WKT string is parsed with shapely. Because it describes a real shape, the polygon is kept as a clip mask, not just its bounds.
wkt = 'POLYGON((-29.5 36.2, -27.7 36.2, -27.7 38.0, -29.5 38.0, -29.5 36.2))'
show('WKT polygon', aoi=wkt)
5. aoi= as GeoJSON¶
A GeoJSON mapping works the same way — a Feature, a bare geometry, or a dict
carrying a bbox key.
geojson = {
'type': 'Polygon',
'coordinates': [
[[-29.5, 36.2], [-27.7, 36.2], [-27.7, 38.0], [-29.5, 38.0], [-29.5, 36.2]]
],
}
show('GeoJSON polygon', aoi=geojson)
6. aoi= as a bbox mapping¶
A mapping that is not GeoJSON is read as the four bbox edges instead, and the
key spelling is flexible — GeoJSON min_lon, eodag lonmin,
shapely/geopandas minx and compass west all work, matched
case-insensitively. This is the form to use when the box arrives as another
tool's JSON rather than as a list.
show(
'bbox mapping (min_lon)',
aoi={'min_lon': -29.5, 'min_lat': 36.2, 'max_lon': -27.7, 'max_lat': 38.0},
)
show(
'bbox mapping (compass)',
aoi={'West': -29.5, 'South': 36.2, 'East': -27.7, 'North': 38.0},
)
7. aoi= as a GeoDataFrame¶
A geopandas GeoDataFrame or GeoSeries is accepted — the check is
duck-typed on total_bounds, so geopandas is never imported just to test a type.
Two things happen for you:
- the frame is reprojected.
total_boundsis in the frame's own CRS, so a UTM or Web Mercator frame would otherwise yield a metre-valued, out-of-range bbox. A frame with no CRS is taken as already lon/lat. - a polygonal frame is kept as a clip mask; a points or lines frame contributes its bounding box only.
import geopandas as gpd
from shapely.geometry import Polygon
azores = Polygon([(-29.5, 36.2), (-27.7, 36.2), (-27.7, 38.0), (-29.5, 38.0)])
gdf = gpd.GeoDataFrame(geometry=[azores], crs='EPSG:4326')
show('GeoDataFrame (4326)', aoi=gdf)
# the same polygon in Web Mercator resolves to the same lon/lat box
show('GeoDataFrame (3857)', aoi=gdf.to_crs(3857))
# a points frame has no area to clip to, so only its bounds are used
points = gpd.GeoDataFrame(
geometry=gpd.points_from_xy([-29.5, -27.7], [36.2, 38.0]), crs='EPSG:4326'
)
show('GeoDataFrame (points)', aoi=points)
8. aoi= as a shapely geometry¶
Anything exposing __geo_interface__ is accepted, which covers shapely
geometries directly.
show('shapely polygon', aoi=azores)
Using one in a real request¶
Any of the forms above can go straight into a backend. GEBCO is anonymous, so this runs without credentials.
from earthlens.core import EarthLens
# load() hands back the fetched rasters as pyramids objects, one per written
# file, so a single-tile request is a one-element list.
[bathymetry] = EarthLens(
data_source='gebco',
dataset='gebco_2020',
aoi=azores,
).load()
stats = bathymetry.stats(approx_ok=False)
print('epsg', bathymetry.epsg, '| shape', bathymetry.shape)
print(
'depth range',
round(float(stats['min'].iloc[0])),
'..',
round(float(stats['max'].iloc[0])),
'm',
)
bathymetry.plot(cmap='gist_earth', title='GEBCO 2020 - Azores')
Two things that will bite you¶
A polygon is only honoured if the backend supports it. resolve_aoi always
returns the mask, but a backend whose SUPPORTS_POLYGON_AOI is False uses the
bounding box and emits a PolygonAoiWarning. That warning's own docstring says
why it exists: the download still succeeds, it just covers the bbox, and that is
"the most dangerous kind of wrong result — a valid raster of the right variable
over roughly the right area". Watch for it whenever you pass a polygon.
Antimeridian crossings are rejected, not guessed. If west ends up east of east, the request raises and tells you to split at ±180 and issue the two halves separately, rather than silently interpreting the box.
try:
resolve_aoi(aoi=[170.0, -20.0, -170.0, -10.0])
except ValueError as exc:
print(f'ValueError: {exc}')