georay

Render georeferenced image, DSM and land-cover references for calibrated camera poses.

georay connects camera geometry with geospatial raster data. Given a calibrated camera, a 6D pose, or an explicit UTM footprint, it computes the needed reference area, selects suitable source tiles, builds mosaics, and returns image/DSM cutouts (and optionally pixel-aligned land-cover masks) together with coordinate mappings and tile metadata.

Features

• Pluggable raster sources: German state data via georaffer, or any local GeoTIFF/COG/VRT files in any projected CRS
• Camera-footprint computation from calibrated pinhole cameras and Pose6D
• DSM-backed ray marching for terrain-aware footprint estimation
• Rendering from explicit BoundingBox2D footprints
• Land-cover masks (e.g. CORINE Land Cover) from vector or class raster data, pixel-aligned with the orthophoto, with CLC class tables and colours
• Orthophoto and DSM tile selection by year, season, and region preference
• Tile retrieval through georaffer, with quiet external output by default
• Mosaic construction, axis-aligned slicing, perspective warping, and metadata
• Rich-based level logging and progress spinners for longer operations

Installation

Install from GitHub:

python3 -m pip install "georay @ git+https://github.com/hlzmnhnry/georay.git"

or in editable mode from a clone:

python3 -m pip install -e .

For development and tests:

python3 -m pip install -e ".[dev]"
python3 -m pytest

For official German state data install the germany extra, which pulls in georaffer:

python3 -m pip install -e ".[germany]"

Land-cover masks from vector data (Shapefile, GeoPackage) need the vector extra (pyogrio, shapely, pyproj):

python3 -m pip install -e ".[vector]"

Raster Sources

georay renders from a regular grid of square tiles in a projected CRS; tile (x, y) covers [x * size, (x + 1) * size) in easting and northing. A raster source provides the tiles:

• GeorafferSource(regions, imagery_from, dsm_from) downloads and converts orthophotos and surface models of NRW, RLP, BW, BY, TH (ETRS89 / UTM 32N) or BB (UTM 33N); TH needs georaffer > 0.2.0 through georaffer. Passing regions=... to the renderers creates it implicitly.
• RasterFileSource(image_files, dsm_files, crs=None, tile_size_m=1000) uses your own files anywhere in the world. They are reprojected into the tile grid (default CRS: that of the first image file) and cached as GeoTIFFs. DSM heights are used as they are, so poses must use the same vertical reference.

from georay.raster.sources import RasterFileSource

source = RasterFileSource(["naip_ortho.tif"], ["dsm.tif"], crs="EPSG:32617")
renderer = ReferenceViewRenderer(camera=camera, source=source)

Any object with crs, tile_size_m, query(coords, tile_type) and get_many(tiles, resolution) works as a source (georay.raster.sources.RasterSource). When mosaicking, georay checks every tile's georeferencing against its grid cell.

Quickstart

The high-level entry point is ReferenceViewRenderer. It can render both an orthophoto reference and a DSM reference from a pose:

from datetime import date

from georay import configure_logging
from georay.geometry.pinhole_camera import PinholeCamera
from georay.raster.reference_tile_selector import ClosestReferenceTileSelector
from georay.reference_view_renderer import ReferenceViewRenderer
from georay.spatial.bounding_box_transform import NorthAlignedBoundingBoxTransform
from georay.types.pose import Pose6D
from georay.types.tile import RankingOrder

configure_logging("INFO")

camera = PinholeCamera(
    image_width=1600,
    image_height=1100,
    fx=1320.0,
    fy=1319.3,
    cx=801.0,
    cy=560.6,
)

selector = ClosestReferenceTileSelector(
    ranking_order=RankingOrder.YEAR_STATE_SEASON,
)

renderer = ReferenceViewRenderer(
    camera=camera,
    regions=["nrw", "rlp"],
    imagery_from=(2015, None),
    dsm_from=(2015, None),
    tile_selector=selector,
)

pose = Pose6D(rotation_matrix=rotation_world_from_camera, translation=utm_xyz)

result = renderer.render(
    pose=pose,
    pixels_per_meter=2,
    output_size=(640, 480),
    bbox_transform=NorthAlignedBoundingBoxTransform(enlargement=0.1),
    preferred_date=date(2022, 2, 23),
    preferred_regions="nrw",
    render_image=True,
    render_dsm=True,
)

image = result.image.rasters[0] if result.image else None
dsm = result.dsm.rasters[0] if result.dsm else None

For a complete runnable example, see examples/simple_reference_from_bbox.py.

Rendering From A Footprint

If the footprint is already known, bypass pose-based footprint computation and render directly from a BoundingBox2D:

reference = renderer.render_from_bbox(
    bbox=bbox,
    pose=pose,
    pixels_per_meter=2,
    bbox_transform=NorthAlignedBoundingBoxTransform(enlargement=0.1),
    preferred_regions="nrw",
    render_image=True,
    render_dsm=True,
    dsm_interpolation="nearest",
)

render_from_bbox() is useful when footprints were precomputed, manually curated, or loaded from a sequence file.

Land-Cover Masks

A label_source provides class codes on the same tile grid as the orthophoto and DSM. With render_labels=True the renderer returns them as a third layer, cut out and warped exactly like the orthophoto (nearest neighbour, so classes are never mixed); pixels without a class are 0.

from georay.landcover import colorize, legend
from georay.raster.vector_source import VectorFileSource

labels = VectorFileSource(["landcover.gpkg"], "class_code", crs="EPSG:25832")
renderer = ReferenceViewRenderer(camera=camera, regions=["nrw"], imagery_from=(2015, None),
    label_source=labels)

reference = renderer.render_from_bbox(bbox, pixels_per_meter=2, render_labels=True)
codes = reference.labels.rasters[0]      # uint16 class codes, aligned with reference.image
preview = colorize(codes)                # RGB
print(legend(codes))                     # classes with names, colours and shares

• VectorFileSource(files, value_field, crs=...) rasterizes polygons of any OGR format (Shapefile, GeoPackage, ...) in any CRS: a pixel gets the class of the polygon containing its centre (all_touched=True for any overlap). Text classes are mapped with value_map. Several files of one dataset can be passed at once.
• RasterFileSource(label_files=[...]) uses class rasters (GeoTIFF/COG/VRT) instead, e.g. CORINE Land Cover 100 m for Europe or ESA WorldCover.
• Rasterized tiles are cached. Reading is much faster from a GeoPackage, which has a spatial index: georay.raster.vector_source.merge_to_geopackage(shapefiles, "landcover.gpkg").

CORINE Land Cover (CLC)

georay.landcover knows the CLC nomenclature: CLC_CLASSES (code, name, level 1 and 2), colorize(codes, palette) with the colours of the LTC dataset ("ltc", default; readable with pygeon) or of the official CLC legend ("eea"), clc_level(codes, 1 | 2) for the CLC hierarchy, aggregate(codes) for the coarse classes of the LTC dataset, class_shares and legend.

Germany, CLC5-2018 (5 ha, vector, BKG): free download from the BKG product page, e.g. directly in ETRS89 / UTM 32N (shapefiles, 1.4 GB): https://daten.gdz.bkg.bund.de/produkte/dlm/clc5_2018/aktuell/clc5_2018.utm32s.shape.zip. The class code is the field CLC18; clc5_source sets this up:

from glob import glob
from georay.landcover import clc5_source

labels = clc5_source(glob("clc5/clc5_class*.shp"), crs="EPSG:25832")

The data are licensed under the Datenlizenz Deutschland – Namensnennung – Version 2.0. Anything you publish or pass on that contains or is derived from them (including rendered masks) must carry the source note "© GeoBasis-DE / BKG (year of your download), dl-de/by-2-0" and a note that the data were modified (e.g. "rasterized").

Europe, CLC (100 m raster or vector, Copernicus): available after free registration from the Copernicus Land Monitoring Service; use the GeoTIFF with RasterFileSource(label_files=[...]). Source note: "© European Union, Copernicus Land Monitoring Service , European Environment Agency (EEA)".

Loading Explicit Tiles

Load one or more tile coordinates directly without building a mosaic:

from georay.reference_tile_renderer import ReferenceTileRenderer

tile_renderer = ReferenceTileRenderer(
    regions=["nrw", "rlp"],
    imagery_from=(2015, None),
    dsm_from=(2015, None),
    tile_selector=selector,
)

tiles = tile_renderer.render_tiles(
    tile_coords=[(357, 5612), (400, 5700)],
    pixels_per_meter=1,
    preferred_regions="nrw",
)

image = tiles[(357, 5612)].image.rasters[0]
dsm = tiles[(357, 5612)].dsm.rasters[0]

Tile coordinates are the 1 km grid indices, so (357, 5612) covers eastings 357000..358000 and northings 5612000..5613000. Explicit tiles are loaded individually and may be disconnected.

Create a colored NumPy point cloud directly from selected image and DSM tiles:

cloud = tile_renderer.render_tile_point_cloud(
    tile_coords=[(357, 5612), (400, 5700)],
    pixels_per_meter=1,
)

The returned array has shape (N, 6) with columns x, y, z, r, g, b. Point coordinates represent raster pixel centers in UTM coordinates. Pixels with non-finite DSM elevations are omitted. Point-cloud generation requires the tile selector to select exactly one image and one DSM variant per tile. ReferenceViewRenderer still exposes render_tiles() and render_tile_point_cloud() as compatibility helpers, but direct tile workflows do not require a camera model.

Logging And Progress

Logging is quiet by default. Enable Rich-formatted logs and progress spinners explicitly:

from georay import configure_logging

configure_logging("INFO")     # DEBUG, INFO, WARNING, ERROR, CRITICAL
configure_logging("DEBUG")    # includes suppressed georaffer output
configure_logging("OFF")      # disables georay logs and spinners

You can also configure logging through the renderer:

renderer = ReferenceViewRenderer(
    camera=camera,
    regions=["nrw", "rlp"],
    imagery_from=(2015, None),
    tile_selector=selector,
    log_level="INFO",
    quiet_georaffer=True,
)

Or from the environment before importing georay:

export GEORAY_LOG_LEVEL=INFO

quiet_georaffer=True is the default. Direct stdout, stderr, and warnings from georaffer calls are suppressed during tile queries and retrievals; with DEBUG logging they are emitted through the normal logger for diagnostics.

Core Concepts

PinholeCamera

Represents an OpenCV-style calibrated pinhole camera. Distortion models are pluggable; no distortion is used when none is supplied.

Pose6D

Stores a 3x3 rotation matrix and 3D translation (camera to world, OpenCV camera axes). Translations are expected in the source's CRS, heights in the DSM's vertical reference.

BoundingBox2D

Stores a quadrilateral footprint in world coordinates, ordered as upper-left, upper-right, lower-right, lower-left. Its utm_zone accepts UTMZone, "utm32" or any "EPSG:<code>"; bbox.crs returns the EPSG code.

ReferenceTileSelector

Chooses source tile variants for each tile coordinate. The default ClosestReferenceTileSelector can rank by acquisition year, season, and region.

BoundingBoxTransform

Adapts footprints before rendering. Common transforms include direct use, proportional enlargement, north-aligned envelopes, and yaw-aligned rectangles.

Outputs

ReferenceViewRenderer returns a ReferenceViewResult with optional image, dsm and labels fields. Each rendered layer contains:

• rasters: one raster per selected reference variant; missing DSM values are NaN
• pixel_to_utm / utm_to_pixel: convert single points between cutout pixels and world coordinates
• pixel_to_utm_many / utm_to_pixel_many: the same for (N, 2) arrays
• pixel_to_utm_matrix: the mapping as a 3x3 homogeneous matrix
• meta: source tile metadata for each variant
• footprint: the rendered footprint
• mosaic_pixel_bbox: the footprint in mosaic pixel coordinates

Pixel coordinates are continuous with pixel (i, j) covering [i, i+1) x [j, j+1), so its value belongs to the centre (i + 0.5, j + 0.5).

By default, raster dimensions are derived from the footprint size and pixels_per_meter. Pass output_size=(width, height) to ReferenceViewRenderer.render() or render_from_bbox() to resize rendered layers to a fixed pixel resolution while keeping the coordinate mappers consistent with the resized output.

Examples

• examples/local_rasters_reference.py: render references from your own orthophoto/DSM files (no German data needed).
• examples/simple_reference_from_bbox.py: render image and DSM references from a known bounding box (German state data via georaffer).
• examples/footprints_csv_reference.py: iterate over a CSV of camera footprints (four UTM corners per image) and render orthophoto and DSM for each, with --labels also land-cover masks (e.g. --labels clc5/clc5_class*.shp); runs out of the box on 10 bundled LTC footprints (examples/data/footprints_sample.csv).

Development

Run tests:

python3 -m pytest

Optional quality checks when the tools are installed:

python3 -m ruff check .
python3 -m black --check .

The package supports Python 3.10 and newer.

License

MIT, see LICENSE. Reference data downloaded through georaffer keep the licenses of their providers.