Helmert tests

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(

ohmg/georeference/georeferencer.py

@@ -5,15 +5,15 @@
 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
 from osgeo import gdal, ogr, osr
 
 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,73 +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
-    ) -> 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:
-            ds = gdal.Open(self.gcps_vrt.get_vsi_url())
-            y = ds.RasterYSize - 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_rotation_from_north(self) -> 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),
-            self._pixel_coords_from_gcp(lower_gcp, cartesian_y=True),
-        )
+        # 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]:
@@ -325,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:
@@ -421,7 +405,7 @@ def make_warped_vrt(
             scale = self._calculate_scale()
 
             ## get rotation
-            rotation = self._calculate_rotation_from_north()
+            rotation = self._calculate_rotation()
             ## adjust and convert to arcseconds
             arcseconds = (rotation + 90) * 3600
 

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.4",
     "Pillow<10.0.0",
     "django_compressor",
     "natsort",

tests/test_warp.py

@@ -0,0 +1,59 @@
+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):
+    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))
+
+        # variables that define positions and expected values to test
+        data_matrix = [
+            (0, 2.5, 0, -0.5, 3),
+            (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),
+        ]
+
+        # 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)
+
+            dx, dy = g._calculate_helmert_offsets(scale=scale, rotation=rotation)
+            self.assertAlmostEqual(dx, x_offset)
+            self.assertAlmostEqual(dy, y_offset)