Data sources

All sources are free and need no API key or account.

USGS 3DEP DEMs

The National Map, United States only. Resolutions 1m, 10m, 30m, and 5m-ak (Alaska IFSAR).

earthfetch.usgs.search_dem

search_dem(bbox: Sequence[float], resolution: str = '10m', max_items: int = 100) -> list[dict]

Find 3DEP DEM tiles intersecting a bbox.

Parameters:

Name	Type	Description	Default
bbox	(min_lon, min_lat, max_lon, max_lat) in WGS84 degrees.		required
resolution	one of ``DEM_DATASETS`` keys ("1m", "10m", "30m", "5m-ak").		'10m'
max_items	cap on returned tiles.		100

Returns:

Type	Description
list of dict	Product dicts with title, downloadURL, sizeInBytes, boundingBox and other TNM metadata.

Source code in src/earthfetch/usgs.py

def search_dem(
    bbox: Sequence[float],
    resolution: str = "10m",
    max_items: int = 100,
) -> list[dict]:
    """Find 3DEP DEM tiles intersecting a bbox.

    Parameters
    ----------
    bbox : (min_lon, min_lat, max_lon, max_lat) in WGS84 degrees.
    resolution : one of ``DEM_DATASETS`` keys ("1m", "10m", "30m", "5m-ak").
    max_items : cap on returned tiles.

    Returns
    -------
    list of dict
        Product dicts with ``title``, ``downloadURL``, ``sizeInBytes``,
        ``boundingBox`` and other TNM metadata.
    """
    if resolution not in DEM_DATASETS:
        raise ValueError(f"resolution must be one of {sorted(DEM_DATASETS)}")
    bbox = validate_bbox(bbox)

    session = get_session()
    items: list[dict] = []
    offset = 0
    while len(items) < max_items:
        params = {
            "datasets": DEM_DATASETS[resolution],
            "bbox": ",".join(str(v) for v in bbox),
            "prodFormats": "GeoTIFF",
            "outputFormat": "JSON",
            "max": min(100, max_items - len(items)),
            "offset": offset,
        }
        resp = session.get(TNM_URL, params=params, timeout=60)
        resp.raise_for_status()
        data = resp.json()
        batch = data.get("items", [])
        items.extend(batch)
        offset += len(batch)
        if not batch or offset >= data.get("total", 0):
            break
    logger.info("TNM: %d %s tile(s) for bbox %s", len(items), resolution, bbox)
    return items[:max_items]

earthfetch.usgs.download_dem

download_dem(bbox: Sequence[float], resolution: str = '10m', out_dir: str | PathLike | None = None, max_items: int = 100, overwrite: bool = False, workers: int = 4, progress: ProgressFn | None = None) -> list[Path]

Download DEM tiles covering a bbox in parallel. Returns local paths.

Tiles already on disk are skipped, so calls are resumable. out_dir defaults to the earthfetch cache. Raises TileNotFoundError when no tiles cover the bbox (e.g. outside the US — try Copernicus instead).

Source code in src/earthfetch/usgs.py

def download_dem(
    bbox: Sequence[float],
    resolution: str = "10m",
    out_dir: str | os.PathLike | None = None,
    max_items: int = 100,
    overwrite: bool = False,
    workers: int = 4,
    progress: ProgressFn | None = None,
) -> list[Path]:
    """Download DEM tiles covering a bbox in parallel. Returns local paths.

    Tiles already on disk are skipped, so calls are resumable. ``out_dir``
    defaults to the earthfetch cache. Raises ``TileNotFoundError`` when no
    tiles cover the bbox (e.g. outside the US — try Copernicus instead).
    """
    urls = dem_tile_urls(bbox, resolution=resolution, max_items=max_items)
    if not urls:
        raise TileNotFoundError(
            f"no USGS {resolution} DEM tiles cover {tuple(bbox)}; "
            "USGS covers the US only — try load_dem(source='copernicus')"
        )
    with ThreadPoolExecutor(max_workers=workers) as pool:
        futures = [
            pool.submit(download_file, url, out_dir=out_dir,
                        overwrite=overwrite, progress=progress)
            for url in urls
        ]
        return [f.result() for f in futures]

Copernicus GLO-30

Global 30 m DEM on AWS Open Data. Used automatically as the fallback outside the US.

earthfetch.copernicus.copernicus_dem_urls

copernicus_dem_urls(bbox: Sequence[float]) -> list[str]

COG URLs for GLO-30 tiles covering a bbox. Ocean-only tiles are absent from the bucket and are silently skipped (checked via HEAD).

Raises TileNotFoundError when no tile exists (open ocean).

Source code in src/earthfetch/copernicus.py

def copernicus_dem_urls(bbox: Sequence[float]) -> list[str]:
    """COG URLs for GLO-30 tiles covering a bbox. Ocean-only tiles are
    absent from the bucket and are silently skipped (checked via HEAD).

    Raises ``TileNotFoundError`` when no tile exists (open ocean).
    """
    min_lon, min_lat, max_lon, max_lat = validate_bbox(bbox)
    session = get_session()
    urls = []
    for lat in range(math.floor(min_lat), math.ceil(max_lat)):
        for lon in range(math.floor(min_lon), math.ceil(max_lon)):
            name = _tile_name(lat, lon)
            url = f"{COP_BUCKET}/{name}/{name}.tif"
            if session.head(url, timeout=30).status_code == 200:
                urls.append(url)
            else:
                logger.debug("no GLO-30 tile at %s (ocean?)", name)
    if not urls:
        raise TileNotFoundError(f"no Copernicus GLO-30 tiles cover {tuple(bbox)}")
    logger.info("Copernicus: %d tile(s) for bbox %s", len(urls), tuple(bbox))
    return urls

Sentinel-2 L2A

Multispectral imagery, global, via the Earth Search STAC API. Bands are addressed by ESA id (B04, B08, SCL, TCI, ...).

earthfetch.sentinel.search_sentinel2

search_sentinel2(bbox: Sequence[float], start: str, end: str, max_cloud: float = 20.0, limit: int = 50) -> list[dict]

Search Sentinel-2 L2A scenes.

Parameters:

Name	Type	Description	Default
bbox	(min_lon, min_lat, max_lon, max_lat) in WGS84 degrees.		required
start	ISO dates, e.g. "2026-05-01" / "2026-06-01".		required
end	ISO dates, e.g. "2026-05-01" / "2026-06-01".		required
max_cloud	maximum scene cloud cover percent.		20.0
limit	cap on returned scenes.		50

Returns:

Type	Description
list of dict	STAC item dicts sorted by cloud cover (clearest first). Each has id, properties (datetime, eo:cloud_cover, ...) and assets.

Source code in src/earthfetch/sentinel.py

def search_sentinel2(
    bbox: Sequence[float],
    start: str,
    end: str,
    max_cloud: float = 20.0,
    limit: int = 50,
) -> list[dict]:
    """Search Sentinel-2 L2A scenes.

    Parameters
    ----------
    bbox : (min_lon, min_lat, max_lon, max_lat) in WGS84 degrees.
    start, end : ISO dates, e.g. "2026-05-01" / "2026-06-01".
    max_cloud : maximum scene cloud cover percent.
    limit : cap on returned scenes.

    Returns
    -------
    list of dict
        STAC item dicts sorted by cloud cover (clearest first). Each has
        ``id``, ``properties`` (datetime, eo:cloud_cover, ...) and ``assets``.
    """
    bbox = validate_bbox(bbox)
    session = get_session()
    body = {
        "collections": [COLLECTION],
        "bbox": list(bbox),
        "datetime": f"{start}T00:00:00Z/{end}T23:59:59Z",
        "limit": min(limit, 100),
        "query": {"eo:cloud_cover": {"lt": max_cloud}},
    }
    items: list[dict] = []
    while len(items) < limit:
        resp = session.post(STAC_URL, json=body, timeout=60)
        resp.raise_for_status()
        page = resp.json()
        items.extend(page.get("features", []))
        nxt = next(
            (lnk for lnk in page.get("links", []) if lnk.get("rel") == "next"), None
        )
        if nxt is None or not page.get("features"):
            break
        body = nxt.get("body", body)
    items = items[:limit]
    items.sort(key=lambda i: i["properties"].get("eo:cloud_cover", 100))
    logger.info("Earth Search: %d scene(s) %s..%s cloud<%s%%",
                len(items), start, end, max_cloud)
    return items

earthfetch.sentinel.clearest_scene

clearest_scene(bbox: Sequence[float], start: str, end: str, max_cloud: float = 20.0) -> dict

The least-cloudy scene in a date range, or raise NoScenesError.

Source code in src/earthfetch/sentinel.py

def clearest_scene(
    bbox: Sequence[float], start: str, end: str, max_cloud: float = 20.0
) -> dict:
    """The least-cloudy scene in a date range, or raise ``NoScenesError``."""
    items = search_sentinel2(bbox, start, end, max_cloud=max_cloud, limit=50)
    if not items:
        raise NoScenesError(
            f"no Sentinel-2 scenes for {tuple(bbox)} in {start}..{end} "
            f"with cloud < {max_cloud}% — widen dates or raise max_cloud"
        )
    return items[0]

earthfetch.sentinel.covering_scenes

covering_scenes(bbox: Sequence[float], start: str, end: str, max_cloud: float = 20.0, limit: int = 100, grid: int = 16) -> list[dict]

Clearest scenes whose footprints together cover the whole bbox.

clearest_scene returns a single scene, which may only partially cover an AOI that straddles Sentinel-2 tile boundaries. This walks scenes clearest-first and keeps each one that adds coverage until the bbox is filled — so the mosaic is built from the clearest available scenes. Feed the result to load_sentinel2(items=...).

Parameters:

Name	Type	Description	Default
bbox	(min_lon, min_lat, max_lon, max_lat) in WGS84 degrees.		required
start	ISO dates bounding the search.		required
end	ISO dates bounding the search.		required
max_cloud	maximum scene cloud cover percent.		20.0
limit	cap on scenes considered.		100
grid	coverage-test resolution (grid x grid sample points over the bbox).		16

Returns:

Type	Description
list of dict	Selected STAC items, clearest first. Raises NoScenesError when no scenes match; returns a partial cover (with a warning) if the scenes found cannot fully cover the bbox.

Source code in src/earthfetch/sentinel.py

def covering_scenes(
    bbox: Sequence[float],
    start: str,
    end: str,
    max_cloud: float = 20.0,
    limit: int = 100,
    grid: int = 16,
) -> list[dict]:
    """Clearest scenes whose footprints together cover the whole bbox.

    ``clearest_scene`` returns a single scene, which may only partially cover
    an AOI that straddles Sentinel-2 tile boundaries. This walks scenes
    clearest-first and keeps each one that adds coverage until the bbox is
    filled — so the mosaic is built from the clearest available scenes. Feed
    the result to ``load_sentinel2(items=...)``.

    Parameters
    ----------
    bbox : (min_lon, min_lat, max_lon, max_lat) in WGS84 degrees.
    start, end : ISO dates bounding the search.
    max_cloud : maximum scene cloud cover percent.
    limit : cap on scenes considered.
    grid : coverage-test resolution (grid x grid sample points over the bbox).

    Returns
    -------
    list of dict
        Selected STAC items, clearest first. Raises ``NoScenesError`` when no
        scenes match; returns a partial cover (with a warning) if the scenes
        found cannot fully cover the bbox.
    """
    items = search_sentinel2(bbox, start, end, max_cloud=max_cloud, limit=limit)
    if not items:
        raise NoScenesError(
            f"no Sentinel-2 scenes for {tuple(bbox)} in {start}..{end} "
            f"with cloud < {max_cloud}%"
        )
    minx, miny, maxx, maxy = validate_bbox(bbox)
    # interior sample points (avoid the exact edges)
    pts = [
        (minx + (maxx - minx) * (i + 0.5) / grid,
         miny + (maxy - miny) * (j + 0.5) / grid)
        for i in range(grid) for j in range(grid)
    ]
    uncovered = set(range(len(pts)))
    # walk scenes clearest-first (then by day, for seam-free mosaics) and keep
    # any that add coverage — so the cover is built from the clearest scenes,
    # never a single cloudy one that happens to span the whole bbox
    candidates = sorted(
        items,
        key=lambda it: (it["properties"].get("eo:cloud_cover", 100),
                        it["properties"].get("datetime", "")),
    )
    chosen: list[dict] = []
    for it in candidates:
        if not uncovered:
            break
        new = {k for k in uncovered if _covers_point(it.get("geometry"), *pts[k])}
        if new:
            chosen.append(it)
            uncovered -= new
    if uncovered:
        logger.warning("covering_scenes: %d/%d sample points uncovered — "
                       "scenes do not fully cover the bbox", len(uncovered), len(pts))
    logger.info("covering_scenes: %d scene(s) cover the bbox", len(chosen))
    return chosen

earthfetch.sentinel.download_sentinel2

download_sentinel2(item: dict, bands: Sequence[str] = ('B04', 'B03', 'B02'), out_dir: str | PathLike | None = None, overwrite: bool = False, workers: int = 4, progress: ProgressFn | None = None) -> dict[str, Path]

Download selected bands of one STAC item as GeoTIFFs, in parallel.

bands accepts ESA ids ("B04", "B08", "TCI", "SCL") or Earth Search asset keys ("red", "nir", "visual"). Files land in out_dir/<scene id>/ (cache dir by default). Returns {band: path}.

Source code in src/earthfetch/sentinel.py

def download_sentinel2(
    item: dict,
    bands: Sequence[str] = ("B04", "B03", "B02"),
    out_dir: str | os.PathLike | None = None,
    overwrite: bool = False,
    workers: int = 4,
    progress: ProgressFn | None = None,
) -> dict[str, Path]:
    """Download selected bands of one STAC item as GeoTIFFs, in parallel.

    ``bands`` accepts ESA ids ("B04", "B08", "TCI", "SCL") or Earth Search
    asset keys ("red", "nir", "visual"). Files land in
    ``out_dir/<scene id>/`` (cache dir by default). Returns {band: path}.
    """
    from .utils import get_cache_dir

    root = Path(out_dir) if out_dir is not None else get_cache_dir() / "downloads"
    base = root / item["id"]
    urls = {band: band_url(item, band) for band in bands}
    with ThreadPoolExecutor(max_workers=workers) as pool:
        futures = {
            band: pool.submit(download_file, url, out_dir=base,
                              overwrite=overwrite, progress=progress)
            for band, url in urls.items()
        }
        return {band: fut.result() for band, fut in futures.items()}

Landsat 8/9

Collection-2 Level-2 surface reflectance, global, via Microsoft Planetary Computer (zero-key). Band labels are normalized to their Sentinel-2 equivalents (red→B04, nir→B08, ...), so the indices and band presets work identically to Sentinel-2. Also available as composite(..., source="landsat").

earthfetch.landsat.search_landsat

search_landsat(bbox: Sequence[float], start: str, end: str, max_cloud: float = 20.0, limit: int = 50) -> list[dict]

Search Landsat 8/9 Collection-2 L2 scenes, clearest first.

Source code in src/earthfetch/landsat.py

def search_landsat(
    bbox: Sequence[float], start: str, end: str,
    max_cloud: float = 20.0, limit: int = 50,
) -> list[dict]:
    """Search Landsat 8/9 Collection-2 L2 scenes, clearest first."""
    bbox = validate_bbox(bbox)
    body = {
        "collections": [COLLECTION],
        "bbox": list(bbox),
        "datetime": f"{start}T00:00:00Z/{end}T23:59:59Z",
        "limit": min(limit, 100),
        "query": {"eo:cloud_cover": {"lt": max_cloud},
                  "platform": {"in": ["landsat-8", "landsat-9"]}},
    }
    resp = get_session().post(PC_STAC_URL, json=body, timeout=60)
    resp.raise_for_status()
    items = resp.json().get("features", [])
    items.sort(key=lambda i: i["properties"].get("eo:cloud_cover", 100))
    logger.info("Planetary Computer: %d Landsat scene(s) %s..%s", len(items),
                start, end)
    return items[:limit]

earthfetch.landsat.load_landsat

load_landsat(bbox: Sequence[float], bands: Sequence[str] = ('red', 'green', 'blue'), crs: str = 'EPSG:4326', res: float | None = None, item: dict | None = None, items: Sequence[dict] | None = None, start: str | None = None, end: str | None = None, max_cloud: float = 20.0, scale: bool = True) -> xarray.DataArray

Load Landsat 8/9 bands for a bbox as one aligned DataArray.

Accepts a single item, a list of items, or start/end dates (the clearest scene is used). Bands may be common names (red, nir, swir1), Landsat ids (B4), Sentinel ids (B04), or a preset ("true_color"). Output band labels are the Sentinel-2 equivalents so indices work unchanged. scale=True returns surface reflectance (0..1).

Source code in src/earthfetch/landsat.py

def load_landsat(
    bbox: Sequence[float],
    bands: Sequence[str] = ("red", "green", "blue"),
    crs: str = "EPSG:4326",
    res: float | None = None,
    item: dict | None = None,
    items: Sequence[dict] | None = None,
    start: str | None = None,
    end: str | None = None,
    max_cloud: float = 20.0,
    scale: bool = True,
) -> xarray.DataArray:
    """Load Landsat 8/9 bands for a bbox as one aligned DataArray.

    Accepts a single ``item``, a list of ``items``, or ``start``/``end`` dates
    (the clearest scene is used). Bands may be common names (``red``, ``nir``,
    ``swir1``), Landsat ids (``B4``), Sentinel ids (``B04``), or a preset
    (``"true_color"``). Output band labels are the Sentinel-2 equivalents so
    indices work unchanged. ``scale=True`` returns surface reflectance (0..1).
    """
    from .aoi import resolve_aoi, resolve_crs
    from .sentinel import resolve_bands

    try:
        import numpy as np

        from .load import _resolve_res, _to_dataarray, _xr
        from .raster import make_grid, warp_into_grid

        _xr()
    except ImportError as exc:
        from .exceptions import MissingDependencyError

        raise MissingDependencyError(
            "'load_landsat' needs the optional 'xarray' dependencies; "
            "install with: pip install 'earthfetch[xarray]'"
        ) from exc

    a = resolve_aoi(bbox)
    bbox = a.bbox
    crs = resolve_crs(crs, bbox)
    bands = resolve_bands(bands)
    if items is None:
        items = [item] if item is not None else [
            clearest_landsat(bbox, start, end, max_cloud=max_cloud)]
    items = list(items)

    res = _resolve_res(res, 30.0, crs)   # Landsat is 30 m
    transform, width, height = make_grid(bbox, crs, res)
    logger.info("load_landsat: %d scene(s) %s -> %dx%d @ %s",
                len(items), list(bands), width, height, crs)

    layers, labels = [], []
    for b in bands:
        layer = warp_into_grid([band_url(it, b) for it in items],
                               transform, width, height, crs)
        if scale:
            sc, off = scale_offset(items[0], b)
            layer = np.where(np.isfinite(layer),
                             np.maximum(layer * sc + off, 0.0), np.nan)
        layers.append(layer)
        labels.append(_resolve_band(b)[1])
    data = np.stack(layers)
    da = _to_dataarray(
        data, transform, width, height, crs, "landsat",
        {"scene_ids": [it["id"] for it in items],
         "datetime": items[0]["properties"].get("datetime"),
         "cloud_cover": items[0]["properties"].get("eo:cloud_cover"),
         "reflectance": scale, "source": "landsat"},
    )
    return da.assign_coords(band=("band", labels))

NAIP aerial imagery

0.6-1 m natural-color/near-infrared aerial photos, US only, via Microsoft Planetary Computer. Anonymous SAS tokens are fetched automatically.

earthfetch.naip.search_naip

search_naip(bbox: Sequence[float], year: int | None = None, limit: int = 100) -> list[dict]

NAIP STAC items covering a bbox, newest first.

year pins a specific survey year; default returns all years (callers usually keep the newest tile per footprint).

Source code in src/earthfetch/naip.py

def search_naip(
    bbox: Sequence[float],
    year: int | None = None,
    limit: int = 100,
) -> list[dict]:
    """NAIP STAC items covering a bbox, newest first.

    ``year`` pins a specific survey year; default returns all years
    (callers usually keep the newest tile per footprint).
    """
    body = {
        "collections": ["naip"],
        "bbox": list(bbox),
        "limit": min(limit, 250),
        "sortby": [{"field": "properties.datetime", "direction": "desc"}],
    }
    if year is not None:
        body["datetime"] = f"{year}-01-01T00:00:00Z/{year}-12-31T23:59:59Z"
    resp = get_session().post(PC_STAC_URL, json=body, timeout=60)
    resp.raise_for_status()
    items = resp.json().get("features", [])
    logger.info("Planetary Computer: %d NAIP item(s) for %s", len(items), tuple(bbox))
    return items