From 0c97110afb1fbca7c710857053d99fc636921e77 Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Tue, 9 Jun 2026 22:03:05 -0500 Subject: [PATCH 1/5] initial framework for helmert transformation tests --- ohmg/georeference/georeferencer.py | 21 ++++---- tests/test_warp.py | 77 ++++++++++++++++++++++++++++++ 2 files changed, 90 insertions(+), 8 deletions(-) create mode 100644 tests/test_warp.py diff --git a/ohmg/georeference/georeferencer.py b/ohmg/georeference/georeferencer.py index e3ebaa58..9864ce7f 100644 --- a/ohmg/georeference/georeferencer.py +++ b/ohmg/georeference/georeferencer.py @@ -5,7 +5,7 @@ import time from io import StringIO from pathlib import Path -from typing import List, Tuple +from typing import List, Tuple, Union from uuid import uuid4 from django.conf import settings @@ -249,7 +249,10 @@ def _geo_coords_from_gcp(self, gcp: gdal.GCP) -> Tuple[float, float]: return (gcp.GCPX, gcp.GCPY) def _pixel_coords_from_gcp( - self, gcp: gdal.GCP, cartesian_y: bool = False + self, + gcp: gdal.GCP, + cartesian_y: bool = False, + img_height: Union[float, None] = None, ) -> Tuple[float, float]: """Return the image pixel coords from the input GCP. @@ -260,8 +263,10 @@ def _pixel_coords_from_gcp( image.""" x, y = gcp.GCPPixel, gcp.GCPLine if cartesian_y: - ds = gdal.Open(self.gcps_vrt.get_vsi_url()) - y = ds.RasterYSize - y + if img_height is None: + ds = gdal.Open(self.gcps_vrt.get_vsi_url()) + img_height = ds.RasterYSize + y = img_height - y return (x, y) def _calculate_scale(self) -> float: @@ -288,7 +293,7 @@ def _calculate_scale(self) -> float: return pt_dist / px_dist - def _calculate_rotation_from_north(self) -> float: + def _calculate_azimuth(self, img_height: Union[float, None] = None) -> float: """Compares two GCPs and calculates the difference in the angles between the geometric points and the pixel points. @@ -306,8 +311,8 @@ def _calculate_rotation_from_north(self) -> float: self._geo_coords_from_gcp(upper_gcp), self._geo_coords_from_gcp(lower_gcp) ) px_degrees = angle_from_coords( - self._pixel_coords_from_gcp(upper_gcp, cartesian_y=True), - self._pixel_coords_from_gcp(lower_gcp, cartesian_y=True), + self._pixel_coords_from_gcp(upper_gcp, cartesian_y=True, img_height=img_height), + self._pixel_coords_from_gcp(lower_gcp, cartesian_y=True, img_height=img_height), ) difference = pt_degrees - px_degrees @@ -421,7 +426,7 @@ def make_warped_vrt( scale = self._calculate_scale() ## get rotation - rotation = self._calculate_rotation_from_north() + rotation = self._calculate_azimuth() ## adjust and convert to arcseconds arcseconds = (rotation + 90) * 3600 diff --git a/tests/test_warp.py b/tests/test_warp.py new file mode 100644 index 00000000..5b38907a --- /dev/null +++ b/tests/test_warp.py @@ -0,0 +1,77 @@ +import math + +from django.test import tag +from osgeo import gdal + +from ohmg.georeference.georeferencer import Georeferencer + +from .base import OHMGTestCase + + +@tag("warp") +class HelmertTransformationTestCase(OHMGTestCase): + """These tests assume an input image with width 5 and height 8. + Two GCPs are used, which are placed vertical from each other + near the left edge of the image.""" + + img_width = 5 + img_height = 8 + + ## Create 2 GCPs for use in each test. These original + ## coordinates will place the image with 0 rotation and .5 scale + ## in the initial test. Later cases will change the X and Y coords + ## of each GCP to test against different angles, etc. + + ## 8 different permutations are tested altogether, the + + # note 1. GCP args are: geo x, geo y, geo z (not used), img x, img y + # note 2: img y uses inverse y axis (per GCP convention) + # note 3: The distance between GCPs is 5 which creates a 3,4,5 + # triangle during rotated permutations of the test (helpful) + gcp1 = gdal.GCP(0, 2.5, 0, 1, 1) + gcp2 = gdal.GCP(0, 0, 0, 1, 6) + + # this is the small angle of a 3/4/5 triangle + + def test_helmert_transformation_A(self): + """This test was hand-created to place the mock 5x8 image + into generic cartesian coordinate space such that: + rotation = 0 + scale = .5 + x, y offset (top left corner) = 4.5, 7.5. + + All other tests were generated by modifying GCP inputs to this one + and recording the output. + """ + + theta_345 = math.degrees(math.asin(3 / 5)) + + ## this is a list of variables that define each position of the test + ## the adjustment is applied to GCP1, which starts with geo position 0, 2.5 + ## x adjustment, y adjustment, target azimuth + data_matrix = [ + (0, 2.5, 0, -0.5, 3), + (2, 1.5, 90 - theta_345, 55, 55), + (2.5, 0, 90, 55, 55), + (2, -1.5, 90 + theta_345, 55, 55), + (0, -2.5, 180, 55, 55), + (-2, -1.5, 270 - theta_345, 55, 55), + (-2.5, 0, 270, 55, 55), + (-2, 1.5, 270 + theta_345, 55, 55), + ] + + for gcpx, gcpy, target_azimuth, x_offset, y_offset in data_matrix: + gcp2 = gdal.GCP(0, 0, 0, 1, 6) + gcp1 = gdal.GCP(gcpx, gcpy, 0, 1, 1) + + g = Georeferencer(crs="EPSG:3857", transformation="helmert", gcps_gdal=[gcp1, gcp2]) + + scale = g._calculate_scale() + self.assertEqual(scale, 0.5) + + azimuth = g._calculate_azimuth(img_height=self.img_height) + self.assertEqual(azimuth, target_azimuth) + + dx, dy = g._calculate_helmert_offsets(scale=scale, rotation=azimuth) + self.assertAlmostEqual(dx, x_offset) + self.assertAlmostEqual(dy, y_offset) From 421ed87e53890b115adbc391481b29318131d2c4 Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Fri, 19 Jun 2026 12:38:22 -0500 Subject: [PATCH 2/5] clarify rotation calc, handle >2 GCPs --- ohmg/georeference/geometry.py | 69 +++++++++++++++++------ ohmg/georeference/georeferencer.py | 89 ++++++++++++------------------ 2 files changed, 86 insertions(+), 72 deletions(-) diff --git a/ohmg/georeference/geometry.py b/ohmg/georeference/geometry.py index cc6471a3..9c7f96dc 100644 --- a/ohmg/georeference/geometry.py +++ b/ohmg/georeference/geometry.py @@ -1,25 +1,60 @@ import math -from typing import Tuple +from typing import List, Tuple +import numpy as np from django.contrib.gis.geos import LineString -def angle_from_coords(pt1: Tuple[float, float], pt2: Tuple[float, float]) -> float: - """ - Calculates the absolute Cartesian angle (in degrees) of the vector - pointing from p1 to p2, relative to the positive Y-axis. - """ - - # Calculate differences - dx = pt2[0] - pt1[0] - dy = pt2[1] - pt1[1] - - # Calculate angle in radians (-pi to pi) - radians = math.atan2(dx, dy) - - degrees = math.degrees(radians) - - return degrees +def azimuth_from_coords(coords: List[Tuple[float, float]]) -> float: + """Calculate the least-squares slope of the input coordinates, return + as degrees relative to the x axis.""" + + x_coords = [i[0] for i in coords] + y_coords = [i[1] for i in coords] + x_diff = x_coords[-1] - x_coords[0] + y_diff = y_coords[-1] - y_coords[0] + + # handle cases where the fit line would be horizontal or vertical, + # as well as an issue where passing all 0s to np.polyfit (e.g. all + # of the x values are 0, even though there are different y values) + # raises an exception + if len(set(x_coords)) == 1: + if y_diff > 0: + azimuth = 0 + else: + azimuth = 180 + elif len(set(y_coords)) == 1: + if x_diff > 0: + azimuth = 90 + else: + azimuth = 270 + # now handle all other cases by calculating the slope + else: + slope, intercept = np.polyfit(x_coords, y_coords, 1) + + # this is the angle from 0 axis + angle = math.degrees(math.atan(slope)) + + # now convert the angle to degrees from north by comparing + # the first and last set of coords to determine the general + # orientation of the fit line + x_diff = coords[-1][0] - coords[0][0] + y_diff = coords[-1][1] - coords[0][1] + + # orientation: ne + if x_diff > 0 and y_diff > 0: + azimuth = 90 - angle + # orientation: se + elif x_diff > 0 and y_diff < 0: + azimuth = 90 + abs(angle) + # orientation: sw + elif x_diff < 0 and y_diff < 0: + azimuth = 270 - angle + # orientation: nw + elif x_diff < 0 and y_diff > 0: + azimuth = 270 + abs(angle) + + return azimuth def extend_vector( diff --git a/ohmg/georeference/georeferencer.py b/ohmg/georeference/georeferencer.py index 9864ce7f..ebe68470 100644 --- a/ohmg/georeference/georeferencer.py +++ b/ohmg/georeference/georeferencer.py @@ -13,7 +13,7 @@ from ohmg.core.utils.srs import retrieve_srs_wkt -from .geometry import angle_from_coords +from .geometry import azimuth_from_coords logger = logging.getLogger(__name__) @@ -244,78 +244,57 @@ def _load_gcps_from_geojson(self, geo_json): ) self.gcps.append(gcp) - def _geo_coords_from_gcp(self, gcp: gdal.GCP) -> Tuple[float, float]: - """Return the geographic x, y coords from the input GCP""" - return (gcp.GCPX, gcp.GCPY) - - def _pixel_coords_from_gcp( - self, - gcp: gdal.GCP, - cartesian_y: bool = False, - img_height: Union[float, None] = None, - ) -> Tuple[float, float]: - """Return the image pixel coords from the input GCP. - - Y coordinate is measured down from the top of the image. - - If cartesian_y=True, then invert the Y coordinate against the height - of the dataset, to match a cartesian plane with 0,0 at the bottom left of the - image.""" - x, y = gcp.GCPPixel, gcp.GCPLine - if cartesian_y: - if img_height is None: - ds = gdal.Open(self.gcps_vrt.get_vsi_url()) - img_height = ds.RasterYSize - y = img_height - y - return (x, y) - def _calculate_scale(self) -> float: """Compares two GCPs and returns a scale factor.""" - if len(self.gcps) != 2: - raise Exception("Two GCPs are needed to calculate a scale factor") + if len(self.gcps) < 2: + raise Exception("At least GCPs are needed to calculate a scale factor") - gcp1, gcp2 = self.gcps + gcp1, gcp2 = self.gcps[0], self.gcps[1] # distance between the geographic coords in each GCP # this distance is absolute so the order of the coords doesn't matter - pt_dist = math.dist( - self._geo_coords_from_gcp(gcp1), - self._geo_coords_from_gcp(gcp2), + geo_dist = math.dist( + (gcp1.GCPX, gcp1.GCPY), + (gcp2.GCPX, gcp2.GCPY), ) # distance between the pixel coords in each GCP # this distance is absolute so the order of the coords doesn't matter - px_dist = math.dist( - self._pixel_coords_from_gcp(gcp1), - self._pixel_coords_from_gcp(gcp2), + img_dist = math.dist( + (gcp1.GCPPixel, gcp1.GCPLine), + (gcp2.GCPPixel, gcp2.GCPLine), ) - return pt_dist / px_dist + return geo_dist / img_dist - def _calculate_azimuth(self, img_height: Union[float, None] = None) -> float: + def _calculate_rotation(self, img_height: Union[float, None] = None) -> float: """Compares two GCPs and calculates the difference in the angles between the geometric points and the pixel points. Returns the angle in degrees from 'north', i.e. the positive Y axis""" - if len(self.gcps) != 2: - raise Exception("Two GCPs are needed to calculate a scale factor") - ## the GCPLine value is the number of pixels DOWN from the top of the - ## image, so we'll call it "upper" because it appears the higher of the - ## two on the page - upper_gcp = min(self.gcps, key=lambda x: x.GCPLine) - lower_gcp = max(self.gcps, key=lambda x: x.GCPLine) + if len(self.gcps) < 2: + raise Exception("At least two GCPs are needed to calculate a scale factor") - pt_degrees = angle_from_coords( - self._geo_coords_from_gcp(upper_gcp), self._geo_coords_from_gcp(lower_gcp) - ) - px_degrees = angle_from_coords( - self._pixel_coords_from_gcp(upper_gcp, cartesian_y=True, img_height=img_height), - self._pixel_coords_from_gcp(lower_gcp, cartesian_y=True, img_height=img_height), - ) + # sort the GCPs so they are ordered from lowest to highest, + # then right to left, as they appear on the source image. + # this allows us to figure out the orientation + self.gcps.sort(key=lambda x: x.GCPPixel) + self.gcps.sort(key=lambda x: x.GCPLine, reverse=True) + + # make sure img height is set because it is needed to properly + # handle the inverted Y coords. + if not img_height: + ds = gdal.Open(self.gcps_vrt.get_vsi_url()) + img_height = ds.RasterYSize + img_coords = [(i.GCPPixel, img_height - i.GCPLine) for i in self.gcps] + img_azimuth = azimuth_from_coords(img_coords) + + geo_coords = [(i.GCPX, i.GCPY) for i in self.gcps] + geo_azimuth = azimuth_from_coords(geo_coords) - difference = pt_degrees - px_degrees + difference = geo_azimuth - img_azimuth return difference def _calculate_helmert_offsets(self, scale: float, rotation: float) -> Tuple[float, float]: @@ -330,8 +309,8 @@ def _calculate_helmert_offsets(self, scale: float, rotation: float) -> Tuple[flo dist_to_page_edge = use_gcp.GCPLine * scale dist_to_page_top = use_gcp.GCPPixel * scale - ## rotation is degrees from positive y-axis, adjust to be degrees from - ## positive x-axis + ## rotation is azimuth, i.e. degrees from positive y-axis, + ## must adjust to be degrees from positive x-axis theta1 = 90 - rotation ## further adjustments to normalize if theta1 < 0: @@ -426,7 +405,7 @@ def make_warped_vrt( scale = self._calculate_scale() ## get rotation - rotation = self._calculate_azimuth() + rotation = self._calculate_rotation() ## adjust and convert to arcseconds arcseconds = (rotation + 90) * 3600 From 3a6c3f31ca1595f5463d8cc5687092bf943667da Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Fri, 19 Jun 2026 13:02:43 -0500 Subject: [PATCH 3/5] update helmert transformation tests --- tests/test_warp.py | 80 ++++++++++++++++++---------------------------- 1 file changed, 31 insertions(+), 49 deletions(-) diff --git a/tests/test_warp.py b/tests/test_warp.py index 5b38907a..16348888 100644 --- a/tests/test_warp.py +++ b/tests/test_warp.py @@ -10,68 +10,50 @@ @tag("warp") class HelmertTransformationTestCase(OHMGTestCase): - """These tests assume an input image with width 5 and height 8. - Two GCPs are used, which are placed vertical from each other - near the left edge of the image.""" - - img_width = 5 - img_height = 8 - - ## Create 2 GCPs for use in each test. These original - ## coordinates will place the image with 0 rotation and .5 scale - ## in the initial test. Later cases will change the X and Y coords - ## of each GCP to test against different angles, etc. - - ## 8 different permutations are tested altogether, the - - # note 1. GCP args are: geo x, geo y, geo z (not used), img x, img y - # note 2: img y uses inverse y axis (per GCP convention) - # note 3: The distance between GCPs is 5 which creates a 3,4,5 - # triangle during rotated permutations of the test (helpful) - gcp1 = gdal.GCP(0, 2.5, 0, 1, 1) - gcp2 = gdal.GCP(0, 0, 0, 1, 6) - - # this is the small angle of a 3/4/5 triangle - - def test_helmert_transformation_A(self): - """This test was hand-created to place the mock 5x8 image - into generic cartesian coordinate space such that: - rotation = 0 - scale = .5 - x, y offset (top left corner) = 4.5, 7.5. - - All other tests were generated by modifying GCP inputs to this one - and recording the output. + def test_helmert_transformation_calculations(self): + """This test runs through 8 permutations of GCPs, and makes + sure that the calculations used to set up the helmert + transformations return the right values for every permutation. """ + # assume an image with these dimensions + img_width, img_height = 5, 8 # noqa: F841 + + # inner (smallest) angle for a 3,4,5 triangle theta_345 = math.degrees(math.asin(3 / 5)) - ## this is a list of variables that define each position of the test - ## the adjustment is applied to GCP1, which starts with geo position 0, 2.5 - ## x adjustment, y adjustment, target azimuth + # variables that define positions and expected values to test data_matrix = [ (0, 2.5, 0, -0.5, 3), - (2, 1.5, 90 - theta_345, 55, 55), - (2.5, 0, 90, 55, 55), - (2, -1.5, 90 + theta_345, 55, 55), - (0, -2.5, 180, 55, 55), - (-2, -1.5, 270 - theta_345, 55, 55), - (-2.5, 0, 270, 55, 55), - (-2, 1.5, 270 + theta_345, 55, 55), + (2, 1.5, 90 - theta_345, 2.1, 2.2), + (2.5, 0, 90, 3, 0.5), + (2, -1.5, 90 + theta_345, 2.7, -1.4), + (0, -2.5, 180, 0.5, -3), + (-2, -1.5, 270 - theta_345, -2.1, -2.2), + (-2.5, 0, 270, -3, -0.5), + (-2, 1.5, 270 + theta_345, -2.7, 1.4), ] - for gcpx, gcpy, target_azimuth, x_offset, y_offset in data_matrix: - gcp2 = gdal.GCP(0, 0, 0, 1, 6) - gcp1 = gdal.GCP(gcpx, gcpy, 0, 1, 1) + # GCP 1 stays constant + # GCP 2's image coords stay constant while the geo coords move + # clockwise around the origin. + gcp1 = gdal.GCP(0, 0, 0, 1, 6) + for gcpx, gcpy, target_rotation, x_offset, y_offset in data_matrix: + # note 1. GCP args are: geo x, geo y, geo z (not used), img x, img y + # note 2: img y uses inverse y axis (per GCP convention) + # note 3: The distance between img GCPs is 5 which creates a 3,4,5 + # triangle during rotated permutations of the test (helpful) + # note 4: All geo GCP coords halve the dimensions of the triangle + # so scale is .5 + gcp2 = gdal.GCP(gcpx, gcpy, 0, 1, 1) g = Georeferencer(crs="EPSG:3857", transformation="helmert", gcps_gdal=[gcp1, gcp2]) scale = g._calculate_scale() self.assertEqual(scale, 0.5) + rotation = g._calculate_rotation(img_height=img_height) + self.assertEqual(rotation, target_rotation) - azimuth = g._calculate_azimuth(img_height=self.img_height) - self.assertEqual(azimuth, target_azimuth) - - dx, dy = g._calculate_helmert_offsets(scale=scale, rotation=azimuth) + dx, dy = g._calculate_helmert_offsets(scale=scale, rotation=rotation) self.assertAlmostEqual(dx, x_offset) self.assertAlmostEqual(dy, y_offset) From 0ad3e382ac8f4084644da5e3c257712f235a54b7 Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Fri, 19 Jun 2026 22:50:44 +0000 Subject: [PATCH 4/5] upgrade gdal python bindings to 3.8 --- pyproject.toml | 2 +- uv.lock | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index d9b78972..5bd0b9f6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,7 +32,7 @@ dependencies = [ "dialogos==0.4", "pinax==0.9a2", "pinax-announcements==4.0.1", - "gdal>=3.5,<3.6", + "gdal>=3.8,<3.8.5", "Pillow<10.0.0", "django_compressor", "natsort", diff --git a/uv.lock b/uv.lock index 3a4cb75b..50af47dd 100644 --- a/uv.lock +++ b/uv.lock @@ -581,9 +581,9 @@ wheels = [ [[package]] name = "gdal" -version = "3.5.3" +version = "3.8.4" source = { registry = "https://pypi.org/simple" } -sdist = { url = "https://files.pythonhosted.org/packages/8b/9f/c9fe38625e8d9ea814d4c3935773341cffb226b12ce2356bad80a9c0e4e8/GDAL-3.5.3.tar.gz", hash = "sha256:ff11ff32b400086930e0e0d15e35da598dfae548423512759e111ea9ad65b19b", size = 756836, upload-time = "2022-10-31T16:40:05.87Z" } +sdist = { url = "https://files.pythonhosted.org/packages/21/6f/7f76d6ff3673fa18e9edd8d76679968d682f937cded93311f82fa172697f/GDAL-3.8.4.tar.gz", hash = "sha256:7c51e0ae7a7ccf43ad9e4bf435176baa9276653dfa16fd167c3632f6e7275207", size = 802459, upload-time = "2024-02-18T15:35:29.047Z" } [[package]] name = "h11" @@ -863,7 +863,7 @@ requires-dist = [ { name = "django-markdownx", specifier = "==4.0.7" }, { name = "django-ninja", specifier = "==0.21.0" }, { name = "django-storages", specifier = "==1.13.2" }, - { name = "gdal", specifier = ">=3.5,<3.6" }, + { name = "gdal", specifier = ">=3.8,<3.8.5" }, { name = "humanize" }, { name = "kombu", specifier = "==5.6.0" }, { name = "lxml" }, From c128ad79ba33ad291d61bf4885104c3d59820b36 Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Sat, 20 Jun 2026 20:55:52 +0000 Subject: [PATCH 5/5] pin gdal python bindings to 3.8.4 --- pyproject.toml | 2 +- uv.lock | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 5bd0b9f6..f3fb89cb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,7 +32,7 @@ dependencies = [ "dialogos==0.4", "pinax==0.9a2", "pinax-announcements==4.0.1", - "gdal>=3.8,<3.8.5", + "gdal==3.8.4", "Pillow<10.0.0", "django_compressor", "natsort", diff --git a/uv.lock b/uv.lock index 50af47dd..d960e9bf 100644 --- a/uv.lock +++ b/uv.lock @@ -863,7 +863,7 @@ requires-dist = [ { name = "django-markdownx", specifier = "==4.0.7" }, { name = "django-ninja", specifier = "==0.21.0" }, { name = "django-storages", specifier = "==1.13.2" }, - { name = "gdal", specifier = ">=3.8,<3.8.5" }, + { name = "gdal", specifier = "==3.8.4" }, { name = "humanize" }, { name = "kombu", specifier = "==5.6.0" }, { name = "lxml" },