SoilGrids tools

`cycles.soilgrids`

`SoilGrids(path, *, maps=ALL_MAPS, crs=None, aggregated=None)`

Read SoilGrids rasters and build Cycles-compatible soil profiles.

Load SoilGrids rasters for selected maps.

Parameters:	`path` (`str \| Path`) – Directory containing downloaded SoilGrids rasters. `maps` (`list[str]`, default: `ALL_MAPS` ) – Map identifiers in `property@layer` format. `crs` (`str \| None`, default: `None` ) – Optional target CRS for reprojection. `aggregated` (`int \| None`, default: `None` ) – Optional aggregated resolution (1000 or 5000 m).

Source code in cycles/soilgrids/soilgrids.py

def __init__(self, path: str | Path, *, maps: list[str]=ALL_MAPS, crs: str | None=None, aggregated: int | None=None) -> None:
    """Load SoilGrids rasters for selected maps.

    Args:
        path: Directory containing downloaded SoilGrids rasters.
        maps: Map identifiers in ``property@layer`` format.
        crs: Optional target CRS for reprojection.
        aggregated: Optional aggregated resolution (1000 or 5000 m).
    """
    if aggregated is not None and aggregated not in [1000, 5000]:
        raise ValueError(f'Invalid value for aggregated: {aggregated}. Supported values are 1000, and 5000.')
    self.crs: str = crs if crs is not None else HOMOLOSINE
    self.maps: dict[str, xarray.DataArray] = _read_maps(Path(path), maps, crs, aggregated)
    self.transformer = Transformer.from_crs('epsg:4326', self.crs, always_xy=True)
    self.matched_maps: pd.DataFrame | None = None

`reproject_match(*, reference_xds, reference_name, boundary)`

Reproject loaded maps to a reference raster grid.

Parameters:	`reference_xds` (`DataArray`) – Reference raster whose CRS, transform, and resolution are used for reprojection. `reference_name` (`str`) – Column name used for the reference raster values in the generated matched table. `boundary` (`GeoDataFrame`) – Boundary geometry used to clip all rasters after reprojection.

Returns:	`None` – None.

Source code in cycles/soilgrids/soilgrids.py

def reproject_match(self, *, reference_xds: xarray.DataArray, reference_name: str, boundary: gpd.GeoDataFrame) -> None:
    """Reproject loaded maps to a reference raster grid.

    Args:
        reference_xds: Reference raster whose CRS, transform, and resolution are used for reprojection.
        reference_name: Column name used for the reference raster values in the generated matched table.
        boundary: Boundary geometry used to clip all rasters after reprojection.

    Returns:
        None.
    """
    reference_xds = reference_xds.rio.clip([boundary], from_disk=True)
    df = pd.DataFrame(reference_xds[0].to_series().rename(reference_name))

    for n, m in self.maps.items():
        reprojected = m.rio.reproject_match(reference_xds, resampling=Resampling.nearest)
        reprojected = reprojected.rio.clip([boundary], from_disk=True)
        multiplier = SOILGRIDS_PROPERTIES[n.split('@')[0]].multiplier
        soil_df = pd.DataFrame(reprojected[0].to_series().rename(n)) * multiplier
        df = pd.concat([df, soil_df], axis=1)

    self.matched_maps = df

`get_soil_profile(lat_lon)`

Sample loaded maps at a coordinate and build a soil profile.

Parameters:	`lat_lon` (`LatLon`) – Latitude and longitude pair in EPSG:4326.

Returns:	`list[SoilLayer]` – Soil layers populated from SoilGrids values at the nearest grid cell.

Source code in cycles/soilgrids/soilgrids.py

def get_soil_profile(self, lat_lon: LatLon) -> list[SoilLayer]:
    """Sample loaded maps at a coordinate and build a soil profile.

    Args:
        lat_lon: Latitude and longitude pair in EPSG:4326.

    Returns:
        Soil layers populated from SoilGrids values at the nearest grid cell.
    """
    values = self._extract_values(lat_lon)

    return [SoilLayer(
        top=layer.top,
        bottom=layer.bottom,
        **{p: values[f'{p}@{key}'] for p in MAPPABLE_PARAMETERS},
    ) for key, layer in SOILGRIDS_LAYERS.items()]

`generate_soil_file(fn, lat_lon, *, desc=None, hsg='', slope=None)`

Generate a Cycles soil file from SoilGrids values.

Parameters:

fn (Path | str) –

Output soil file path.
lat_lon (LatLon) –

Latitude and longitude pair in EPSG:4326.
desc (str | None, default: None ) –

Optional custom header text for the output file.
hsg (str, default: '' ) –

Optional hydrologic soil group code used for curve-number mapping.
slope (float | None, default: None ) –

Optional slope value written to the soil file header.

Returns:	`None` – None.

Source code in cycles/soilgrids/soilgrids.py

def generate_soil_file(self, fn: Path | str, lat_lon: LatLon, *, desc: str | None=None, hsg: str='', slope: float | None=None) -> None:
    """Generate a Cycles soil file from SoilGrids values.

    Args:
        fn: Output soil file path.
        lat_lon: Latitude and longitude pair in EPSG:4326.
        desc: Optional custom header text for the output file.
        hsg: Optional hydrologic soil group code used for curve-number mapping.
        slope: Optional slope value written to the soil file header.

    Returns:
        None.
    """
    profile: list[SoilLayer] = self.get_soil_profile(lat_lon)
    desc = desc if desc is not None else _build_desc(lat_lon, hsg)

    _generate_soil_file(Path(fn), profile, desc=desc, hsg=hsg, slope=slope)

`download_soilgrids_data(path, *, maps=ALL_MAPS, boundary=None, bbox=None, crs='epsg:4326')`

Download SoilGrids raster layers via WCS.

You can provide either a boundary polygon or an explicit bounding box. Bounding boxes are expected in (west, south, east, north) order.

Parameters:

path (str | Path) –

Directory where downloaded GeoTIFF files are written.
maps (list[str], default: ALL_MAPS ) –

Map identifiers in property@layer format.
boundary (Polygon | None, default: None ) –

Optional polygon used to derive a buffered bounding box.
bbox (tuple[float, float, float, float] | None, default: None ) –

Optional explicit bounding box as (west, south, east, north).
crs (str, default: 'epsg:4326' ) –

CRS of the provided bbox coordinates.

Returns:	`None` – None.

Raises:	`ValueError` – If neither `boundary` nor `bbox` is provided.

Source code in cycles/soilgrids/soilgrids.py

def download_soilgrids_data(path: str | Path, *,
    maps: list[str]=ALL_MAPS, boundary: Polygon | None=None, bbox: tuple[float, float, float, float] | None=None, crs: str='epsg:4326') -> None:
    """Download SoilGrids raster layers via WCS.

    You can provide either a boundary polygon or an explicit bounding box.  Bounding boxes are expected in
    ``(west, south, east, north)`` order.

    Args:
        path: Directory where downloaded GeoTIFF files are written.
        maps: Map identifiers in ``property@layer`` format.
        boundary: Optional polygon used to derive a buffered bounding box.
        bbox: Optional explicit bounding box as ``(west, south, east, north)``.
        crs: CRS of the provided ``bbox`` coordinates.

    Returns:
        None.

    Raises:
        ValueError: If neither ``boundary`` nor ``bbox`` is provided.
    """
    if boundary is not None and bbox is None:
        bbox = _bbox_from_boundary(boundary)
    if bbox is None:
        raise ValueError("Either boundary or bbox must be provided.")

    bbox = _get_bounding_box(bbox, crs)

    for m in maps:
        parameter, layer = m.split('@')
        v = SOILGRIDS_PROPERTIES[parameter].soilgrids_name
        wcs = WebCoverageService(f'http://maps.isric.org/mapserv?map=/map/{v}.map', version='1.0.0')

        while True:
            try:
                response = wcs.getCoverage(     # type: ignore
                    identifier=f'{v}_{layer}_mean',
                    crs='urn:ogc:def:crs:EPSG::152160',
                    bbox=bbox,
                    resx=250,
                    resy=250,
                    format='GEOTIFF_INT16',
                )
                (Path(path) / f'{v}_{layer}.tif').write_bytes(response.read())
                break
            except Exception:
                continue