diff --git a/orbit-georef/README.md b/orbit-georef/README.md index 11d2dda..7e92885 100644 --- a/orbit-georef/README.md +++ b/orbit-georef/README.md @@ -1,6 +1,6 @@ # orbit-georef -Standalone Python library for pixel↔geo coordinate transformation using georeferencing parameters exported from ORBIT. +Standalone Python library for pixel↔geo coordinate transformation using georeferencing parameters exported from ORBIT. Supports transformation to any coordinate system that can be expressed via [PROJ strings](https://proj.org/en/stable/usage/quickstart.html). ## Installation @@ -15,7 +15,7 @@ pip install -e . ## Usage -### Load from ORBIT export +### Load from ORBIT export and convert coordinates to desired coordinate system ```python from orbit_georef import load_georef @@ -26,10 +26,13 @@ georef = load_georef("georef_params.json") # Convert pixel to geographic coordinates lon, lat = georef.pixel_to_geo(1234.5, 567.8) -# Convert geographic to pixel coordinates +# Convert pixel coordinates to a projected coordinate system via a PROJ string (https://proj.org/en/stable/usage/quickstart.html), e.g. WGS84 +x, y = georef.pixel_to_geo(1234.5, 567.8, proj_string="+proj=longlat +datum=WGS84 +no_defs") + +# Convert geographic to pixel coordinates (assumes WGS84 as input) px, py = georef.geo_to_pixel(12.945, 57.720) -# Batch conversion +# Batch conversion (as `pixel_to_geo()`, this supports also the optional `proj_string` argument) geo_coords = georef.pixels_to_geo_batch([(100, 200), (300, 400), (500, 600)]) # Get scale factors (meters per pixel) diff --git a/orbit-georef/pyproject.toml b/orbit-georef/pyproject.toml index f754fee..cc5f276 100644 --- a/orbit-georef/pyproject.toml +++ b/orbit-georef/pyproject.toml @@ -26,6 +26,7 @@ classifiers = [ ] dependencies = [ "numpy>=1.20.0", + "pyproj>=3.0.0", ] [project.optional-dependencies] diff --git a/orbit-georef/src/orbit_georef/transformer.py b/orbit-georef/src/orbit_georef/transformer.py index 479b5f2..63fa193 100644 --- a/orbit-georef/src/orbit_georef/transformer.py +++ b/orbit-georef/src/orbit_georef/transformer.py @@ -9,6 +9,7 @@ from typing import Any, Dict, List, Optional, Tuple import numpy as np +import pyproj from .models import ControlPoint @@ -131,31 +132,67 @@ def from_control_points( control_points=control_points, ) - def pixel_to_geo(self, pixel_x: float, pixel_y: float) -> Tuple[float, float]: + def pixel_to_geo(self, pixel_x: float, pixel_y: float, proj_string: str = None) -> Tuple[float, float]: """ Convert pixel coordinates to geographic coordinates. Args: pixel_x: X coordinate in pixels pixel_y: Y coordinate in pixels + proj_string: Projection string for the output. If `None`, assumes WGS84 (lat/lon) as output. Returns: - Tuple of (longitude, latitude) in decimal degrees + Tuple of (x, y) in the target coordinate system. If `proj_string` is `None`, returns (longitude, latitude). """ if self.method == "homography": # Homography: pixel → meters → geo pixel_homo = np.array([pixel_x, pixel_y, 1.0]) ground_homo = self.transform_matrix @ pixel_homo - east = ground_homo[0] / ground_homo[2] - north = ground_homo[1] / ground_homo[2] - lat, lon = self._meters_to_latlon(east, north) - return lon, lat + mx = ground_homo[0] / ground_homo[2] + my = ground_homo[1] / ground_homo[2] + + if proj_string: + return self._meters_to_proj(mx, my, proj_string) + + else: + lat, lon = self._meters_to_latlon(mx, my) + return lon, lat else: # Affine: direct transformation + if proj_string is not None: + raise NotImplementedError("Affine transformation not supported in combination with `proj` string.") point = np.array([pixel_x, pixel_y, 1.0]) result = self.transform_matrix @ point return result[0], result[1] # lon, lat + def _meters_to_proj(self, mx: float, my: float, proj_string: str) -> Tuple[float, float]: + """ + Convert meter offsets to a projected coordinate system. + + Args: + mx: Meter offset x (easting) + my: Meter offset y (northing) + proj_string: PROJ string for the target coordinate system + + Returns: + Tuple of (x, y) in the target coordinate system + """ + # The local meter-based system used by ORBIT is an Equirectangular projection + # centered at the reference lat/lon, with the latitude of true scale + # set to the reference latitude. This is the correct projection to use + # as it precisely matches how the (mx, my) coordinates were calculated. + # See: https://proj.org/operations/projections/eqc.html + source_proj_string = ( + f"+proj=eqc +lat_ts={self.reference_lat} +lat_0={self.reference_lat} " + f"+lon_0={self.reference_lon} +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" + ) + + # Create a transformer from the local CRS to the target CRS + transformer = pyproj.Transformer.from_crs( + source_proj_string, proj_string, always_xy=True + ) + return transformer.transform(mx, my) + def geo_to_pixel(self, longitude: float, latitude: float) -> Tuple[float, float]: """ Convert geographic coordinates to pixel coordinates. @@ -180,18 +217,21 @@ def geo_to_pixel(self, longitude: float, latitude: float) -> Tuple[float, float] return result[0], result[1] def pixels_to_geo_batch( - self, points: List[Tuple[float, float]] + self, + points: List[Tuple[float, float]], + proj_string: str = None, ) -> List[Tuple[float, float]]: """ Convert multiple pixel coordinates to geographic coordinates. Args: points: List of (pixel_x, pixel_y) tuples + proj_string: Optional PROJ string for output coordinate system. If none, assumes WGS84 for output. Returns: - List of (longitude, latitude) tuples + List of (longitude, latitude) or (x, y) tuples in the target coordinate system. """ - return [self.pixel_to_geo(x, y) for x, y in points] + return [self.pixel_to_geo(x, y, proj_string=proj_string) for x, y in points] def geo_to_pixels_batch( self, points: List[Tuple[float, float]]