OrthomosaicReader

OrthomosaicReader#

class rapidtools.data_sources.OrthomosaicReader(dataset_path: str | Path)[source]#

Bases: object

Client for extracting high-resolution imagery from local orthomosaic rasters.

This class is designed to be used as a context manager for efficient batch processing. It opens the underlying TIFF file once, allowing you to rapidly extract thousands of image patches without the I/O bottleneck of opening and closing the file repeatedly.

On entering the context the reader also samples a downsampled thumbnail (or the smallest suitable overview) to compute robust global 1st/99th percentile statistics. These are used to radiometrically scale high-bit-depth imagery consistently across all extracted patches.

Variables:

  • dataset_path (Path) – Absolute path to the orthomosaic on disk.

  • width (int) – Raster width in pixels.

  • height (int) – Raster height in pixels.

  • global_min (float | None) – Robust global minimum pixel value, set on __enter__.

  • global_max (float | None) – Robust global maximum pixel value, set on __enter__.

  • dataset_extent (tuple | None) – WGS84 bounds (min_lon, min_lat, max_lon, max_lat), set on __enter__.

Example

>>> from rapidtools.data_sources import OrthomosaicReader
>>>
>>> with OrthomosaicReader('path/to/map.tif') as reader:
...     result = reader.get_image_patch(asset_coords)
...     if result is not None:
...         patch, pixel_coords, bbox = result

Methods

generate_tiles

Scan through the raster and yield populated image patches.

get_image_patch

Extract a bounded image patch for a geometric asset.

get_mosaic_bbox_wgs84

Calculate the bounding box of the entire raster dataset in WGS84.

get_raster_dimensions

Determine the width and height of the raster.

haversine_distance

Calculate great-circle distance between two points on Earth's surface.
get_mosaic_bbox_wgs84() → tuple[float, float, float, float][source]#

Calculate the bounding box of the entire raster dataset in WGS84.

Works both inside and outside the context manager. Outside the context, the file is briefly opened to read its metadata.

Returns:

tuple[float, float, float, float] – The geographic bounds as (min_lon, min_lat, max_lon, max_lat).

Example

>>> reader = OrthomosaicReader('path/to/map.tif')
>>> min_lon, min_lat, max_lon, max_lat = reader.get_mosaic_bbox_wgs84()
>>> min_lon < max_lon and min_lat < max_lat
True
get_image_patch(asset_geometry: list[tuple[float, float]], max_missing_data_ratio: float = 0.2, buffer: float | str = '20%', force_square: bool = True, min_footprint_coverage: float | None = None, pad_edges: bool | None = None, return_georef: Literal[False] = False) → tuple[Image, list[tuple[int, int]], tuple[float, float, float, float]] | None[source]#
get_image_patch(asset_geometry: list[tuple[float, float]], max_missing_data_ratio: float = 0.2, buffer: float | str = '20%', force_square: bool = True, min_footprint_coverage: float | None = None, pad_edges: bool | None = None, *, return_georef: Literal[True]) → tuple[Image, list[tuple[int, int]], tuple[float, float, float, float], PatchGeoref] | None

Extract a bounded image patch for a geometric asset.

The geometry type is inferred from the number of coordinates: one pair is a Point, two pairs a LineString, and four or more pairs whose first and last coordinates coincide form a Polygon. Any other open sequence is treated as a LineString.

Parameters:

  • asset_geometry – A list of (longitude, latitude) coordinates defining the asset (Point, Line, or Polygon).

  • max_missing_data_ratio – The maximum acceptable fraction (0.0 to 1.0) of missing pixels in the crop, judged from the dataset mask (nodata value, alpha band or internal mask). If exceeded, extraction is aborted.

  • buffer – The padding applied around the geometry. Can be a percentage (e.g., ‘20%’), an absolute CRS distance (e.g., 50.0), or a real-world distance (e.g., ‘50 ft’, ‘10 m’).

  • force_square – If True, guarantees the extracted image patch covers a perfectly square area in the real world.

  • min_footprint_coverage – Minimum fraction (0.0 to 1.0) of the asset footprint that must contain valid (non-nodata) pixels. When given, this check replaces the patch-level max_missing_data_ratio test so assets on the edge of the imagery are kept. Defaults to None.

  • pad_edges – Pad crop windows that extend beyond the raster with nodata so they keep their full requested size and shape. Defaults to True when min_footprint_coverage is given, False otherwise (window clipped to the raster).

  • return_georef – If True, a PatchGeoref with the patch’s native affine transform and CRS is returned as a fourth element.

Returns:

A tuple containing – - The extracted PIL Image in 8-bit RGB format.

  • The local pixel coordinates of the asset relative to the new image patch.

  • The WGS84 bounding box of the extracted image patch: (min_lon, min_lat, max_lon, max_lat).

  • With return_georef=True, the PatchGeoref of the patch.

Returns None if the extracted patch contains too much missing data or falls completely outside the raster bounds.

Raises:

  • RuntimeError – If called outside of a context manager block.

  • ValueError – If the asset geometry is empty or the buffer string cannot be parsed.

Example

>>> footprint = [
...     (-118.198, 34.098), (-118.197, 34.098),
...     (-118.197, 34.097), (-118.198, 34.097),
...     (-118.198, 34.098),
... ]
>>> with OrthomosaicReader('path/to/map.tif') as reader:
...     patch, pixel_coords, bbox = reader.get_image_patch(
...         footprint, buffer='5 m', force_square=False
...     )
>>> len(pixel_coords)
5
generate_tiles(patch_size: float = 512, unit: str = 'pixels', overlap_ratio: float = 0.2, max_missing_data_ratio: float = 0.8, pad_edge_tiles: bool = True, return_georef: Literal[False] = False) → Generator[tuple[Image, tuple[float, float, float, float]], None, None][source]#
generate_tiles(patch_size: float = 512, unit: str = 'pixels', overlap_ratio: float = 0.2, max_missing_data_ratio: float = 0.8, pad_edge_tiles: bool = True, *, return_georef: Literal[True]) → Generator[tuple[Image, tuple[float, float, float, float], PatchGeoref], None, None]

Scan through the raster and yield populated image patches.

The raster is traversed with a sliding window from the top-left corner. Each candidate window is first evaluated using the dataset mask so that expensive pixel reads are only performed for tiles with an acceptable amount of valid data.

Parameters:

  • patch_size – The dimension of the square patch in the specified unit (e.g., 100 for a 100x100 meter patch).

  • unit – The spatial unit of the patch size. Supported units include: ‘pixels’, ‘meters’, ‘feet’, ‘yards’, ‘kilometers’, ‘miles’.

  • overlap_ratio – Float between 0.0 and 1.0 representing the percentage of overlap.

  • max_missing_data_ratio – The maximum allowable ratio of empty/nodata pixels.

  • pad_edge_tiles – If True, pads tiles at the edges of the map with nodata to maintain a strictly uniform output array size.

  • return_georef – If True, each tile also carries a PatchGeoref with the tile’s native affine transform and CRS, which is what mask vectorization should use for anything that is not a north-up WGS84 raster.

Yields:

tuple[Image.Image, tuple[float, float, float, float]] – An 8-bit RGB PIL image and its WGS84 bounding box (min_lon, min_lat, max_lon, max_lat); a third PatchGeoref element when return_georef is set.

Raises:

  • RuntimeError – If called outside of a context manager block.

  • ValueError – If overlap_ratio is outside [0.0, 1.0) or the unit is unsupported.

Example

>>> with OrthomosaicReader('path/to/map.tif') as reader:
...     tiles = list(
...         reader.generate_tiles(patch_size=32, overlap_ratio=0.0)
...     )
>>> tiles[0][0].size
(32, 32)
get_raster_dimensions(unit: str = 'pixels') → tuple[float, float][source]#

Determine the width and height of the raster.

For non-pixel units the dimensions are computed as great-circle distances across the WGS84 extent: width along the central latitude and height along the central longitude.

Parameters:

unit – The desired unit of measurement (e.g., ‘pixels’, ‘m’, ‘feet’, ‘km’). Defaults to ‘pixels’.

Returns:

tuple[float, float] – The (width, height) of the raster.

Raises:

  • RuntimeError – If called outside of a context manager block.

  • ValueError – If an unsupported unit is provided.

Example

>>> with OrthomosaicReader('path/to/map.tif') as reader:
...     width_px, height_px = reader.get_raster_dimensions()
...     width_m, height_m = reader.get_raster_dimensions('m')
>>> width_px
64.0
static haversine_distance(lon1: float, lat1: float, lon2: float, lat2: float) → float[source]#

Calculate great-circle distance between two points on Earth’s surface.

Parameters:

  • lon1 – Longitude of the first point in degrees.

  • lat1 – Latitude of the first point in degrees.

  • lon2 – Longitude of the second point in degrees.

  • lat2 – Latitude of the second point in degrees.

Returns:

float – The geographic distance in meters.

Example

>>> d = OrthomosaicReader.haversine_distance(0.0, 0.0, 1.0, 0.0)
>>> round(d / 1000.0, 1)
111.2