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 e3ebaa58..ebe68470 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 @@ -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,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 diff --git a/pyproject.toml b/pyproject.toml index d9b78972..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.5,<3.6", + "gdal==3.8.4", "Pillow<10.0.0", "django_compressor", "natsort", diff --git a/tests/test_warp.py b/tests/test_warp.py new file mode 100644 index 00000000..16348888 --- /dev/null +++ b/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) diff --git a/uv.lock b/uv.lock index 3a4cb75b..d960e9bf 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.4" }, { name = "humanize" }, { name = "kombu", specifier = "==5.6.0" }, { name = "lxml" },