API reference¶

Array API¶

One call, bbox in, aligned xarray out. Requires earthfetch[xarray].

earthfetch.load.load_dem ¶

load_dem(bbox: Sequence[float], resolution: str = '10m', crs: str = 'EPSG:4326', res: float | None = None, source: str = 'auto') -> xarray.DataArray

Load a DEM for a bbox as an xarray.DataArray — no files, no keys.

Parameters:

Name	Type	Description	Default
`bbox`	`(min_lon, min_lat, max_lon, max_lat) in WGS84 degrees.`		required
`resolution`	`USGS dataset ("1m", "10m", "30m", "5m-ak"); ignored for`	the Copernicus source (always 30 m).	`'10m'`
`crs`	`output CRS ("EPSG:32612", "EPSG:5070", ...).`		`'EPSG:4326'`
`res`	output pixel size in ``crs`` units; defaults to the native	resolution (converted to degrees for geographic CRSs).	`None`
`source`	`"usgs", "copernicus", or "auto" (USGS first, Copernicus`	fallback outside the US).	`'auto'`

Returns:

Type	Description
`DataArray`	float32 (y, x) elevation in meters, NaN nodata, with `crs`, `transform` and `sources` attrs.

Source code in src/earthfetch/load.py

def load_dem(
    bbox: Sequence[float],
    resolution: str = "10m",
    crs: str = "EPSG:4326",
    res: float | None = None,
    source: str = "auto",
) -> xarray.DataArray:
    """Load a DEM for a bbox as an ``xarray.DataArray`` — no files, no keys.

    Parameters
    ----------
    bbox : (min_lon, min_lat, max_lon, max_lat) in WGS84 degrees.
    resolution : USGS dataset ("1m", "10m", "30m", "5m-ak"); ignored for
        the Copernicus source (always 30 m).
    crs : output CRS ("EPSG:32612", "EPSG:5070", ...).
    res : output pixel size in ``crs`` units; defaults to the native
        resolution (converted to degrees for geographic CRSs).
    source : "usgs", "copernicus", or "auto" (USGS first, Copernicus
        fallback outside the US).

    Returns
    -------
    xarray.DataArray
        float32 (y, x) elevation in meters, NaN nodata, with ``crs``,
        ``transform`` and ``sources`` attrs.
    """
    _xr()  # fail fast on missing xarray, before any network work
    a = resolve_aoi(bbox)
    bbox = a.bbox
    crs = resolve_crs(crs, bbox)
    urls: list = []
    used = source
    if source in ("usgs", "auto"):
        try:
            urls = dem_tile_urls(bbox, resolution=resolution)
            used = "usgs"
        except Exception:
            if source == "usgs":
                raise
    if not urls:
        if source == "usgs":
            raise TileNotFoundError(f"no USGS {resolution} tiles cover {bbox}")
        urls = copernicus_dem_urls(bbox)
        used = "copernicus"
        resolution = "30m"

    native = 30.0 if used == "copernicus" else _DEM_NATIVE_M[resolution]
    res = _resolve_res(res, native, crs)
    transform, width, height = make_grid(bbox, crs, res)
    logger.info("load_dem: %s %s -> %dx%d @ %s", used, resolution, width, height, crs)
    data = warp_into_grid(urls, transform, width, height, crs)
    return _to_dataarray(
        data, transform, width, height, crs, "dem",
        {"units": "m", "source": used, "resolution": resolution, "sources": urls},
    )

earthfetch.load.load_sentinel2 ¶

load_sentinel2(bbox: Sequence[float], bands: Sequence[str] = ('B04', 'B03', 'B02'), 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 Sentinel-2 L2A bands for a bbox as one aligned DataArray.

Pass one of: a single STAC item; a list of items (from covering_scenes — mosaicked together to fill an AOI that spans tile boundaries); or start/end dates, in which case the single clearest scene in the range is used.

res defaults to the finest native resolution among bands. scale=True (default) converts DNs to surface reflectance (0..1) using STAC metadata — matching composite/time_series so indices behave the same on every source; pass scale=False for raw DNs. Returns a float32 DataArray (band, y, x), NaN nodata, with scene metadata in attrs. Multi-band assets like TCI are not supported here — use download_sentinel2 for those.

Source code in src/earthfetch/load.py

def load_sentinel2(
    bbox: Sequence[float],
    bands: Sequence[str] = ("B04", "B03", "B02"),
    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 Sentinel-2 L2A bands for a bbox as one aligned DataArray.

    Pass one of: a single STAC ``item``; a list of ``items`` (from
    ``covering_scenes`` — mosaicked together to fill an AOI that spans
    tile boundaries); or ``start``/``end`` dates, in which case the single
    clearest scene in the range is used.

    ``res`` defaults to the finest native resolution among ``bands``.
    ``scale=True`` (default) converts DNs to surface reflectance (0..1)
    using STAC metadata — matching ``composite``/``time_series`` so indices
    behave the same on every source; pass ``scale=False`` for raw DNs.
    Returns a float32 DataArray (band, y, x), NaN nodata, with scene
    metadata in attrs. Multi-band assets like TCI are not supported here —
    use ``download_sentinel2`` for those.
    """
    _xr()  # fail fast on missing xarray, before any network work
    a = resolve_aoi(bbox)
    bbox = a.bbox
    crs = resolve_crs(crs, bbox)
    bands = resolve_bands(bands)
    if items is None:
        if item is None:
            if start is None or end is None:
                raise ValueError("pass item=..., items=..., or start=/end= dates")
            item = clearest_scene(bbox, start, end, max_cloud=max_cloud)
        items = [item]
    items = list(items)
    if not items:
        raise ValueError("items is empty")

    native = min(BAND_RESOLUTION.get(b.upper(), 10) for b in bands)
    res = _resolve_res(res, float(native), crs)
    transform, width, height = make_grid(bbox, crs, res)
    logger.info("load_sentinel2: %d scene(s) %s -> %dx%d @ %s",
                len(items), list(bands), width, height, crs)

    layers = []
    for b in bands:
        # mosaic this band across every scene (warp_into_grid merges a list)
        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)
    data = np.stack(layers)
    da = _to_dataarray(
        data, transform, width, height, crs, "sentinel2",
        {
            "scene_id": items[0]["id"] if len(items) == 1 else "",
            "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,
            "sources": [band_url(items[0], b) for b in bands],
        },
    )
    return da.assign_coords(band=("band", [b.upper() for b in bands]))

earthfetch.load.stack ¶

stack(bbox: Sequence[float], crs: str, res: float, bands: Sequence[str] = ('B04', 'B08'), start: str | None = None, end: str | None = None, max_cloud: float = 20.0, item: dict | None = None, dem_resolution: str = '10m', dem_source: str = 'auto', scale: bool = True) -> xarray.Dataset

DEM + Sentinel-2 bands on one pixel-aligned grid — ML-ready.

Sentinel-2 bands are surface reflectance (0..1) by default, matching the rest of the array API; pass scale=False for raw DNs. Returns an xarray.Dataset with a dem variable plus one variable per band (B04, B08, ...), all float32 on the same (y, x) grid.

Source code in src/earthfetch/load.py

def stack(
    bbox: Sequence[float],
    crs: str,
    res: float,
    bands: Sequence[str] = ("B04", "B08"),
    start: str | None = None,
    end: str | None = None,
    max_cloud: float = 20.0,
    item: dict | None = None,
    dem_resolution: str = "10m",
    dem_source: str = "auto",
    scale: bool = True,
) -> xarray.Dataset:
    """DEM + Sentinel-2 bands on one pixel-aligned grid — ML-ready.

    Sentinel-2 bands are surface reflectance (0..1) by default, matching the
    rest of the array API; pass ``scale=False`` for raw DNs. Returns an
    ``xarray.Dataset`` with a ``dem`` variable plus one variable per band
    (``B04``, ``B08``, ...), all float32 on the same (y, x) grid.
    """
    xr = _xr()
    a = resolve_aoi(bbox)
    bbox = a.bbox
    crs = resolve_crs(crs, bbox)
    dem = load_dem(bbox, resolution=dem_resolution, crs=crs, res=res,
                   source=dem_source)
    s2 = load_sentinel2(bbox, bands=bands, crs=crs, res=res, item=item,
                        start=start, end=end, max_cloud=max_cloud, scale=scale)
    ds = xr.Dataset({"dem": dem})
    for b in s2.band.values:
        ds[str(b)] = s2.sel(band=b).drop_vars("band")
    ds.attrs = {**s2.attrs, "dem_source": dem.attrs["source"], "crs": str(crs)}
    return ds

earthfetch._composite.composite ¶

composite(aoi, bands: Sequence[str] = ('B04', 'B03', 'B02'), start: str = None, end: str = None, crs: str = 'utm', res: float | None = None, method: str = 'median', mask_clouds: bool = True, max_cloud: float = 60.0, max_scenes: int = 8, scale: bool = True, clip: bool = None, source: str = 'sentinel2') -> xarray.DataArray

Cloud-free composite of Sentinel-2 or Landsat bands over a date range.

Searches every scene touching the AOI, keeps the max_scenes clearest acquisition days (all MGRS tiles of each day, so seams disappear), masks invalid pixels with the SCL layer, and reduces per pixel.

Parameters:

Name	Type	Description	Default
`aoi`	`bbox tuple, GeoJSON, .geojson path, shapely geometry, or place name.`		required
`bands`	`ESA band ids. SCL/TCI are not composable inputs here.`		`('B04', 'B03', 'B02')`
`start`	`ISO dates bounding the search.`		`None`
`end`	`ISO dates bounding the search.`		`None`
`crs`	`output CRS; "utm" (default) picks the AOI's UTM zone.`		`'utm'`
`res`	`pixel size in CRS units; defaults to the finest native band res.`		`None`
`method`	`"median" (robust, default), "mean", or "first" (first valid,`	clearest day first — fastest).	`'median'`
`mask_clouds`	`mask SCL classes {0,1,3,8,9,10} before compositing.`		`True`
`max_cloud`	`scene-level cloud prefilter for the search.`		`60.0`
`max_scenes`	`acquisition days to blend.`		`8`
`scale`	`convert DNs to surface reflectance (0..1) using STAC metadata.`		`True`
`clip`	`NaN-out pixels outside a polygon AOI. Default (None): clip`	polygons you passed explicitly, but keep the full rectangle for geocoded place names (pass clip=True to cut to a city boundary).	`None`

Returns:

Type	Description
`DataArray`	float32 (band, y, x), NaN nodata, with the scene ids and dates used in `attrs`.

Source code in src/earthfetch/_composite.py

def composite(
    aoi,
    bands: Sequence[str] = ("B04", "B03", "B02"),
    start: str = None,
    end: str = None,
    crs: str = "utm",
    res: float | None = None,
    method: str = "median",
    mask_clouds: bool = True,
    max_cloud: float = 60.0,
    max_scenes: int = 8,
    scale: bool = True,
    clip: bool = None,
    source: str = "sentinel2",
) -> xarray.DataArray:
    """Cloud-free composite of Sentinel-2 or Landsat bands over a date range.

    Searches every scene touching the AOI, keeps the ``max_scenes`` clearest
    acquisition days (all MGRS tiles of each day, so seams disappear), masks
    invalid pixels with the SCL layer, and reduces per pixel.

    Parameters
    ----------
    aoi : bbox tuple, GeoJSON, .geojson path, shapely geometry, or place name.
    bands : ESA band ids. SCL/TCI are not composable inputs here.
    start, end : ISO dates bounding the search.
    crs : output CRS; "utm" (default) picks the AOI's UTM zone.
    res : pixel size in CRS units; defaults to the finest native band res.
    method : "median" (robust, default), "mean", or "first" (first valid,
        clearest day first — fastest).
    mask_clouds : mask SCL classes {0,1,3,8,9,10} before compositing.
    max_cloud : scene-level cloud prefilter for the search.
    max_scenes : acquisition days to blend.
    scale : convert DNs to surface reflectance (0..1) using STAC metadata.
    clip : NaN-out pixels outside a polygon AOI. Default (None): clip
        polygons you passed explicitly, but keep the full rectangle for
        geocoded place names (pass clip=True to cut to a city boundary).

    Returns
    -------
    xarray.DataArray
        float32 (band, y, x), NaN nodata, with the scene ids and dates
        used in ``attrs``.
    """
    if method not in ("median", "mean", "first"):
        raise ValueError("method must be 'median', 'mean', or 'first'")
    _xr()  # fail fast on missing xarray, before any network work
    search, url_of, scale_of, validity, native_res, label_of = \
        _source_adapter(source)
    bands = resolve_bands(bands)
    a = resolve_aoi(aoi)
    crs = resolve_crs(crs, a.bbox)
    items = search(a.bbox, start, end, max_cloud=max_cloud, limit=100)
    if not items:
        raise NoScenesError(
            f"no {source} scenes for {a.bbox} in {start}..{end} "
            f"with cloud < {max_cloud}%"
        )
    groups = _select_day_groups(items, max_scenes)
    scenes = [it for g in groups for it in g]
    logger.info("composite: %d %s scene(s) over %d day(s), method=%s",
                len(scenes), source, len(groups), method)

    from .load import _resolve_res, _to_dataarray

    res = _resolve_res(res, native_res(bands), crs)
    transform, width, height = make_grid(a.bbox, crs, res)

    acc = np.full((len(bands), height, width), np.nan, dtype="float32")
    stacks: list = []
    for item in scenes:
        valid = validity(item, transform, width, height, crs) if mask_clouds else None
        layer = np.full((len(bands), height, width), np.nan, dtype="float32")
        for i, b in enumerate(bands):
            data = warp_into_grid([url_of(item, b)], transform, width,
                                  height, crs)
            if scale:
                sc, off = scale_of(item, b)
                data = np.where(np.isfinite(data),
                                np.maximum(data * sc + off, 0.0), np.nan)
            if valid is not None:
                data = np.where(valid, data, np.nan)
            layer[i] = data
        if method == "first":
            hole = np.isnan(acc)
            acc[hole] = layer[hole]
            if not np.isnan(acc).any():
                break
        else:
            stacks.append(layer)

    if method == "median" and stacks:
        acc = np.nanmedian(np.stack(stacks), axis=0)
    elif method == "mean" and stacks:
        acc = np.nanmean(np.stack(stacks), axis=0)

    if clip is None:
        clip = a.clip_default
    if clip and a.geometry is not None:
        mask_to_geometry(acc, a.geometry, transform, crs)

    da = _to_dataarray(
        acc.astype("float32"), transform, width, height, crs, "composite",
        {
            "method": method,
            "scenes": [it["id"] for it in scenes],
            "dates": sorted({it["properties"]["datetime"][:10] for it in scenes}),
            "cloud_masked": mask_clouds,
            "reflectance": scale,
            "aoi_name": a.name or "",
            "source": source,
        },
    )
    return da.assign_coords(band=("band", [label_of(b) for b in bands]))

earthfetch._terrain.terrain ¶

terrain(aoi, products=('dem', 'slope', 'aspect', 'hillshade'), resolution: str = '10m', crs: str = 'utm', res: float = None, source: str = 'auto', azimuth: float = 315.0, altitude: float = 45.0, clip: bool = None) -> xarray.Dataset

DEM + terrain derivatives for any AOI as an aligned Dataset.

products chooses among "dem" (meters), "slope" (degrees), "aspect" (compass degrees), "hillshade" (0-255). Other args follow load_dem.

Source code in src/earthfetch/_terrain.py

def terrain(
    aoi,
    products=("dem", "slope", "aspect", "hillshade"),
    resolution: str = "10m",
    crs: str = "utm",
    res: float = None,
    source: str = "auto",
    azimuth: float = 315.0,
    altitude: float = 45.0,
    clip: bool = None,
) -> xarray.Dataset:
    """DEM + terrain derivatives for any AOI as an aligned Dataset.

    ``products`` chooses among "dem" (meters), "slope" (degrees), "aspect"
    (compass degrees), "hillshade" (0-255). Other args follow ``load_dem``.
    """
    unknown = set(products) - set(PRODUCTS)
    if unknown:
        raise ValueError(f"unknown products {sorted(unknown)}; pick from {PRODUCTS}")
    from .load import _xr, load_dem

    xr = _xr()
    a = resolve_aoi(aoi)
    crs = resolve_crs(crs, a.bbox)
    dem = load_dem(a.bbox, resolution=resolution, crs=crs, res=res, source=source)
    pixel = abs(dem.attrs["transform"][0])
    logger.info("terrain: %s at %.1f-unit pixels in %s", list(products), pixel, crs)

    layers = {"dem": dem.values}
    if "slope" in products or "aspect" in products:
        layers["slope"], layers["aspect"] = slope_aspect(dem.values, pixel)
    if "hillshade" in products:
        layers["hillshade"] = hillshade(dem.values, pixel, azimuth, altitude)

    if clip is None:
        clip = a.clip_default
    if clip and a.geometry is not None:
        transform, width, height = make_grid(a.bbox, crs, pixel)
        for arr in layers.values():
            mask_to_geometry(arr, a.geometry, transform, crs)

    ds = xr.Dataset(
        {name: dem.copy(data=layers[name]).rename(name)
         for name in products},
        attrs={**dem.attrs, "products": list(products)},
    )
    return ds

earthfetch.naip.load_naip ¶

load_naip(aoi, bands: Sequence[str] = ('R', 'G', 'B'), crs: str = 'utm', res: float | None = None, year: int | None = None, clip: bool | None = None) -> xarray.DataArray

US aerial imagery for any AOI as an xarray.DataArray.

Parameters:

Name	Type	Description	Default
`aoi`	`bbox, GeoJSON, .geojson path, shapely geometry, or place name.`		required
`bands`	`subset of R, G, B, N (near-infrared).`		`('R', 'G', 'B')`
`crs`	`output CRS; "utm" picks the AOI's zone.`		`'utm'`
`res`	`pixel size in CRS units; defaults to 1 m (native is 0.6-1 m).`		`None`
`year`	`pin a survey year; default mosaics the newest available.`		`None`
`clip`	`NaN-out pixels outside a polygon AOI. Defaults to True for`	polygons you pass explicitly, False for geocoded place names.	`None`

Returns:

Type	Description
`DataArray`	float32 (band, y, x) of 0-255 values, NaN nodata.

Source code in src/earthfetch/naip.py

def load_naip(
    aoi,
    bands: Sequence[str] = ("R", "G", "B"),
    crs: str = "utm",
    res: float | None = None,
    year: int | None = None,
    clip: bool | None = None,
) -> xarray.DataArray:
    """US aerial imagery for any AOI as an ``xarray.DataArray``.

    Parameters
    ----------
    aoi : bbox, GeoJSON, .geojson path, shapely geometry, or place name.
    bands : subset of R, G, B, N (near-infrared).
    crs : output CRS; "utm" picks the AOI's zone.
    res : pixel size in CRS units; defaults to 1 m (native is 0.6-1 m).
    year : pin a survey year; default mosaics the newest available.
    clip : NaN-out pixels outside a polygon AOI. Defaults to True for
        polygons you pass explicitly, False for geocoded place names.

    Returns
    -------
    xarray.DataArray
        float32 (band, y, x) of 0-255 values, NaN nodata.
    """
    from .aoi import resolve_aoi, resolve_crs

    # heavy deps are imported here (not at module top) so the requests-only
    # search_naip/naip_tile_urls stay importable on a core install; guard
    # them so a missing dep names the right extra and fails before any network
    try:
        import numpy as np

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

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

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

    a = resolve_aoi(aoi)
    crs = resolve_crs(crs, a.bbox)
    urls = naip_tile_urls(a.bbox, year=year)
    res = _resolve_res(res, 1.0, crs)
    transform, width, height = make_grid(a.bbox, crs, res)
    logger.info("load_naip: %d tile(s) -> %dx%d @ %s", len(urls), width, height, crs)

    layers = [
        warp_into_grid(urls, transform, width, height, crs,
                       band=NAIP_BANDS[b.upper()])
        for b in bands
    ]
    data = np.stack(layers)
    if clip is None:
        clip = a.clip_default
    if clip and a.geometry is not None:
        mask_to_geometry(data, a.geometry, transform, crs)
    da = _to_dataarray(
        data, transform, width, height, crs, "naip",
        {"source": "naip", "year": year or "newest", "n_tiles": len(urls)},
    )
    return da.assign_coords(band=("band", [b.upper() for b in bands]))

earthfetch.timeseries.time_series ¶

time_series(aoi, bands: Sequence[str] = ('B04', 'B03', 'B02'), start: str = None, end: str = None, crs: str = 'utm', res: float | None = None, max_cloud: float = 60.0, mask_clouds: bool = True, scale: bool = True, max_steps: int = 50, freq: str | None = None, clip: bool = None) -> xarray.DataArray

Cloud-masked Sentinel-2 time series over a date range.

Each clear acquisition day becomes a time step; the day's MGRS tiles are mosaicked together (first valid pixel wins), cloud/shadow pixels are masked with SCL, and values are reflectance-scaled.

Parameters:

Name	Type	Description	Default
`aoi`	`bbox tuple, GeoJSON, .geojson path, shapely geometry, or place name.`		required
`bands`	`ESA band ids. SCL/TCI are not composable inputs here.`		`('B04', 'B03', 'B02')`
`start`	`ISO dates bounding the search.`		`None`
`end`	`ISO dates bounding the search.`		`None`
`crs`	`output CRS; "utm" (default) picks the AOI's UTM zone.`		`'utm'`
`res`	`pixel size in CRS units; defaults to the finest native band res.`		`None`
`max_cloud`	`scene-level cloud prefilter for the search.`		`60.0`
`mask_clouds`	`mask SCL classes {0,1,3,8,9,10} before stacking.`		`True`
`scale`	`convert DNs to surface reflectance (0..1) using STAC metadata.`		`True`
`max_steps`	`cap on acquisition days (time steps) returned.`		`50`
`freq`	`optional pandas offset alias ("MS", "W", "QS", ...) to resample`	the time axis with a per-period median (e.g. monthly composites).	`None`
`clip`	NaN-out pixels outside a polygon AOI (see ``composite``).		`None`

Returns:

Type	Description
`DataArray`	float32 (time, band, y, x), NaN nodata, `time` as datetime64 and `band` as ESA ids, with `crs`/`transform` in `attrs`.

Source code in src/earthfetch/timeseries.py

def time_series(
    aoi,
    bands: Sequence[str] = ("B04", "B03", "B02"),
    start: str = None,
    end: str = None,
    crs: str = "utm",
    res: float | None = None,
    max_cloud: float = 60.0,
    mask_clouds: bool = True,
    scale: bool = True,
    max_steps: int = 50,
    freq: str | None = None,
    clip: bool = None,
) -> xarray.DataArray:
    """Cloud-masked Sentinel-2 time series over a date range.

    Each clear acquisition day becomes a time step; the day's MGRS tiles are
    mosaicked together (first valid pixel wins), cloud/shadow pixels are
    masked with SCL, and values are reflectance-scaled.

    Parameters
    ----------
    aoi : bbox tuple, GeoJSON, .geojson path, shapely geometry, or place name.
    bands : ESA band ids. SCL/TCI are not composable inputs here.
    start, end : ISO dates bounding the search.
    crs : output CRS; "utm" (default) picks the AOI's UTM zone.
    res : pixel size in CRS units; defaults to the finest native band res.
    max_cloud : scene-level cloud prefilter for the search.
    mask_clouds : mask SCL classes {0,1,3,8,9,10} before stacking.
    scale : convert DNs to surface reflectance (0..1) using STAC metadata.
    max_steps : cap on acquisition days (time steps) returned.
    freq : optional pandas offset alias ("MS", "W", "QS", ...) to resample
        the time axis with a per-period median (e.g. monthly composites).
    clip : NaN-out pixels outside a polygon AOI (see ``composite``).

    Returns
    -------
    xarray.DataArray
        float32 (time, band, y, x), NaN nodata, ``time`` as datetime64 and
        ``band`` as ESA ids, with ``crs``/``transform`` in ``attrs``.
    """
    from .load import _coords, _resolve_res, _xr

    xr = _xr()
    bands = resolve_bands(bands)
    a = resolve_aoi(aoi)
    crs = resolve_crs(crs, a.bbox)
    items = search_sentinel2(a.bbox, start, end, max_cloud=max_cloud, limit=100)
    if not items:
        raise NoScenesError(
            f"no scenes for {a.bbox} in {start}..{end} with cloud < {max_cloud}%"
        )
    groups = _select_day_groups(items, max_steps)
    groups.sort(key=lambda g: g[0]["properties"]["datetime"][:10])  # time order
    logger.info("time_series: %d day(s), bands=%s", len(groups), list(bands))

    native = min(BAND_RESOLUTION.get(b.upper(), 10) for b in bands)
    res = _resolve_res(res, float(native), crs)
    transform, width, height = make_grid(a.bbox, crs, res)
    do_clip = a.clip_default if clip is None else clip

    times, slices = [], []
    for g in groups:
        day = g[0]["properties"]["datetime"][:10]
        layer = np.full((len(bands), height, width), np.nan, dtype="float32")
        for scene in g:  # fill holes across the day's tiles
            s = _render_scene(scene, bands, transform, width, height, crs,
                              mask_clouds, scale)
            hole = np.isnan(layer)
            layer[hole] = s[hole]
        if do_clip and a.geometry is not None:
            mask_to_geometry(layer, a.geometry, transform, crs)
        times.append(np.datetime64(day))
        slices.append(layer)

    xs, ys = _coords(transform, width, height)
    da = xr.DataArray(
        np.stack(slices),
        dims=("time", "band", "y", "x"),
        coords={"time": times, "band": [b.upper() for b in bands],
                "y": ys, "x": xs},
        name="sentinel2",
        attrs={
            "crs": str(crs),
            "transform": tuple(transform)[:6],
            "reflectance": scale,
            "cloud_masked": mask_clouds,
            "aoi_name": a.name or "",
        },
    )
    if freq:
        da = da.resample(time=freq).median()
        da.attrs.update({"crs": str(crs), "transform": tuple(transform)[:6]})
    return da

earthfetch.load.elevation ¶

elevation(points, source: str = 'auto', resolution: str = '10m', with_source: bool = False)

Elevation in meters at one or more (lon, lat) points.

The most direct question a DEM answers. Elevation is mixed-source: with source="auto" (default) it uses USGS 3DEP inside the United States and falls back to Copernicus GLO-30 everywhere else, so it works worldwide. Pass with_source=True to also get which source was used.

Parameters:

Name	Type	Description	Default
`points`	a single ``(lon, lat)``, a list of them, or GeoJSON	Point/MultiPoint (WGS84).	required
`source`	`"usgs", "copernicus", or "auto".`		`'auto'`
`resolution`	`DEM resolution when the USGS source is used.`		`'10m'`
`with_source`	`also return the DEM source ("usgs" or "copernicus").`		`False`

Returns:

Type	Description
`float or ndarray`	A float for a single point, else an array of elevations (NaN where outside coverage). If `with_source` is True, returns `(elevation, source)`. Best for nearby points — far-apart points load a DEM spanning their whole bounding box.

Source code in src/earthfetch/load.py

def elevation(points, source: str = "auto", resolution: str = "10m",
              with_source: bool = False):
    """Elevation in meters at one or more ``(lon, lat)`` points.

    The most direct question a DEM answers. Elevation is **mixed-source**:
    with ``source="auto"`` (default) it uses USGS 3DEP inside the United
    States and falls back to Copernicus GLO-30 everywhere else, so it works
    worldwide. Pass ``with_source=True`` to also get which source was used.

    Parameters
    ----------
    points : a single ``(lon, lat)``, a list of them, or GeoJSON
        Point/MultiPoint (WGS84).
    source : "usgs", "copernicus", or "auto".
    resolution : DEM resolution when the USGS source is used.
    with_source : also return the DEM source ("usgs" or "copernicus").

    Returns
    -------
    float or numpy.ndarray
        A float for a single point, else an array of elevations (NaN where
        outside coverage). If ``with_source`` is True, returns
        ``(elevation, source)``. Best for nearby points — far-apart points
        load a DEM spanning their whole bounding box.
    """
    from .zonal import _as_points, sample

    pts = _as_points(points)
    lons = [p[0] for p in pts]
    lats = [p[1] for p in pts]
    pad = 0.002
    bbox = (min(lons) - pad, min(lats) - pad, max(lons) + pad, max(lats) + pad)
    dem = load_dem(bbox, resolution=resolution, crs="utm", source=source)
    vals = sample(dem, pts)
    result = float(vals[0]) if len(pts) == 1 else vals
    return (result, dem.attrs["source"]) if with_source else result

Analysis¶

Extract values from any earthfetch raster at points or over polygons. Requires the raster extra.

earthfetch.zonal.sample ¶

sample(obj, points)

Read raster values at points (nearest pixel).

Parameters:

Name	Type	Description	Default
`obj`	an earthfetch DataArray (2D ``(y, x)`` or 3D ``(band, y, x)``).		required
`points`	a single ``(lon, lat)``, a list of them, or GeoJSON	Point/MultiPoint/Feature/FeatureCollection (WGS84).	required

Returns:

Type	Description
`ndarray`	Shape `(n_points,)` for a 2D input, or `(n_bands, n_points)` for a banded input. Points outside the raster are NaN.

Source code in src/earthfetch/zonal.py

def sample(obj, points):
    """Read raster values at points (nearest pixel).

    Parameters
    ----------
    obj : an earthfetch DataArray (2D ``(y, x)`` or 3D ``(band, y, x)``).
    points : a single ``(lon, lat)``, a list of them, or GeoJSON
        Point/MultiPoint/Feature/FeatureCollection (WGS84).

    Returns
    -------
    numpy.ndarray
        Shape ``(n_points,)`` for a 2D input, or ``(n_bands, n_points)``
        for a banded input. Points outside the raster are NaN.
    """
    transform, crs = _georef(obj)
    pts = _as_points(points)
    lons = [p[0] for p in pts]
    lats = [p[1] for p in pts]
    xs, ys = _warp_xy("EPSG:4326", crs, lons, lats)
    rows, cols = rowcol(transform, xs, ys)
    rows = np.atleast_1d(np.asarray(rows, dtype=int))
    cols = np.atleast_1d(np.asarray(cols, dtype=int))
    data = np.asarray(obj.values, dtype="float32")
    h, w = data.shape[-2], data.shape[-1]
    inb = (rows >= 0) & (rows < h) & (cols >= 0) & (cols < w)

    if data.ndim == 2:
        out = np.full(len(rows), np.nan, dtype="float32")
        out[inb] = data[rows[inb], cols[inb]]
        return out
    out = np.full((data.shape[0], len(rows)), np.nan, dtype="float32")
    out[:, inb] = data[:, rows[inb], cols[inb]]
    return out

earthfetch.zonal.zonal_stats ¶

zonal_stats(obj, polygons, stats=('mean', 'min', 'max', 'std', 'count'))

Aggregate raster values within each polygon.

Parameters:

Name	Type	Description	Default
`obj`	`an earthfetch DataArray (2D or banded 3D).`		required
`polygons`	`GeoJSON Polygon/MultiPolygon/Feature/FeatureCollection, a`	shapely geometry, or a list of any of these (WGS84).	required
`stats`	`any of "mean", "min", "max", "std", "median", "sum", plus`	"count" (number of valid pixels).	`('mean', 'min', 'max', 'std', 'count')`

Returns:

Type	Description
`list of dict`	One dict per polygon. For a 2D input each dict maps stat name to value; for a banded input it maps `band -> {stat: value}`.

Source code in src/earthfetch/zonal.py

def zonal_stats(obj, polygons, stats=("mean", "min", "max", "std", "count")):
    """Aggregate raster values within each polygon.

    Parameters
    ----------
    obj : an earthfetch DataArray (2D or banded 3D).
    polygons : GeoJSON Polygon/MultiPolygon/Feature/FeatureCollection, a
        shapely geometry, or a list of any of these (WGS84).
    stats : any of "mean", "min", "max", "std", "median", "sum", plus
        "count" (number of valid pixels).

    Returns
    -------
    list of dict
        One dict per polygon. For a 2D input each dict maps stat name to
        value; for a banded input it maps ``band -> {stat: value}``.
    """
    from .raster import mask_to_geometry

    transform, crs = _georef(obj)
    geoms = _as_geometries(polygons)
    data = np.asarray(obj.values, dtype="float32")
    banded = data.ndim == 3
    names = _band_names(obj, data.shape[0]) if banded else None

    results = []
    for geom in geoms:
        masked = mask_to_geometry(data.copy(), geom, transform, crs)
        if banded:
            results.append({names[i]: _reduce(masked[i], stats)
                            for i in range(masked.shape[0])})
        else:
            results.append(_reduce(masked, stats))
    return results

Interoperability¶

Bridge earthfetch results into the rioxarray / rasterio ecosystem. Requires the interop extra (pip install earthfetch[interop]).

earthfetch.interop.to_rioxarray ¶

to_rioxarray(obj)

Attach the earthfetch CRS/transform so the .rio accessor works.

Returns the same object with its CRS written; rioxarray infers the transform from the x/y coordinates. After this you can call any .rio method::

da = ef.to_rioxarray(ef.composite("Moab, Utah", ...))
da.rio.to_raster("moab.tif")
da.rio.reproject("EPSG:4326")

Source code in src/earthfetch/interop.py

def to_rioxarray(obj):
    """Attach the earthfetch CRS/transform so the ``.rio`` accessor works.

    Returns the same object with its CRS written; rioxarray infers the
    transform from the ``x``/``y`` coordinates. After this you can call any
    ``.rio`` method::

        da = ef.to_rioxarray(ef.composite("Moab, Utah", ...))
        da.rio.to_raster("moab.tif")
        da.rio.reproject("EPSG:4326")
    """
    _rioxarray()
    return obj.rio.write_crs(_crs_of(obj))

earthfetch.interop.reproject ¶

reproject(obj, dst_crs: str, resolution: float | None = None, resampling: str = 'bilinear')

Reproject a DataArray/Dataset to dst_crs (via rioxarray).

Parameters:

Name	Type	Description	Default
`obj`	an earthfetch result (CRS taken from ``attrs``).		required
`dst_crs`	`target CRS, e.g. "EPSG:4326" or "EPSG:5070".`		required
`resolution`	output pixel size in ``dst_crs`` units; None keeps the	source resolution reprojected.	`None`
`resampling`	`"nearest", "bilinear", "cubic", "average", ... (rasterio).`		`'bilinear'`

Returns:

Type	Description
`DataArray or Dataset`	The reprojected object, with a working `.rio` accessor.

Source code in src/earthfetch/interop.py

def reproject(obj, dst_crs: str, resolution: float | None = None,
              resampling: str = "bilinear"):
    """Reproject a DataArray/Dataset to ``dst_crs`` (via rioxarray).

    Parameters
    ----------
    obj : an earthfetch result (CRS taken from ``attrs``).
    dst_crs : target CRS, e.g. "EPSG:4326" or "EPSG:5070".
    resolution : output pixel size in ``dst_crs`` units; None keeps the
        source resolution reprojected.
    resampling : "nearest", "bilinear", "cubic", "average", ... (rasterio).

    Returns
    -------
    xarray.DataArray or xarray.Dataset
        The reprojected object, with a working ``.rio`` accessor.
    """
    _rioxarray()
    from rasterio.enums import Resampling

    how = getattr(Resampling, resampling)
    src = obj.rio.write_crs(_crs_of(obj))
    return src.rio.reproject(dst_crs, resolution=resolution, resampling=how)

Spectral indices¶

Each accepts an xarray.Dataset (band-named variables) or a DataArray with a band coordinate, and returns a named DataArray.

earthfetch.indices.normalized_difference ¶

normalized_difference(obj, band_a: str, band_b: str, name: str = 'nd')

Generic normalized difference (A - B)/(A + B) for any two bands.

Escape hatch beyond the twelve named indices, e.g. normalized_difference(ds, "B03", "B08", name="ndwi").