From b5a8bd69d57948b3169f942050cff9b93b7b5bb0 Mon Sep 17 00:00:00 2001 From: Bastian Havers-Zulka Date: Tue, 17 Feb 2026 14:40:13 +0100 Subject: [PATCH 1/7] add pyproj --- orbit-georef/pyproject.toml | 1 + orbit-georef/src/orbit_georef/test.py | 32 ++++++++++++++ orbit-georef/src/orbit_georef/transformer.py | 45 +++++++++++++++++--- 3 files changed, 72 insertions(+), 6 deletions(-) create mode 100644 orbit-georef/src/orbit_georef/test.py 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/test.py b/orbit-georef/src/orbit_georef/test.py new file mode 100644 index 0000000..60f6a79 --- /dev/null +++ b/orbit-georef/src/orbit_georef/test.py @@ -0,0 +1,32 @@ +from orbit_georef import GeoTransformer, ControlPoint +import random +import numpy as np + +# Define control points +control_points = [ + ControlPoint(pixel_x=100, pixel_y=100, longitude=12.94, latitude=57.72), + ControlPoint(pixel_x=200, pixel_y=100, longitude=12.95, latitude=57.72), + ControlPoint(pixel_x=100, pixel_y=200, longitude=12.94, latitude=57.71), + ControlPoint(pixel_x=200, pixel_y=200, longitude=12.95, latitude=57.71), +] + +# Create transformer (uses homography by default) +georef = GeoTransformer.from_control_points(control_points, method="homography") + +# Test pixel to geo +lon, lat = georef.pixel_to_geo(100, 100) + +print(f"Pixel (100, 100) -> Geo ({lon}, {lat})") + +lon, lat = georef.pixel_to_geo(100, 100, proj_string="+proj=longlat +datum=WGS84 +no_defs") + +print(f"Pixel (100, 100) -> Geo ({lon}, {lat}) with explicit proj string") + + +for i in range(10): + x, y = random.randint(0, 300), random.randint(0, 300) + without_proj = georef.pixel_to_geo(x, y) + with_proj = georef.pixel_to_geo(x, y, proj_string="+proj=longlat +datum=WGS84 +no_defs") + diffs = np.array(without_proj) - np.array(with_proj) + + print(f"Pixel ({x}, {y}) -> diffs: {diffs}") \ No newline at end of file diff --git a/orbit-georef/src/orbit_georef/transformer.py b/orbit-georef/src/orbit_georef/transformer.py index 479b5f2..d05eab4 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,63 @@ 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 + my: Meter offset y + proj_string: PROJ string for the target coordinate system + Returns: + Tuple of (x, y) in the target coordinate system + """ + # The local meter-based system is an equirectangular projection + # centered at the reference lat/lon. We can define a PROJ string for it. + # See: https://proj.org/operations/projections/eqc.html + source_proj_string = ( + f"+proj=eqc +lat_ts=0 +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. From ee4fb0cbe4100c8447a8434fa308af9df9dbd833 Mon Sep 17 00:00:00 2001 From: Bastian Havers-Zulka Date: Tue, 17 Feb 2026 15:18:35 +0100 Subject: [PATCH 2/7] set latitude of true scale to gain much higher precision --- orbit-georef/src/orbit_georef/transformer.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/orbit-georef/src/orbit_georef/transformer.py b/orbit-georef/src/orbit_georef/transformer.py index d05eab4..8fd813c 100644 --- a/orbit-georef/src/orbit_georef/transformer.py +++ b/orbit-georef/src/orbit_georef/transformer.py @@ -168,18 +168,22 @@ def pixel_to_geo(self, pixel_x: float, pixel_y: float, proj_string: str = None) 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 - my: Meter offset y + 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 is an equirectangular projection - # centered at the reference lat/lon. We can define a PROJ string for it. + # 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=0 +lat_0={self.reference_lat} " + 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" ) From df1bff5c07b92e66582d0bc5a7e31907cff805db Mon Sep 17 00:00:00 2001 From: Bastian Havers-Zulka Date: Tue, 17 Feb 2026 16:56:45 +0100 Subject: [PATCH 3/7] remove test file --- orbit-georef/src/orbit_georef/test.py | 32 --------------------------- 1 file changed, 32 deletions(-) delete mode 100644 orbit-georef/src/orbit_georef/test.py diff --git a/orbit-georef/src/orbit_georef/test.py b/orbit-georef/src/orbit_georef/test.py deleted file mode 100644 index 60f6a79..0000000 --- a/orbit-georef/src/orbit_georef/test.py +++ /dev/null @@ -1,32 +0,0 @@ -from orbit_georef import GeoTransformer, ControlPoint -import random -import numpy as np - -# Define control points -control_points = [ - ControlPoint(pixel_x=100, pixel_y=100, longitude=12.94, latitude=57.72), - ControlPoint(pixel_x=200, pixel_y=100, longitude=12.95, latitude=57.72), - ControlPoint(pixel_x=100, pixel_y=200, longitude=12.94, latitude=57.71), - ControlPoint(pixel_x=200, pixel_y=200, longitude=12.95, latitude=57.71), -] - -# Create transformer (uses homography by default) -georef = GeoTransformer.from_control_points(control_points, method="homography") - -# Test pixel to geo -lon, lat = georef.pixel_to_geo(100, 100) - -print(f"Pixel (100, 100) -> Geo ({lon}, {lat})") - -lon, lat = georef.pixel_to_geo(100, 100, proj_string="+proj=longlat +datum=WGS84 +no_defs") - -print(f"Pixel (100, 100) -> Geo ({lon}, {lat}) with explicit proj string") - - -for i in range(10): - x, y = random.randint(0, 300), random.randint(0, 300) - without_proj = georef.pixel_to_geo(x, y) - with_proj = georef.pixel_to_geo(x, y, proj_string="+proj=longlat +datum=WGS84 +no_defs") - diffs = np.array(without_proj) - np.array(with_proj) - - print(f"Pixel ({x}, {y}) -> diffs: {diffs}") \ No newline at end of file From 313d1bfcb6de5984cfb0b21ca229ed31059fa76c Mon Sep 17 00:00:00 2001 From: Bastian Havers-Zulka Date: Tue, 17 Feb 2026 17:01:49 +0100 Subject: [PATCH 4/7] adapt README to changes --- orbit-georef/README.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/orbit-georef/README.md b/orbit-georef/README.md index 11d2dda..71d04c6 100644 --- a/orbit-georef/README.md +++ b/orbit-georef/README.md @@ -26,7 +26,10 @@ 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 (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 From a7f327e71fd1abd8ff9c03bf09c40f4c04351e5c Mon Sep 17 00:00:00 2001 From: Bastian Havers-Zulka Date: Tue, 17 Feb 2026 20:02:11 +0100 Subject: [PATCH 5/7] improve README to mention PROJ strings --- orbit-georef/README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/orbit-georef/README.md b/orbit-georef/README.md index 71d04c6..3d76c63 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 @@ -26,7 +26,7 @@ georef = load_georef("georef_params.json") # Convert pixel to geographic coordinates lon, lat = georef.pixel_to_geo(1234.5, 567.8) -# Convert pixel coordinates to a projected coordinate system (e.g., WGS84) +# 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) From 7ad0ab03df5d42b67b6a475edc5864a5ab029529 Mon Sep 17 00:00:00 2001 From: Bastian Havers-Zulka Date: Tue, 17 Feb 2026 20:09:56 +0100 Subject: [PATCH 6/7] add support for batch conversion via proj strings --- orbit-georef/src/orbit_georef/transformer.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/orbit-georef/src/orbit_georef/transformer.py b/orbit-georef/src/orbit_georef/transformer.py index 8fd813c..63fa193 100644 --- a/orbit-georef/src/orbit_georef/transformer.py +++ b/orbit-georef/src/orbit_georef/transformer.py @@ -217,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]] From 10dc2bc9c820e8730bff6b6e200b3ca619d8eefb Mon Sep 17 00:00:00 2001 From: Bastian Havers-Zulka Date: Tue, 17 Feb 2026 20:10:13 +0100 Subject: [PATCH 7/7] mention proj_strings for batch conversion --- orbit-georef/README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/orbit-georef/README.md b/orbit-georef/README.md index 3d76c63..7e92885 100644 --- a/orbit-georef/README.md +++ b/orbit-georef/README.md @@ -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 @@ -32,7 +32,7 @@ x, y = georef.pixel_to_geo(1234.5, 567.8, proj_string="+proj=longlat +datum=WGS8 # 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)