Skip to content

Areas of interest

Every public function takes its area of interest in any of these forms:

ef.composite((-111.9, 40.7, -111.8, 40.8), ...)   # bbox tuple (WGS84)
ef.composite("Yosemite National Park", ...)        # place name (geocoded)
ef.composite("watershed.geojson", ...)             # GeoJSON file
ef.composite("/data/fields.shp", ...)              # shapefile (or .gpkg, .kml)
ef.composite(geojson_dict, ...)                    # GeoJSON dict
ef.composite(shapely_polygon, ...)                 # __geo_interface__
  • bbox — (min_lon, min_lat, max_lon, max_lat) in WGS84 degrees.
  • place name — resolved with OpenStreetMap Nominatim (free, no key).
  • GeoJSON / shapely — Polygon, MultiPolygon, Point, LineString, and collections. Degenerate geometries (a Point, or an axis-aligned line) resolve to a small valid bbox around the feature.

Vector files (shapefile, GeoPackage, KML)

Point any function at a local vector file. GeoJSON works out of the box; shapefiles, GeoPackages, and KML need the vector extra (pip install "earthfetch[vector]"):

ef.load_dem("/data/basin.shp", ...)        # shapefile
ef.composite("/data/parcels.gpkg", ...)    # GeoPackage

All features in the file are unioned into one AOI, and — importantly — the file's CRS is read and reprojected to WGS84 automatically, so a shapefile stored in a UTM or State-Plane projection lands in the right place without any manual step.

Clipping

crs="utm" (the default for composite, terrain, load_naip) picks the right UTM zone for metric pixels. Explicit polygons clip results to their boundary; geocoded place names return the full rectangle (pass clip=True to cut to the boundary).

earthfetch.aoi.resolve_aoi 

resolve_aoi(aoi) -> AOI

Normalize any supported AOI input to an AOI (see module docs).

Source code in src/earthfetch/aoi.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
def resolve_aoi(aoi) -> AOI:
    """Normalize any supported AOI input to an ``AOI`` (see module docs)."""
    if isinstance(aoi, AOI):
        return aoi
    if hasattr(aoi, "__geo_interface__"):
        return _from_geojson(dict(aoi.__geo_interface__))
    if isinstance(aoi, dict):
        return _from_geojson(aoi)
    if isinstance(aoi, str) or isinstance(aoi, os.PathLike):
        text = str(aoi)
        low = text.lower()
        if low.endswith((".geojson", ".json")):
            return _from_geojson(json.loads(Path(text).read_text()))
        from .vector import VECTOR_EXTENSIONS
        if low.endswith(VECTOR_EXTENSIONS):
            from .vector import read_vector
            return read_vector(text)
        if Path(text).is_file():
            # a local file with an unfamiliar name — try GeoJSON, else vector
            try:
                return _from_geojson(json.loads(Path(text).read_text()))
            except (ValueError, UnicodeDecodeError):
                from .vector import read_vector
                return read_vector(text)
        return geocode(text)
    if isinstance(aoi, Sequence) and len(aoi) == 4:
        return AOI(bbox=validate_bbox(aoi))
    raise EarthfetchError(
        f"cannot interpret AOI of type {type(aoi).__name__}: pass a bbox "
        "tuple, GeoJSON, a .geojson/.shp/.gpkg path, a shapely geometry, "
        "or a place name"
    )

earthfetch.aoi.geocode 

geocode(place: str) -> AOI

Resolve a place name to a bbox via OpenStreetMap Nominatim (free).

Source code in src/earthfetch/aoi.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def geocode(place: str) -> AOI:
    """Resolve a place name to a bbox via OpenStreetMap Nominatim (free)."""
    resp = get_session().get(
        NOMINATIM_URL,
        params={"q": place, "format": "json", "limit": 1, "polygon_geojson": 1},
        timeout=30,
    )
    resp.raise_for_status()
    results = resp.json()
    if not results:
        raise EarthfetchError(f"place not found: {place!r}")
    hit = results[0]
    south, north, west, east = (float(v) for v in hit["boundingbox"])
    geometry = hit.get("geojson")
    if geometry and geometry.get("type") not in ("Polygon", "MultiPolygon"):
        geometry = None  # points/lines can't mask an area
    logger.info("geocoded %r -> %s (%s)", place, hit.get("display_name"),
                (west, south, east, north))
    return AOI(bbox=validate_bbox((west, south, east, north)),
               geometry=geometry, name=hit.get("display_name"),
               clip_default=False)

earthfetch.aoi.utm_crs 

utm_crs(bbox: Sequence[float]) -> str

EPSG code of the UTM zone at a bbox center, e.g. 'EPSG:32612'.

Source code in src/earthfetch/aoi.py 
151
152
153
154
155
156
157
def utm_crs(bbox: Sequence[float]) -> str:
    """EPSG code of the UTM zone at a bbox center, e.g. 'EPSG:32612'."""
    min_lon, min_lat, max_lon, max_lat = validate_bbox(bbox)
    lon = (min_lon + max_lon) / 2
    lat = (min_lat + max_lat) / 2
    zone = min(60, max(1, math.floor((lon + 180) / 6) + 1))
    return f"EPSG:{(32600 if lat >= 0 else 32700) + zone}"