.. _mapper_raster: ============= Raster Mapper ============= The raster mapper :py:class:`weitsicht.MappingRaster` uses a raster DEM/DSM to: - sample heights for given ground coordinates, and - intersect rays (e.g. from a perspective image) with the raster surface to obtain 3D hit points. Internally it uses `rasterio `__ for IO and sampling. Memory and IO behavior ====================== By default, :py:class:`~weitsicht.MappingRaster` **does not load the raster into memory**. It keeps a ``rasterio.DatasetReader`` open and reads only the height values needed at runtime (e.g. via sampling). Only when you explicitly enable window/full preloading (``preload_full_raster``, ``preload_window``) or call :py:meth:`~weitsicht.MappingRaster.load_window` is raster data read into a numpy array and stored in :attr:`~weitsicht.MappingRaster.georef_array` (a :py:class:`~weitsicht.MappingGeorefArray` instance). When :attr:`~weitsicht.MappingRaster.georef_array` is present, :py:meth:`~weitsicht.MappingRaster.map_coordinates_from_rays` and :py:meth:`~weitsicht.MappingRaster.map_heights_from_coordinates` automatically delegate to that in-memory backend. This avoids disk IO and uses the (more accurate) bilinear-patch ray intersection model implemented by :py:class:`~weitsicht.MappingGeorefArray`. When to use this mapper ======================= Use :py:class:`~weitsicht.MappingRaster` when you have a DEM/DSM and want terrain-aware pixel-to-ground mapping: - drone imagery over terrain - wildlife surveys over variable relief - any use-case where a flat plane is not a good approximation Initialization requirements ============================ You must provide a readable raster file path. Common configuration parameters: - ``raster_path``: path to a file supported by rasterio (GeoTIFF is typical). - ``index_band``: required if the dataset has more than one band (bands start at 1). - ``crs``: optional override of the raster CRS. - ``force_no_crs``: force CRS to ``None`` (advanced; only safe if you know everything is in the same CRS). - ``preload_full_raster`` / ``preload_window``: optional in-memory loading for faster repeated operations. If a mapper CRS is set (either from the raster metadata or via ``crs=...``), it must be 3D (have a Z axis). CRS and projections =================== The raster itself can be stored in many different coordinate systems. Typical examples are: - projected CRSs in meters (UTM, Lambert Conformal Conic, national grids, ...), - geographic (ellipsoidal) CRSs in degrees (lon/lat) with a height axis. ``weitsicht`` can work across different CRSs by transforming *points* between the input CRS (``crs_s``) and the mapper CRS. See :doc:`../transformation` for details. .. important:: Ray/line intersection math assumes an approximately Euclidean coordinate system with linear units in XY. For reliable ray intersections, provide rays in a projected (meter-based) 3D CRS where Z represents height. Do not run ray intersections directly in lon/lat degrees; if your poses are in geographic coordinates, transform them to a suitable projected CRS first, and transform results back afterwards if needed. Core methods ============ Height sampling (coordinates → heights) ------------------------------------------ ``map_heights_from_coordinates(...)`` takes (x, y) (or (x, y, z)) coordinates and returns 3D coordinates with an interpolated height from the raster. Notes: - bilinear interpolation is used (a 2×2 window around the target pixel), - coordinates outside the raster are masked and return ``Issue.OUTSIDE_RASTER``. Ray intersections (rays → surface hits) ------------------------------------------ ``map_coordinates_from_rays(...)`` intersects rays with the raster surface and returns 3D hit points. High-level behavior: - rays are tested per input element, - rays that never intersect the raster surface in a valid direction are masked and may add ``Issue.WRONG_DIRECTION``, - rays leaving the raster extent are masked and add ``Issue.OUTSIDE_RASTER``. - rays touching raster no-data cells (holes) are masked and may add ``Issue.RASTER_NO_DATA``, - rays where the iterative intersection does not converge are masked and may add ``Issue.MAX_ITTERATION``. How ray intersections work (iterative ray–plane intersection) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The disk-backed raster mapper uses an *iterative* approach to find an intersection between a 3D ray and a raster height field. Conceptually, the raster is treated as a function: .. math:: z = f(x, y) and for each ray the mapper searches for a parameter :math:`t` such that: .. math:: \text{ray}_z(t) = f(\text{ray}_x(t), \text{ray}_y(t)) Implementation idea (in words): 1. Choose an initial plane height guess :math:`z_i` between the camera height and the raster height below the ray start. 2. Intersect the ray with the horizontal plane :math:`z = z_i` (ray–plane intersection). 3. Transform that intersection point into the raster CRS and sample the raster height at its :math:`(x, y)`. 4. Transform the sampled height back to the source CRS (important when vertical transformations are involved). 5. Update the plane height guess toward the sampled height (a damped update is used to reduce oscillation) and repeat. Stopping criteria: - the height difference between the current plane guess and the sampled raster height falls below ~2 cm, or - a maximum iteration count is reached (then the ray is treated as not intersecting reliably). The following figure shows the idea in a vertical profile (X–Z slice): each iteration intersects the ray with a horizontal plane :math:`z = z_i`, samples the raster height at the resulting :math:`(x, y)` location, and updates the next plane height guess. .. figure:: /_static/raster_ray_iteration.svg :align: center :width: 800 :alt: Iterative ray intersection on a raster height field Iterative ray–raster intersection as used by :py:class:`weitsicht.MappingRaster`. .. note:: This approach keeps the ray direction vector in the source CRS and transforms *points* into the raster CRS for sampling. This avoids transforming direction vectors through potentially non-linear CRS operations, but it requires multiple point transformations and raster samples per ray. In-memory window (optional) =========================== If you repeatedly intersect many rays within a limited raster extent, you can load a raster subset into memory. ``MappingRaster.load_window(...)`` populates ``mapper.georef_array`` (and also returns it) as a :py:class:`~weitsicht.MappingGeorefArray` instance. Once loaded, ``MappingRaster.map_*`` methods automatically use this in-memory backend. Convenience helpers: - ``mapper.backend`` returns the active backend (``MappingGeorefArray`` if loaded, otherwise ``MappingRaster``). - ``mapper.georef_mapper`` returns the cached backend or raises if no window/full raster is loaded. .. important:: If you cache only a window, mapping operations are limited to that cached extent. Points/rays outside the cached window will be reported as ``Issue.OUTSIDE_RASTER``. For details about the in-memory backend (including its bilinear patch ray intersection model), see :doc:`georef_array`. Minimal example =============== .. code-block:: python import numpy as np from pyproj import CRS from weitsicht import MappingRaster mapper = MappingRaster( raster_path="dem.tif", crs=CRS.from_epsg(32633), # optional override; otherwise read from the file index_band=1, # required if the raster has multiple bands ) # Height sampling xy = np.array([[601932.0, 5340348.0], [601948.0, 5340359.0]]) res_h = mapper.map_heights_from_coordinates(xy, crs_s=mapper.crs) # Ray intersections (origins + direction vectors must be in the same source CRS) # hits = mapper.map_coordinates_from_rays(ray_vectors, ray_starts, crs_s=image_crs) Related pages ============= - :doc:`horizontal_plane` (simpler / flat model) - :doc:`georef_array` (in-memory raster window backend) - :doc:`mesh` (complex surfaces via triangle meshes) - :doc:`../transformation` (CRS and vertical transformation notes)