Quickstart#

A short, step-by-step path: choose a mapper, build a camera, create a geo-referenced image, project 3D points, map pixels.

What the code does#

  • Enables pyproj network so compound CRSs (with geoid/vertical grids) can resolve required PROJ data.

  • Sets a horizontal-plane mapper at altitude 250 m in UTM 33 with vertical reference (EGM2008) for mapper and image.

  • Defines a calibrated OpenCV pinhole camera.

  • Specifies image pose (position + orientation) and its CRS.

  • Instantiates ImagePerspective combining camera, CRS, pose, and mapper.

  • Calls project to project 3D coordinates into the image.

  • Calls map_points to cast rays through pixels and intersect them with the mapper, returning 3D coords, normals, and GSD estimates.

Hints#

  • You can specify different CRSs for image and mapper. weitsicht will estimate transformations via pyproj.

  • Add distortion parameters (k1-k4, p1, p2) to CameraOpenCVPerspective if known; otherwise they default to zero.

  • For terrain-aware mapping, use a DEM/mesh mapper instead of a plane.

  • Replace hardcoded pose with values from metadata (see Metadata for Camera and Pose Extraction) for drone imagery workflows.

1) Imports#

import numpy as np
import pyproj.network
from pyproj import CRS

from weitsicht import CameraOpenCVPerspective, ImagePerspective, MappingHorizontalPlane, Rotation

# PROJ activate network grid/data
pyproj.network.set_network_enabled(True)  # pyright: ignore[reportPrivateImportUsage]

2) Coordinate system#

# Coordinate System we want to use
crs = CRS("EPSG:25833+3855")
# This will be UTM/Zone 33 North with vertical reference EGM2008

3) Ground model (mapper)#

Use a horizontal plane at a constant altitude as the intersection surface.

mapper = MappingHorizontalPlane(plane_altitude=250.0, crs=crs)

4) Camera model#

Provide intrinsics and distortion. Width/height belong to the image.

cam = CameraOpenCVPerspective(
    width=11648,
    height=8736,
    fx=12906.9238,
    fy=12906.9238,
    cx=5858.437,
    cy=4400.116,
    k1=-0.00565896,
    k2=0.0334872,
    k3=-0.0683847,
    p1=0.000933359,
    p2=0.00065228,
)

5) Geo-referenced image#

Attach pose, CRS, camera, and mapper. Once set, image.is_geo_referenced is True.

image = ImagePerspective(
    width=11648,
    height=8736,
    camera=cam,
    crs=crs,
    position=np.array([601968.0, 5340368.0, 320.0]),  # in image CRS
    orientation=Rotation.from_opk_degree(omega=10, phi=10, kappa=90.0),
    mapper=mapper,
)

6) Project world points to pixels#

project transforms coordinates to the image CRS, applies pose + intrinsics, and returns pixels + a validity mask.

pts_world = np.array([[601932.34, 5340348.43, 254.1], [601948.43, 5340359.13, 251.3], [601969.76, 5340392.42, 250.1]])

proj = image.project(pts_world, crs_s=image.crs)

print("\nProjection of 3D coordinates into image")
if proj.ok is True:
    pixels_in_frame = proj.pixels[proj.mask]
    print("Valid pixels: ", pixels_in_frame)
    if not np.all(proj.mask):
        print("Outside image: ", proj.pixels[~proj.mask])
else:
    print("Projecting coordinates did not work")
    print(proj.error)

7) Map pixels back to 3D#

map_points casts rays through the pixels and intersects them with the mapper. The mapping result includes coordinates, normals, gsd and gsd_per_point.

print("\nMapping image points")
result_mapping = image.map_points([[1000, 800], [2500, 1600]])
if result_mapping.ok is True:
    coords_3d = result_mapping.coordinates[result_mapping.mask]
    normals = result_mapping.normals[result_mapping.mask]
    gsd = result_mapping.gsd  # mean ground sampling distance in mapper CRS units
    gsd_per_point = (
        result_mapping.gsd_per_point[result_mapping.mask] if result_mapping.gsd_per_point is not None else None
    )
    print("coordinates mapped:", coords_3d)
    print("normals:", normals)
    print("gsd:", gsd)
    print("gsd_per_point:", gsd_per_point)

8) Inspect issues#

Both proj and the mapping result carry issues and an ok flag—log them to catch partial failures.

print("\n Checks:")
if not result_mapping.ok:
    print("Projection problems:", result_mapping.issues)
if proj.issues:
    print("Mapping warnings:", proj.issues)

What to try next#

  • Swap the mapper for MappingRaster or MappingTrimesh for real terrain/mesh intersections.

  • Batch process many images with ImageBatch and a shared mapper.

  • See Results, Masks, and Errors for return-object fields and issue codes.