Projections and Mapping#

ImagePerspective and ImageOrtho wrap geo-referenced imagery and expose projection + mapping helpers.

Core concepts#

  • Projection = 3D world -> image pixels (project).

  • Mapping = image pixels -> 3D coordinates (map_points, map_center_point, map_footprint).

  • Monoplotting (see Monoplotting) = single-image mapping via ray/surface intersection (pixel -> ray -> ground model/mesh hit). In weitsicht this is exactly what the map_* methods do when a mapper (plane/raster/mesh/array) is provided.

  • CRS matters: both directions transform inputs into the image CRS before optics/geo math.

  • Mapper matters: mapping calls need a ground model (plane, raster, mesh, etc.).

Projection functions#

  • ImagePerspective.project: 3D coordinates -> pixels. Does CRS transform, applies pose, projects with camera model, re-applies distortion (only where inside distortion border). Needs camera, position, orientation, crs set (is_geo_referenced).

  • ImageOrtho.project: 3D coordinates -> pixels via affine inverse transform; validity check is raster bounds.

  • ImagePerspective.pixel_to_ray_vector: pixels -> unit rays in camera CRS (optionally undistorted). Useful for custom ray tracing or mapper reuse.

  • Validity: both return a mask marking in-frame pixels; outside-frame points stay in the array but are flagged.

Mapping functions#

  • map_points (both classes): pixels -> 3D using the attached/provided mapper. Perspective computes GSD from mapped range + camera model; ortho uses its resolution.

  • map_footprint (both): builds image footprint polygon in 3D; perspective can densify edges, ortho follows raster corners.

  • map_center_point (both): maps principal point (perspective) or raster center (ortho) to 3D.

  • Mapper requirements: a MappingBase instance must be present on the image or passed into the call. Typical choices: horizontal plane, raster DEM, mesh.

  • CRS handling: pixels/rays are intersected in mapper CRS; inputs are transformed automatically via CoordinateTransformer.

How projection works (perspective)#

  1. Transform coordinates from crs_s to image CRS.

  2. Move into camera space with pose (orientation + position).

  3. Project through camera intrinsics; undistort if requested.

  4. Mark pixels outside distortion border as invalid.

How mapping works (pixels -> 3D)#

  1. (Perspective only) Pixels -> rays via pixel_to_ray_vector; origins = image position.

  2. (Ortho) Affine transform pixels to ground 2D, then mapper supplies heights if available.

  3. Mapper intersects rays/points with its surface (plane/raster/mesh/array).

  4. GSD derived from range (perspective) or stored resolution (ortho).

Returned surface normals (.normals)#

All mapping helpers return a MappingResultSuccess on success. Besides coordinates and mask, the result also provides surface normals via normals:

  • shape: (N, 3) (aligned with the input order),

  • CRS: same CRS as coordinates,

  • magnitude: unit vectors for valid entries (length ~= 1),

  • invalid entries: where mask is False, normals are nan.

Normals are mainly used to improve perspective-image GSD estimation (incidence-angle correction), but they can also be useful to derive viewing angles, local surface orientation, or to post-process footprints.

The exact meaning/quality depends on the mapper backend (plane/raster/mesh/etc.). See Surface Normals for details.

mp = img.map_points([[1000, 900], [1200, 1100]])
if mp.ok:
    coords = mp.coordinates[mp.mask]
    normals = mp.normals[mp.mask]

GSD estimation for perspective images#

For perspective images, the ground sampling distance (GSD) depends on the distance to the mapped 3D point and varies across the image. See GSD estimation (mapping results) for details on how gsd and gsd_per_point are estimated.

weitsicht reports GSD in mapping results via:

  • gsd: mean GSD over valid points (scalar)

  • gsd_per_point: per input pixel (length N)

Computation outline (per mapped point):

  1. Compute the range \(R = \|P - C\|\) from the camera position \(C\) to the mapped 3D point \(P\).

  2. Use the camera focal length in pixels (focal_length_for_gsd_in_pixel) to estimate the small-angle scale, roughly \(\mathrm{GSD} \approx R / f_{px}\).

  3. If a surface normal is available (see Surface Normals), apply an incidence-angle correction so oblique views yield larger GSD on the surface.

  4. Refine the estimate using neighbouring pixel rays (one-pixel step in x/y): if normals are valid, neighbour rays are intersected with the local tangent plane at \(P\); otherwise a range-sphere chord-length estimate is used.

Note

GSD is expressed in the linear unit of the image/mapping CRS. Use a metric CRS (meters) for physically meaningful values; in geographic CRSs (degrees) the returned values are degrees/pixel.

Orthophoto GSD and area#

For orthophotos, weitsicht assumes a constant pixel size and uses resolution:

  • gsd is set to resolution

  • gsd_per_point is filled with resolution (per input pixel)

  • area (for map_footprint) is computed from the mapped footprint polygon in the orthophoto CRS

Warning

These values are expressed in the orthophoto CRS units. If the orthophoto is in a non-metric CRS (e.g. lon/lat degrees), resolution is degrees/pixel and the polygon area is degrees^2, which is not a meaningful physical area. Reproject to a projected metric CRS before using GSD/area for measurements.

Example: project then map#

import numpy as np
from pyproj import CRS
from weitsicht.image.perspective import ImagePerspective
from weitsicht.mapping.horizontal_plane import MappingHorizontalPlane

# assume camera is already built
mapper = MappingHorizontalPlane(plane_altitude=5.0, crs=CRS.from_epsg(4979))
img = ImagePerspective(width=4000, height=3000,
                       camera=camera,
                       crs=CRS("EPSG:31256+5778"),
                       position=np.array([0, 0, 120]),
                       orientation=Rotation.from_opk_degree(omega=0.0, phi=0.1, kappa=90.0),
                       mapper=mapper)

pts_world = np.array([[5, 20, 0], [10, 25, 0]])
proj = img.project(pts_world, crs_s=img.crs)
pixels = proj.pixels[proj.mask]

mp = img.map_points([[1000, 900], [1200, 1100]])
coords_3d = mp.coordinates[mp.mask]
gsd = mp.gsd

Tips#

  • project returns undistorted pixels when to_distorted is False; distortion is re-applied only where valid_mask is True.

  • Use image_points_inside (perspective) or ImageOrtho.image_points_inside to pre-filter coordinates.

  • For single points, both projection and mapping accept 1D arrays; they are internally promoted to Nx2 / Nx3.

  • Keep mapper CRS aligned with image CRS to avoid extra transforms (or always pass crs_s).

  • For large point sets, batch your calls per image type to reuse transformers and mappers efficiently.