Python API

The public API is exposed at the top level of the osm_rasterizer package:

from osm_rasterizer import rasterize, RasterizeResult, fetch_features, get_utm_crs

rasterize

RasterizeResult dataclass

Result of a rasterization operation.

When single_layer=False (default): array shape is (n_bands, H, W), values are 0/1 per-feature presence masks.

When single_layer=True: array shape is (1, H, W), values are 1-based category indices into categories (0 = no data). The priority order is determined by the feature list order passed to rasterize(): the last feature has the highest priority and wins conflicts.

Source code in osm_rasterizer/rasterize.py

@dataclasses.dataclass
class RasterizeResult:
    """Result of a rasterization operation.

    When ``single_layer=False`` (default): ``array`` shape is ``(n_bands, H, W)``,
    values are 0/1 per-feature presence masks.

    When ``single_layer=True``: ``array`` shape is ``(1, H, W)``, values are
    1-based category indices into ``categories`` (0 = no data).  The priority
    order is determined by the feature list order passed to ``rasterize()``:
    the **last** feature has the highest priority and wins conflicts.
    """

    array: np.ndarray      # shape: (n_bands, height, width), dtype uint8
    transform: Affine
    crs: rasterio.CRS
    band_names: list[str]
    nodata: int = 0
    categories: list[str] | None = None  # populated in single_layer mode

rasterize(bbox, features, resolution=10.0, single_layer=False, fill_nodata=False, fill_nodata_distance=None, output_path=None, transform=None, crs=None, date=None, provider='osm')

Rasterize OSM features into a GeoTIFF or return a RasterizeResult.

Parameters:

Name	Type	Description	Default
bbox	tuple[float, float, float, float]	(minx, miny, maxx, maxy) in WGS84 degrees.	required
features	list	List of OSM tag dicts, (name, source) tuples, or (name, source, options) tuples. source is either an OSM tag dict (fetched from the provider) or a pre-fetched GeoDataFrame used directly — handy for fetching a tag class once and splitting it into several bands, or for any filtering pandas can express. A GeoDataFrame source must be named (bare GeoDataFrames are rejected). The options dict controls per-feature behaviour: filter (dict or callable): keep only some of the fetched rows, enabling AND-on-attribute splits that a tag query cannot express (osmnx ORs tag keys). A dict {column: value | [values]} keeps a row when every column matches (AND); a cell matches when its value, split on ; (so multi-values like "soccer;basketball" work), intersects the allowed set — an absent column or missing cell matches nothing. A callable gdf -> gdf may instead return a filtered GeoDataFrame or a boolean row mask. Example: split leisure=pitch by surface with {"filter": {"surface": ["grass"], "sport": ["soccer"]}}. line_width (float): real-world width in metres; lines are buffered by line_width / 2 in the projected CRS before rasterization. width_from_tags (bool): derive a per-geometry width from the feature's own OSM tags — width in metres, else lanes × 3.5 m — falling back to line_width (if given) when neither tag is present or parseable. Without options, lines burn as 1-pixel-wide traces. Polygons and points are never buffered.	required
resolution	float	Pixel size in metres (ignored when transform is supplied).	10.0
single_layer	bool	If True, merge all features into a single categorical band. Pixel values are 1-based indices into the feature list (0 = no data). Conflicts are resolved by priority: features listed later win. Reorder your feature list from least to most important category.	False
fill_nodata	bool	If True, replace empty (zero) pixels with the value of their nearest non-zero neighbour. Applied per-band in multi-band mode and on the merged array in single-layer mode.	False
fill_nodata_distance	Union[float, None]	Maximum distance in pixels within which an empty pixel will be filled. Pixels farther than this from any labelled pixel are left as 0, preventing border areas outside OSM coverage from being flooded with a nearby label. Ignored when fill_nodata is False. None (default) fills all empty pixels regardless of distance.	None
output_path	Union[str, Path, None]	Write a GeoTIFF here; return None instead of RasterizeResult.	None
transform	Union[Affine, None]	Explicit affine transform (overrides resolution).	None
crs	Union[CRS, str, None]	Output CRS; auto-detected from bbox if None.	None
date	Union[str, None]	Optional ISO 8601 date string (e.g. "2020-01-01"). With provider="osm" this queries OSM data as it existed at that point in time; with provider="ohm" it selects features that existed in the real world at that date via their start_date/end_date tags (partial dates like "1900" and BCE years like "-0500" are supported).	None
provider	str	Data provider: "osm" (OpenStreetMap, default) or "ohm" (OpenHistoricalMap, CC0-licensed historical data).	'osm'

Returns:

Type	Description
RasterizeResult or None	None when output_path is given; RasterizeResult otherwise.

Source code in osm_rasterizer/rasterize.py

def rasterize(
    bbox: tuple[float, float, float, float],
    features: list,
    resolution: float = 10.0,
    single_layer: bool = False,
    fill_nodata: bool = False,
    fill_nodata_distance: Union[float, None] = None,
    output_path: Union[str, Path, None] = None,
    transform: Union[Affine, None] = None,
    crs: Union[rasterio.CRS, str, None] = None,
    date: Union[str, None] = None,
    provider: str = "osm",
) -> Union[RasterizeResult, None]:
    """Rasterize OSM features into a GeoTIFF or return a RasterizeResult.

    Parameters
    ----------
    bbox:
        ``(minx, miny, maxx, maxy)`` in WGS84 degrees.
    features:
        List of OSM tag dicts, ``(name, source)`` tuples, or
        ``(name, source, options)`` tuples.  ``source`` is either an OSM tag
        dict (fetched from the provider) or a **pre-fetched GeoDataFrame**
        used directly — handy for fetching a tag class once and splitting it
        into several bands, or for any filtering pandas can express.  A
        GeoDataFrame source must be named (bare GeoDataFrames are rejected).

        The options dict controls per-feature behaviour:

        - ``filter`` (dict or callable): keep only some of the fetched rows,
          enabling AND-on-attribute splits that a tag query cannot express
          (osmnx ORs tag keys).  A dict ``{column: value | [values]}`` keeps a
          row when *every* column matches (AND); a cell matches when its value,
          split on ``;`` (so multi-values like ``"soccer;basketball"`` work),
          intersects the allowed set — an absent column or missing cell matches
          nothing.  A callable ``gdf -> gdf`` may instead return a filtered
          GeoDataFrame or a boolean row mask.  Example: split ``leisure=pitch``
          by surface with ``{"filter": {"surface": ["grass"], "sport":
          ["soccer"]}}``.
        - ``line_width`` (float): real-world width in metres; lines are
          buffered by ``line_width / 2`` in the projected CRS before
          rasterization.
        - ``width_from_tags`` (bool): derive a per-geometry width from the
          feature's own OSM tags — ``width`` in metres, else ``lanes`` ×
          3.5 m — falling back to ``line_width`` (if given) when neither
          tag is present or parseable.

        Without options, lines burn as 1-pixel-wide traces.  Polygons and
        points are never buffered.
    resolution:
        Pixel size in metres (ignored when *transform* is supplied).
    single_layer:
        If True, merge all features into a single categorical band.  Pixel
        values are 1-based indices into the feature list (0 = no data).
        Conflicts are resolved by priority: features listed **later** win.
        Reorder your feature list from least to most important category.
    fill_nodata:
        If True, replace empty (zero) pixels with the value of their nearest
        non-zero neighbour.  Applied per-band in multi-band mode and on the
        merged array in single-layer mode.
    fill_nodata_distance:
        Maximum distance in pixels within which an empty pixel will be
        filled.  Pixels farther than this from any labelled pixel are left
        as 0, preventing border areas outside OSM coverage from being
        flooded with a nearby label.  Ignored when *fill_nodata* is False.
        ``None`` (default) fills all empty pixels regardless of distance.
    output_path:
        Write a GeoTIFF here; return None instead of RasterizeResult.
    transform:
        Explicit affine transform (overrides *resolution*).
    crs:
        Output CRS; auto-detected from *bbox* if None.
    date:
        Optional ISO 8601 date string (e.g. ``"2020-01-01"``).  With
        ``provider="osm"`` this queries OSM data as it existed at that point
        in time; with ``provider="ohm"`` it selects features that existed in
        the real world at that date via their ``start_date``/``end_date``
        tags (partial dates like ``"1900"`` and BCE years like ``"-0500"``
        are supported).
    provider:
        Data provider: ``"osm"`` (OpenStreetMap, default) or ``"ohm"``
        (OpenHistoricalMap, CC0-licensed historical data).

    Returns
    -------
    RasterizeResult or None
        None when *output_path* is given; RasterizeResult otherwise.
    """
    # 1. Validate
    minx, miny, maxx, maxy = bbox
    if minx >= maxx or miny >= maxy:
        raise ValueError(
            f"Invalid bbox {bbox}: minx must be < maxx and miny must be < maxy"
        )
    if not features:
        raise ValueError("features list must not be empty")

    # 2. Normalize features
    named_features = _normalize_features(features)

    # 3. Determine output CRS
    if crs is None:
        out_crs_pyproj = get_utm_crs(bbox)
        out_crs = rasterio.CRS.from_wkt(out_crs_pyproj.to_wkt())
    elif isinstance(crs, str):
        out_crs = rasterio.CRS.from_string(crs)
    else:
        out_crs = crs

    # 4. Reproject bbox corners to UTM
    transformer = Transformer.from_crs("EPSG:4326", out_crs.to_wkt(), always_xy=True)
    corners_x, corners_y = transformer.transform(
        [minx, maxx, minx, maxx],
        [miny, miny, maxy, maxy],
    )
    utm_minx = min(corners_x)
    utm_maxx = max(corners_x)
    utm_miny = min(corners_y)
    utm_maxy = max(corners_y)

    # 5. Compute affine transform and grid size
    if transform is None:
        out_transform = Affine.translation(utm_minx, utm_maxy) * Affine.scale(resolution, -resolution)
        width = ceil((utm_maxx - utm_minx) / resolution)
        height = ceil((utm_maxy - utm_miny) / resolution)
    else:
        out_transform = transform
        px_width = abs(transform.a)
        px_height = abs(transform.e)
        width = ceil((utm_maxx - utm_minx) / px_width)
        height = ceil((utm_maxy - utm_miny) / px_height)

    # 6. Per-feature rasterization
    bands: list[np.ndarray] = []
    band_names: list[str] = []

    for name, source, options in named_features:
        if isinstance(source, gpd.GeoDataFrame):
            gdf = source
        else:
            gdf = fetch_features(bbox, source, date=date, provider=provider)

        if "filter" in options and not gdf.empty:
            gdf = _apply_filter(gdf, options["filter"])

        if gdf.empty:
            spec_desc = "pre-fetched GeoDataFrame" if isinstance(source, gpd.GeoDataFrame) else repr(source)
            filter_desc = " after filter" if "filter" in options else ""
            warnings.warn(
                f"No features found for feature spec {spec_desc}{filter_desc} "
                f"(band '{name}'); writing zero band.",
                stacklevel=2,
            )
            bands.append(np.zeros((height, width), dtype=np.uint8))
            band_names.append(name)
            continue

        gdf_utm = gdf.to_crs(out_crs.to_wkt())
        geometries = _apply_line_widths(gdf_utm, options)

        shapes = (
            (geom, 1)
            for geom in geometries
            if geom is not None and not geom.is_empty
        )
        burned = rio_rasterize(
            shapes=shapes,
            out_shape=(height, width),
            transform=out_transform,
            fill=0,
            dtype=np.uint8,
        )
        bands.append(burned)
        band_names.append(name)

    # 7. Handle single_layer
    if single_layer:
        # Priority-based merge: later features overwrite earlier ones.
        # Pixel value = 1-based category index; 0 = no data.
        merged = np.zeros((height, width), dtype=np.uint8)
        for i, band in enumerate(bands):
            merged = np.where(band > 0, np.uint8(i + 1), merged)
        if fill_nodata:
            merged = _fill_nodata_consensus(merged, max_distance=fill_nodata_distance)
        array = merged[np.newaxis, ...]
        final_band_names = ["landcover"]
        categories = band_names
    else:
        if fill_nodata:
            bands = [_fill_nodata_consensus(b, max_distance=fill_nodata_distance) for b in bands]
        array = np.stack(bands, axis=0)
        final_band_names = band_names
        categories = None

    result = RasterizeResult(
        array=array,
        transform=out_transform,
        crs=out_crs,
        band_names=final_band_names,
        categories=categories,
    )

    # 8. Write or return
    if output_path is None:
        return result

    output_path = Path(output_path)
    n_bands, h, w = array.shape
    with rasterio.open(
        output_path,
        "w",
        driver="GTiff",
        height=h,
        width=w,
        count=n_bands,
        dtype=np.uint8,
        crs=out_crs,
        transform=out_transform,
        nodata=0,
        compress="lzw",
        tiled=True,
        blockxsize=256,
        blockysize=256,
    ) as dst:
        dst.write(array)
        dst.update_tags(BAND_NAMES=",".join(final_band_names))
        if categories is not None:
            # single_layer mode: record the category→value mapping.
            # Value 1 = categories[0], value 2 = categories[1], etc.
            dst.update_tags(CATEGORIES=",".join(categories))
        for i, band_name in enumerate(final_band_names, start=1):
            dst.update_tags(i, name=band_name)

    return None

RasterizeResult

Result of a rasterization operation.

When single_layer=False (default): array shape is (n_bands, H, W), values are 0/1 per-feature presence masks.

When single_layer=True: array shape is (1, H, W), values are 1-based category indices into categories (0 = no data). The priority order is determined by the feature list order passed to rasterize(): the last feature has the highest priority and wins conflicts.

Source code in osm_rasterizer/rasterize.py

@dataclasses.dataclass
class RasterizeResult:
    """Result of a rasterization operation.

    When ``single_layer=False`` (default): ``array`` shape is ``(n_bands, H, W)``,
    values are 0/1 per-feature presence masks.

    When ``single_layer=True``: ``array`` shape is ``(1, H, W)``, values are
    1-based category indices into ``categories`` (0 = no data).  The priority
    order is determined by the feature list order passed to ``rasterize()``:
    the **last** feature has the highest priority and wins conflicts.
    """

    array: np.ndarray      # shape: (n_bands, height, width), dtype uint8
    transform: Affine
    crs: rasterio.CRS
    band_names: list[str]
    nodata: int = 0
    categories: list[str] | None = None  # populated in single_layer mode

fetch_features

Fetch OSM or OpenHistoricalMap features for a WGS84 bounding box.

Parameters:

Name	Type	Description	Default
bbox	tuple[float, float, float, float]	(minx, miny, maxx, maxy) in WGS84 degrees.	required
tags	dict	Tag dict in osmnx convention, e.g. {"building": True}. OHM uses the OSM tag vocabulary, so the same dicts work for both providers.	required
date	str | None	Optional ISO 8601 date string (e.g. "2020-01-01" or "2020-01-01T00:00:00Z"). With provider="osm" this queries the OSM database as it existed at that point in time via the Overpass [date:"..."] filter. With provider="ohm" it selects features that existed in the real world at that date, by filtering on their start_date/end_date tags (features without a start_date are treated as always existing); partial dates ("1900") and BCE years ("-0500") are supported.	None
provider	str	"osm" (OpenStreetMap, default) or "ohm" (OpenHistoricalMap).	'osm'

Returns:

Type	Description
GeoDataFrame	Features clipped to the exact bbox polygon. Returns an empty GeoDataFrame (with a geometry column) if no features are found.

Source code in osm_rasterizer/fetch.py

def fetch_features(
    bbox: tuple[float, float, float, float],
    tags: dict,
    date: str | None = None,
    provider: str = "osm",
) -> gpd.GeoDataFrame:
    """Fetch OSM or OpenHistoricalMap features for a WGS84 bounding box.

    Parameters
    ----------
    bbox:
        ``(minx, miny, maxx, maxy)`` in WGS84 degrees.
    tags:
        Tag dict in osmnx convention, e.g. ``{"building": True}``.  OHM uses
        the OSM tag vocabulary, so the same dicts work for both providers.
    date:
        Optional ISO 8601 date string (e.g. ``"2020-01-01"`` or
        ``"2020-01-01T00:00:00Z"``).  With ``provider="osm"`` this queries
        the OSM database as it existed at that point in time via the
        Overpass ``[date:"..."]`` filter.  With ``provider="ohm"`` it selects
        features that existed in the real world at that date, by filtering
        on their ``start_date``/``end_date`` tags (features without a
        ``start_date`` are treated as always existing); partial dates
        (``"1900"``) and BCE years (``"-0500"``) are supported.
    provider:
        ``"osm"`` (OpenStreetMap, default) or ``"ohm"`` (OpenHistoricalMap).

    Returns
    -------
    geopandas.GeoDataFrame
        Features clipped to the exact bbox polygon.  Returns an empty
        GeoDataFrame (with a geometry column) if no features are found.
    """
    if provider not in VALID_PROVIDERS:
        raise ValueError(f"provider must be one of {VALID_PROVIDERS}, got {provider!r}")

    minx, miny, maxx, maxy = bbox

    original_settings = ox.settings.overpass_settings
    original_url = ox.settings.overpass_url
    if provider == "ohm":
        # OHM's [date:] selects database edit history, not historical
        # reality, so the date is applied post-fetch via start_date/end_date.
        ox.settings.overpass_url = OHM_OVERPASS_URL
    elif date:
        dt = date if "T" in date else f"{date}T00:00:00Z"
        ox.settings.overpass_settings = f'[out:json][timeout:180][date:"{dt}"]'

    # osmnx 2.x uses (west, south, east, north) = (minx, miny, maxx, maxy),
    # which matches our convention directly. osmnx 1.x used (north, south, east, west).
    try:
        gdf = ox.features_from_bbox(bbox=(minx, miny, maxx, maxy), tags=tags)
    except InsufficientResponseError:
        return gpd.GeoDataFrame(geometry=[], crs="EPSG:4326")
    finally:
        ox.settings.overpass_settings = original_settings
        ox.settings.overpass_url = original_url

    if provider == "ohm" and date:
        gdf = filter_by_date(gdf, date)

    if gdf.empty:
        return gpd.GeoDataFrame(geometry=[], crs="EPSG:4326")

    # Clip to exact bbox to remove any features partially outside
    clip_poly = box(minx, miny, maxx, maxy)
    return gdf.clip(clip_poly)

get_utm_crs

Auto-detect best UTM CRS for a WGS84 bbox.

Parameters:

Name	Type	Description	Default
bbox	tuple[float, float, float, float]	(minx, miny, maxx, maxy) in WGS84 (EPSG:4326).	required

Returns:

Type	Description
CRS	The best-fit UTM CRS for the given bounding box.

Source code in osm_rasterizer/crs.py

def get_utm_crs(bbox: tuple[float, float, float, float]) -> CRS:
    """Auto-detect best UTM CRS for a WGS84 bbox.

    Parameters
    ----------
    bbox:
        ``(minx, miny, maxx, maxy)`` in WGS84 (EPSG:4326).

    Returns
    -------
    pyproj.CRS
        The best-fit UTM CRS for the given bounding box.
    """
    minx, miny, maxx, maxy = bbox
    aoi = AreaOfInterest(
        west_lon_degree=minx,
        south_lat_degree=miny,
        east_lon_degree=maxx,
        north_lat_degree=maxy,
    )
    results = query_utm_crs_info(datum_name="WGS 84", area_of_interest=aoi)
    if not results:
        raise ValueError(f"No UTM CRS found for bbox {bbox}")
    best = results[0]
    return CRS.from_authority(best.auth_name, best.code)