From 8ca1447f9939a1f4f937cd556589e487f825e6a3 Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Wed, 3 Jun 2026 16:56:22 -0500 Subject: [PATCH 01/10] move some geometry utils to new file --- ohmg/georeference/geometry.py | 63 +++++++++++++++++++++++++++++++++++ ohmg/georeference/splitter.py | 40 ++-------------------- 2 files changed, 66 insertions(+), 37 deletions(-) create mode 100644 ohmg/georeference/geometry.py diff --git a/ohmg/georeference/geometry.py b/ohmg/georeference/geometry.py new file mode 100644 index 00000000..37ae3d50 --- /dev/null +++ b/ohmg/georeference/geometry.py @@ -0,0 +1,63 @@ +import math +from typing import Tuple + +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 X-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 extend_vector( + p1: Tuple[float, float], + p2: Tuple[float, float], + distance: float, +) -> Tuple[float, float]: + """https://math.stackexchange.com/a/3346108 (credit to Oliver Roche) + takes the two input points, which represent a vector, and creates a + third point that would extend that vector by the given distance.""" + + x1, y1 = p1 + x2, y2 = p2 + rise = y2 - y1 + run = x2 - x1 + + norm = math.sqrt((run**2) + (rise**2)) + + # if negative coords are used norm will be 0.0, silently return original point + if norm == 0.0: + return (x2, y2) + + x3 = x2 + distance * (run / norm) + y3 = y2 + distance * (rise / norm) + + return (x3, y3) + + +def extend_linestring(linestring: LineString, distance: int = 10) -> LineString: + """takes the input GEOS LineString and extends it in both directions + (following the trajectory of each end segment) by the given distance.""" + + coord_list = list(linestring.coords) + + new_start = extend_vector(coord_list[1], coord_list[0], distance) + new_end = extend_vector(coord_list[-2], coord_list[-1], distance) + + coord_list.insert(0, new_start) + coord_list.append(new_end) + + return LineString(coord_list) diff --git a/ohmg/georeference/splitter.py b/ohmg/georeference/splitter.py index 81854551..37785d09 100644 --- a/ohmg/georeference/splitter.py +++ b/ohmg/georeference/splitter.py @@ -1,5 +1,4 @@ import logging -import math import os import time @@ -9,6 +8,8 @@ from django.db.models import FileField from PIL import Image, ImageDraw, ImageFilter +from .geometry import extend_linestring + logger = logging.getLogger(__name__) @@ -39,41 +40,6 @@ def transform_coordinates(self, shape, img_height): coords.append((x, img_height - y)) return coords - def extend_vector(self, p1, p2, distance): - """https://math.stackexchange.com/a/3346108 (credit to Oliver Roche) - takes the two input points, which represent a vector, and creates a - third point that would extend that vector by the given distance.""" - - x1, y1 = p1 - x2, y2 = p2 - rise = y2 - y1 - run = x2 - x1 - - norm = math.sqrt((run**2) + (rise**2)) - - # if negative coords are used norm will be 0.0, silently return original point - if norm == 0.0: - return (x2, y2) - - x3 = x2 + distance * (run / norm) - y3 = y2 + distance * (rise / norm) - - return (x3, y3) - - def extend_linestring(self, linestring, distance=10): - """takes the input GEOS LineString and extends it in both directions - (following the trajectory of each end segment) by the given distance.""" - - coord_list = list(linestring.coords) - - new_start = self.extend_vector(coord_list[1], coord_list[0], distance) - new_end = self.extend_vector(coord_list[-2], coord_list[-1], distance) - - coord_list.insert(0, new_start) - coord_list.append(new_end) - - return LineString(coord_list) - def generate_divisions(self, cutlines): """takes the input border and then tries to cut it with the cutlines. any sub polygons resulting from the cut are also compared to the cutlines, @@ -86,7 +52,7 @@ def generate_divisions(self, cutlines): for line in cutlines: ## this function extends each end of the original line by 10 pixels. ## this facilitates a more robust splitting process. - ls_extended = self.extend_linestring(LineString(line)) + ls_extended = extend_linestring(LineString(line)) cut_shapes.append({"geom": ls_extended, "used": False}) ## candidates is a list of polygons that may be the final polygons From 85321666e0cb4aede288a5728cfa27004f32f9cc Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Wed, 3 Jun 2026 16:58:52 -0500 Subject: [PATCH 02/10] initial addition of helmert transformation logic (in-progress) --- ohmg/georeference/georeferencer.py | 210 ++++++++++++++++++++++++----- 1 file changed, 180 insertions(+), 30 deletions(-) diff --git a/ohmg/georeference/georeferencer.py b/ohmg/georeference/georeferencer.py index f15c7952..20e7316a 100644 --- a/ohmg/georeference/georeferencer.py +++ b/ohmg/georeference/georeferencer.py @@ -1,9 +1,11 @@ import logging +import math import os import sys import time from io import StringIO from pathlib import Path +from typing import List, Tuple from uuid import uuid4 from django.conf import settings @@ -11,11 +13,19 @@ from ohmg.core.utils.srs import retrieve_srs_wkt +from .geometry import angle_from_coords + logger = logging.getLogger(__name__) gdal.SetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS") TRANSFORMATION_LOOKUP = { + "helmert": { + "id": "helmert", + "gdal_code": None, + "name": "Helmert", + "desc": "implemented as four-parameter transformation", + }, "tps": { "id": "tps", "gdal_code": -1, @@ -139,7 +149,7 @@ def __init__( crs="EPSG:3857", transformation="poly1", # three different ways to add GCPs, one must be provided - gcps_gdal=None, + gcps_gdal: List[gdal.GCP] = [], gcps_geojson=None, gcps_points_file=None, verbose=False, @@ -164,7 +174,8 @@ def __init__( # handle the input GCPs to GDAL GCPs, method depends on the input # format. self.gcps should be a list of gdal.GCP objects. - if gcps_gdal: + self.gcps = [] + if len(gcps_gdal) > 0: self.gcps = gcps_gdal elif gcps_geojson: self._load_gcps_from_geojson(gcps_geojson) @@ -233,6 +244,111 @@ 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) -> Tuple[float, float]: + """Return the image pixel coords from the input GCP""" + return (gcp.GCPPixel, gcp.GCPLine) + + def _calculate_scale(self) -> float: + """Compares two GCPs and returns a scale factor.""" + print("CALCULATE SCALE") + if len(self.gcps) < 2: + raise Exception("Two GCPs are needed to calculate a scale factor") + + gcp1_pt = self._geo_coords_from_gcp(self.gcps[0]) + gcp2_pt = self._geo_coords_from_gcp(self.gcps[1]) + gcp1_px = self._pixel_coords_from_gcp(self.gcps[0]) + gcp2_px = self._pixel_coords_from_gcp(self.gcps[1]) + + # distance between the geographic coords in each GCP + pt_dist = math.dist(gcp1_pt, gcp2_pt) + + # distance between the pixel coords in each GCP + px_dist = math.dist(gcp1_px, gcp2_px) + + return pt_dist / px_dist + + def _calculate_rotation(self) -> float: + """Compares two GCPs and calculates the difference in the angles + between the geometric points and the pixel points. + + Returns the angle in degrees""" + print("CALCULATE ROTATION") + if len(self.gcps) < 2: + raise Exception("Two GCPs are needed to calculate a scale factor") + + gcp1_pt = self._geo_coords_from_gcp(self.gcps[0]) + gcp2_pt = self._geo_coords_from_gcp(self.gcps[1]) + gcp1_px = self._pixel_coords_from_gcp(self.gcps[0]) + gcp2_px = self._pixel_coords_from_gcp(self.gcps[1]) + + print("gcp1_pt", gcp1_pt) + print("gcp2_pt", gcp2_pt) + print("gcp1_px", gcp1_px) + print("gcp2_px", gcp2_px) + + pt_degrees = angle_from_coords(gcp1_pt, gcp2_pt) + px_degrees = -angle_from_coords(gcp1_px, gcp2_px) + + print("pt_degrees", pt_degrees) + print("px_degrees", px_degrees) + + difference = pt_degrees - px_degrees + print("difference", difference) + return difference + + def _calculate_offsets(self, scale: float, rotation: float) -> Tuple[float, float]: + """Calculates the x and y offsets needed for helmert transformation, based + on GCP + """ + print("CALCULATE OFFSETS") + gcp1 = self.gcps[0] + print("offset gcp", gcp1.GCPLine, gcp1.GCPPixel) + + geo_dist_to_page_edge = gcp1.GCPLine * scale + geo_dist_to_page_top = gcp1.GCPPixel * scale + geo_dist_to_page_corner = math.sqrt(geo_dist_to_page_edge**2 + geo_dist_to_page_top**2) + + print("geo_dist_to_page_edge", geo_dist_to_page_edge) + print("geo_dist_to_page_top", geo_dist_to_page_top) + print("geo_dist_to_page_corner", geo_dist_to_page_corner) + + print("rotation", rotation) + print("rotation used", rotation + 270) + theta1 = 90 - rotation + print("theta1", theta1) + if theta1 < 0: + theta1 += 180 + if theta1 > 180: + theta1 -= 180 + print("theta1", theta1) + + m_degrees = 180 - (90 + theta1) + m_radians = math.radians(m_degrees) + n_radians = math.atan(geo_dist_to_page_edge / geo_dist_to_page_top) + theta3 = n_radians + m_radians + + print("m_radians", m_radians) + print("n_radians", n_radians) + print("theta3", theta3) + print("m_degrees", math.degrees(m_radians)) + print("n_degrees", math.degrees(n_radians)) + print("theta3_degrees", math.degrees(theta3)) + + a2 = math.cos(theta3) * geo_dist_to_page_corner + o2 = math.sin(theta3) * geo_dist_to_page_corner + + print("a2", a2) + print("o2", o2) + + dx = gcp1.GCPX - a2 + dy = gcp1.GCPY + o2 + + return (dx, dy) + def cleanup_files(self): for vrt in [ self.gcps_vrt, @@ -288,34 +404,68 @@ def make_warped_vrt( self.make_gcps_vrt(src_path, out_name) self.warped_vrt = VRTHandler(out_name, as_variant="modified") - wo = gdal.WarpOptions( - creationOptions=[ - # "NUM_THREADS=ALL_CPUS", - # ## originally used this set of flags used - # # "COMPRESS=DEFLATE", - # ## should have been used PREDICTOR=2 with DEFLATE but didn't know about it - # # "PREDICTOR=2" - # ## useful in general but not needed when using COG driver - "BLOCKXSIZE=512", - "BLOCKYSIZE=512", - # ## advisable if using JPEG with GTiff, but not supported in COG - # # "JPEG_QUALITY=75", - # # "PHOTOMETRIC=YCBCR", - # ## Use JPEG, as recommended by Paul Ramsey article: - # ## https://blog.cleverelephant.ca/2015/02/geotiff-compression-for-dummies.html - # "COMPRESS=JPEG", - ], - transformerOptions=[ - f"DST_SRS={self.crs_wkt}", - f'MAX_GCP_ORDER={self.transformation["gdal_code"]}', - ], - format="VRT", - dstSRS=f"{self.crs_code}", - # srcNodata=src_nodata, - # srcAlpha=True, - dstAlpha=True, - resampleAlg="nearest", - ) + if len(self.gcps) == 2 and self.transformation["id"] == "helmert": + ## get scale factor + scale = self._calculate_scale() + + ## get rotation + rotation_degrees = self._calculate_rotation() + ## adjust to north = 0 and convert to arcseconds + arcseconds = (rotation_degrees + 270) * 3600 + + ## get the x y offsets + dx, dy = self._calculate_offsets(scale, rotation_degrees) + + pipeline = ( + "+proj=pipeline " + "+step +proj=axisswap +order=2,1 " + f"+step +proj=helmert +x={dx} +y={dy} +theta={arcseconds} +s={scale}" + ) + + print(pipeline) + + wo = gdal.WarpOptions( + creationOptions=[ + "BLOCKXSIZE=512", + "BLOCKYSIZE=512", + ], + coordinateOperation=pipeline, + format="VRT", + dstSRS="EPSG:3857", + transformerOptions=["SRC_METHOD=NO_GEOTRANSFORM"], + dstAlpha=True, + resampleAlg="nearest", + ) + else: + wo = gdal.WarpOptions( + creationOptions=[ + # "NUM_THREADS=ALL_CPUS", + # ## originally used this set of flags used + # # "COMPRESS=DEFLATE", + # ## should have been used PREDICTOR=2 with DEFLATE but didn't know about it + # # "PREDICTOR=2" + # ## useful in general but not needed when using COG driver + "BLOCKXSIZE=512", + "BLOCKYSIZE=512", + # ## advisable if using JPEG with GTiff, but not supported in COG + # # "JPEG_QUALITY=75", + # # "PHOTOMETRIC=YCBCR", + # ## Use JPEG, as recommended by Paul Ramsey article: + # ## https://blog.cleverelephant.ca/2015/02/geotiff-compression-for-dummies.html + # "COMPRESS=JPEG", + ], + transformerOptions=[ + f"DST_SRS={self.crs_wkt}", + f'MAX_GCP_ORDER={self.transformation["gdal_code"]}', + ], + format="VRT", + dstSRS=f"{self.crs_code}", + # srcNodata=src_nodata, + # srcAlpha=True, + dstAlpha=True, + resampleAlg="nearest", + ) + try: gdal.Warp(str(self.warped_vrt.get_path()), self.gcps_vrt.get_vsi_url(), options=wo) except Exception as e: From 2fb1bee3d16bf3d94df3ebd80aa82cff88cd3ab0 Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Wed, 3 Jun 2026 16:59:13 -0500 Subject: [PATCH 03/10] add generic confirmation modal --- .../src/components/modals/ConfirmModal.svelte | 32 +++++++++++++++++++ 1 file changed, 32 insertions(+) create mode 100644 ohmg/frontend/svelte_components/src/components/modals/ConfirmModal.svelte diff --git a/ohmg/frontend/svelte_components/src/components/modals/ConfirmModal.svelte b/ohmg/frontend/svelte_components/src/components/modals/ConfirmModal.svelte new file mode 100644 index 00000000..d130c282 --- /dev/null +++ b/ohmg/frontend/svelte_components/src/components/modals/ConfirmModal.svelte @@ -0,0 +1,32 @@ + + + + +

"Confirm this action?"

+ + + + From baddb60358db81940ecdb531b9cc867a855adfde Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Wed, 3 Jun 2026 17:00:00 -0500 Subject: [PATCH 04/10] allow 2 GCP submit after prompt --- .../src/components/Georeferencer.svelte | 26 +++++++++++++++---- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/ohmg/frontend/svelte_components/src/components/Georeferencer.svelte b/ohmg/frontend/svelte_components/src/components/Georeferencer.svelte index ff7e680b..f323fc1b 100644 --- a/ohmg/frontend/svelte_components/src/components/Georeferencer.svelte +++ b/ohmg/frontend/svelte_components/src/components/Georeferencer.svelte @@ -54,6 +54,8 @@ import InfoModalButton from './buttons/InfoModalButton.svelte'; import SigninReminder from './common/SigninReminder.svelte'; + import ConfirmModal from './modals/ConfirmModal.svelte'; + export let CONTEXT; export let REGION; export let MAP; @@ -119,7 +121,7 @@ leaveOkay = false; enableButtons = true; } - $: enableSave = gcpList.length >= 3 && enableButtons; + $: enableSave = gcpList.length >= 2 && enableButtons; let countdown = 10; let timer; @@ -150,9 +152,11 @@ const noteInputElId = 'note-input'; let currentTransformation = 'poly1'; + $: currentTransformation = gcpList.length == 2 ? 'helmert' : currentTransformation const transformations = [ { id: 'poly1', name: 'Polynomial' }, { id: 'tps', name: 'Thin Plate Spline' }, + { id: 'helmert', name: 'Helmert' }, ]; let currentTargetProjection = 'EPSG:3857'; @@ -703,7 +707,7 @@ }; function getPreview() { - if (gcpList.length < 3) { + if (gcpList.length < 2) { previewMode = 'n/a'; return; } @@ -713,7 +717,7 @@ 'preview', { gcp_geojson: asGeoJSON(), - transformation: currentTransformation, + transformation: gcpList.length == 2 ? 'helmert' : currentTransformation, projection: currentTargetProjection, sesh_id: sessionId, last_preview_id: currentPreviewId, @@ -729,7 +733,7 @@ } function submitSession() { - if (gcpList.length < 3) { + if (gcpList.length < 2) { previewMode = 'n/a'; return; } @@ -869,6 +873,12 @@ > + +

+ You have only two GCPs, but it is highly advisable to have three or more. + Please add more, unless you are unable to do so. +

+ @@ -944,7 +954,13 @@
- + { + if (gcpList.length == 2) { + getModal('modal-submit-helmert').open() + } else { + submitSession() + } + }} title="Save control points" disabled={!enableSave}> { getModal('modal-cancel').open(); From 19a013fe10cb39eef676d45cd8eb6156f47dff14 Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Thu, 4 Jun 2026 15:33:55 -0500 Subject: [PATCH 05/10] complete offset calcuation --- ohmg/georeference/geometry.py | 2 +- ohmg/georeference/georeferencer.py | 149 ++++++++++++++++------------- 2 files changed, 86 insertions(+), 65 deletions(-) diff --git a/ohmg/georeference/geometry.py b/ohmg/georeference/geometry.py index 37ae3d50..cc6471a3 100644 --- a/ohmg/georeference/geometry.py +++ b/ohmg/georeference/geometry.py @@ -7,7 +7,7 @@ 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 X-axis. + pointing from p1 to p2, relative to the positive Y-axis. """ # Calculate differences diff --git a/ohmg/georeference/georeferencer.py b/ohmg/georeference/georeferencer.py index 20e7316a..edd3d7ae 100644 --- a/ohmg/georeference/georeferencer.py +++ b/ohmg/georeference/georeferencer.py @@ -248,50 +248,72 @@ 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) -> Tuple[float, float]: - """Return the image pixel coords from the input GCP""" - return (gcp.GCPPixel, gcp.GCPLine) + 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.""" - print("CALCULATE SCALE") - if len(self.gcps) < 2: + + if len(self.gcps) != 2: raise Exception("Two GCPs are needed to calculate a scale factor") - gcp1_pt = self._geo_coords_from_gcp(self.gcps[0]) - gcp2_pt = self._geo_coords_from_gcp(self.gcps[1]) - gcp1_px = self._pixel_coords_from_gcp(self.gcps[0]) - gcp2_px = self._pixel_coords_from_gcp(self.gcps[1]) + gcp1, gcp2 = self.gcps # distance between the geographic coords in each GCP - pt_dist = math.dist(gcp1_pt, gcp2_pt) + # 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), + ) # distance between the pixel coords in each GCP - px_dist = math.dist(gcp1_px, gcp2_px) + # 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), + ) return pt_dist / px_dist - def _calculate_rotation(self) -> float: + def _calculate_rotation_from_north(self) -> float: """Compares two GCPs and calculates the difference in the angles between the geometric points and the pixel points. - Returns the angle in degrees""" - print("CALCULATE ROTATION") - if len(self.gcps) < 2: + 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") - gcp1_pt = self._geo_coords_from_gcp(self.gcps[0]) - gcp2_pt = self._geo_coords_from_gcp(self.gcps[1]) - gcp1_px = self._pixel_coords_from_gcp(self.gcps[0]) - gcp2_px = self._pixel_coords_from_gcp(self.gcps[1]) + ## 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) - print("gcp1_pt", gcp1_pt) - print("gcp2_pt", gcp2_pt) - print("gcp1_px", gcp1_px) - print("gcp2_px", gcp2_px) + upper_gcp_pt = self._geo_coords_from_gcp(upper_gcp) + lower_gcp_pt = self._geo_coords_from_gcp(lower_gcp) + upper_gcp_px = self._pixel_coords_from_gcp(upper_gcp, cartesian_y=True) + lower_gcp_px = self._pixel_coords_from_gcp(lower_gcp, cartesian_y=True) - pt_degrees = angle_from_coords(gcp1_pt, gcp2_pt) - px_degrees = -angle_from_coords(gcp1_px, gcp2_px) + print("upper_gcp_pt", upper_gcp_pt) + print("lower_gcp_pt", lower_gcp_pt) + print("upper_gcp_px", upper_gcp_px) + print("lower_gcp_px", lower_gcp_px) + + pt_degrees = angle_from_coords(upper_gcp_pt, lower_gcp_pt) + px_degrees = angle_from_coords(upper_gcp_px, lower_gcp_px) print("pt_degrees", pt_degrees) print("px_degrees", px_degrees) @@ -300,52 +322,51 @@ def _calculate_rotation(self) -> float: print("difference", difference) return difference - def _calculate_offsets(self, scale: float, rotation: float) -> Tuple[float, float]: + def _calculate_helmert_offsets(self, scale: float, rotation: float) -> Tuple[float, float]: """Calculates the x and y offsets needed for helmert transformation, based on GCP """ - print("CALCULATE OFFSETS") - gcp1 = self.gcps[0] - print("offset gcp", gcp1.GCPLine, gcp1.GCPPixel) - geo_dist_to_page_edge = gcp1.GCPLine * scale - geo_dist_to_page_top = gcp1.GCPPixel * scale - geo_dist_to_page_corner = math.sqrt(geo_dist_to_page_edge**2 + geo_dist_to_page_top**2) + ## get the gcp that is closest to the top of the page + use_gcp = min(self.gcps, key=lambda x: x.GCPPixel) - print("geo_dist_to_page_edge", geo_dist_to_page_edge) - print("geo_dist_to_page_top", geo_dist_to_page_top) - print("geo_dist_to_page_corner", geo_dist_to_page_corner) + ## use the scale factor to calculate the real distance to page extent + dist_to_page_edge = use_gcp.GCPLine * scale + dist_to_page_top = use_gcp.GCPPixel * scale - print("rotation", rotation) - print("rotation used", rotation + 270) + ## rotation is degrees from positive y-axis, adjust to be degrees from + ## positive x-axis theta1 = 90 - rotation - print("theta1", theta1) + ## further adjustments to normalize if theta1 < 0: theta1 += 180 if theta1 > 180: theta1 -= 180 - print("theta1", theta1) - - m_degrees = 180 - (90 + theta1) - m_radians = math.radians(m_degrees) - n_radians = math.atan(geo_dist_to_page_edge / geo_dist_to_page_top) - theta3 = n_radians + m_radians - print("m_radians", m_radians) - print("n_radians", n_radians) - print("theta3", theta3) - print("m_degrees", math.degrees(m_radians)) - print("n_degrees", math.degrees(n_radians)) - print("theta3_degrees", math.degrees(theta3)) - - a2 = math.cos(theta3) * geo_dist_to_page_corner - o2 = math.sin(theta3) * geo_dist_to_page_corner - - print("a2", a2) - print("o2", o2) - - dx = gcp1.GCPX - a2 - dy = gcp1.GCPY + o2 + ## calculate the target angle q, the angle from the GCP coord + ## to the top left corner of the image, relative to negative + ## x-axis + r_degrees = 180 - (90 + theta1) + r_radians = math.radians(r_degrees) + s_radians = math.atan(dist_to_page_edge / dist_to_page_top) + q_radians = r_radians + s_radians + + ## get the real distance between the GCP coord and the page corner + dist_to_page_corner = math.sqrt(dist_to_page_edge**2 + dist_to_page_top**2) + + ## calculate the offset distances using q angle and distance from + ## GCP to page corner as hypotenuse + x_offset = math.cos(q_radians) * dist_to_page_corner + y_offset = math.sin(q_radians) * dist_to_page_corner + + ## adjust the GCP's coords by adding/subtracting the offsets based + ## on whether the top of the page is above or below the x-axis + if abs(rotation) > 90: + dx = use_gcp.GCPX + x_offset + dy = use_gcp.GCPY - y_offset + else: + dx = use_gcp.GCPX - x_offset + dy = use_gcp.GCPY + y_offset return (dx, dy) @@ -409,12 +430,12 @@ def make_warped_vrt( scale = self._calculate_scale() ## get rotation - rotation_degrees = self._calculate_rotation() - ## adjust to north = 0 and convert to arcseconds - arcseconds = (rotation_degrees + 270) * 3600 + rotation = self._calculate_rotation_from_north() + ## adjust and convert to arcseconds + arcseconds = (rotation + 90) * 3600 ## get the x y offsets - dx, dy = self._calculate_offsets(scale, rotation_degrees) + dx, dy = self._calculate_helmert_offsets(scale, rotation) pipeline = ( "+proj=pipeline " @@ -422,7 +443,7 @@ def make_warped_vrt( f"+step +proj=helmert +x={dx} +y={dy} +theta={arcseconds} +s={scale}" ) - print(pipeline) + logger.debug(f"applying: {pipeline}") wo = gdal.WarpOptions( creationOptions=[ From 54b1e12de2fa15b6333b19a46ecb44e9dda1ddd1 Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Thu, 4 Jun 2026 16:30:38 -0500 Subject: [PATCH 06/10] conditionally allow user to switch to helmert only if 2 GCPs exist --- .../src/components/Georeferencer.svelte | 38 ++++++++++++++----- 1 file changed, 29 insertions(+), 9 deletions(-) diff --git a/ohmg/frontend/svelte_components/src/components/Georeferencer.svelte b/ohmg/frontend/svelte_components/src/components/Georeferencer.svelte index f323fc1b..e9a3ad37 100644 --- a/ohmg/frontend/svelte_components/src/components/Georeferencer.svelte +++ b/ohmg/frontend/svelte_components/src/components/Georeferencer.svelte @@ -121,7 +121,7 @@ leaveOkay = false; enableButtons = true; } - $: enableSave = gcpList.length >= 2 && enableButtons; + $: enableSave = gcpList.length >= minGCPs && enableButtons; let countdown = 10; let timer; @@ -152,12 +152,27 @@ const noteInputElId = 'note-input'; let currentTransformation = 'poly1'; - $: currentTransformation = gcpList.length == 2 ? 'helmert' : currentTransformation - const transformations = [ - { id: 'poly1', name: 'Polynomial' }, - { id: 'tps', name: 'Thin Plate Spline' }, - { id: 'helmert', name: 'Helmert' }, + let minGCPs = 3; + $: transformations = [ + { id: 'poly1', name: 'Polynomial', enabled: true }, + { id: 'tps', name: 'Thin Plate Spline', enabled: true }, + { id: 'helmert', name: 'Helmert 4-param', enabled: gcpList.length == 2 }, ]; + $: { + if (gcpList.length == 3 && currentTransformation == "helmert") { + currentTransformation = "poly1" + } + if (currentTransformation == "helmert") { + minGCPs = 2; + } else { + minGCPs = 3; + } + } + $: { + if (gcpList.length <= 2 && currentTransformation != "helmert") { + previewMode = "n/a" + } + } let currentTargetProjection = 'EPSG:3857'; const availableProjections = [ @@ -707,7 +722,12 @@ }; function getPreview() { - if (gcpList.length < 2) { + if (currentTransformation == "helmert") { + minGCPs = 2; + } else { + minGCPs = 3; + } + if (gcpList.length < minGCPs) { previewMode = 'n/a'; return; } @@ -717,7 +737,7 @@ 'preview', { gcp_geojson: asGeoJSON(), - transformation: gcpList.length == 2 ? 'helmert' : currentTransformation, + transformation: currentTransformation, projection: currentTargetProjection, sesh_id: sessionId, last_preview_id: currentPreviewId, @@ -1026,7 +1046,7 @@ Transformation: From f4a4ccd93cb4227f37a5fdd532020a40d0082a42 Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Thu, 4 Jun 2026 16:34:08 -0500 Subject: [PATCH 07/10] cleanup print statements --- ohmg/georeference/georeferencer.py | 23 +++++++---------------- 1 file changed, 7 insertions(+), 16 deletions(-) diff --git a/ohmg/georeference/georeferencer.py b/ohmg/georeference/georeferencer.py index edd3d7ae..1c3a9712 100644 --- a/ohmg/georeference/georeferencer.py +++ b/ohmg/georeference/georeferencer.py @@ -302,24 +302,15 @@ def _calculate_rotation_from_north(self) -> float: upper_gcp = min(self.gcps, key=lambda x: x.GCPLine) lower_gcp = max(self.gcps, key=lambda x: x.GCPLine) - upper_gcp_pt = self._geo_coords_from_gcp(upper_gcp) - lower_gcp_pt = self._geo_coords_from_gcp(lower_gcp) - upper_gcp_px = self._pixel_coords_from_gcp(upper_gcp, cartesian_y=True) - lower_gcp_px = self._pixel_coords_from_gcp(lower_gcp, cartesian_y=True) - - print("upper_gcp_pt", upper_gcp_pt) - print("lower_gcp_pt", lower_gcp_pt) - print("upper_gcp_px", upper_gcp_px) - print("lower_gcp_px", lower_gcp_px) - - pt_degrees = angle_from_coords(upper_gcp_pt, lower_gcp_pt) - px_degrees = angle_from_coords(upper_gcp_px, lower_gcp_px) - - print("pt_degrees", pt_degrees) - print("px_degrees", px_degrees) + 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), + ) difference = pt_degrees - px_degrees - print("difference", difference) return difference def _calculate_helmert_offsets(self, scale: float, rotation: float) -> Tuple[float, float]: From 39bd4d8b0bbe70d5fd6109c99cc87061bb25f98d Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Fri, 5 Jun 2026 10:17:31 -0500 Subject: [PATCH 08/10] create model for custom permissions objects --- ohmg/conf/http.py | 3 +++ ohmg/core/migrations/0013_permissions.py | 24 ++++++++++++++++++++++++ ohmg/core/models/__init__.py | 1 + ohmg/core/models/permissions.py | 19 +++++++++++++++++++ 4 files changed, 47 insertions(+) create mode 100644 ohmg/core/migrations/0013_permissions.py create mode 100644 ohmg/core/models/permissions.py diff --git a/ohmg/conf/http.py b/ohmg/conf/http.py index 92cb8eb3..0dc3ba6f 100644 --- a/ohmg/conf/http.py +++ b/ohmg/conf/http.py @@ -38,10 +38,13 @@ def user_info_from_request(request): user_info["api_keys"] = user.api_keys user_info["is_authenticated"] = True user_info["is_staff"] = user.is_staff + user_info["perms"] = list(user.get_all_permissions()) + print(json.dumps(user_info, indent=1)) else: user_info = { "is_authenticated": False, "is_staff": False, + "perms": [], } return user_info diff --git a/ohmg/core/migrations/0013_permissions.py b/ohmg/core/migrations/0013_permissions.py new file mode 100644 index 00000000..91531e4e --- /dev/null +++ b/ohmg/core/migrations/0013_permissions.py @@ -0,0 +1,24 @@ +# Generated by Django 4.2.27 on 2026-06-05 09:10 + +from django.db import migrations, models + + +class Migration(migrations.Migration): + + dependencies = [ + ('core', '0012_layerset_xyz_tiles_prefix'), + ] + + operations = [ + migrations.CreateModel( + name='Permissions', + fields=[ + ('id', models.BigAutoField(auto_created=True, primary_key=True, serialize=False, verbose_name='ID')), + ], + options={ + 'permissions': (('use_helmert', 'Can use Helmert transformation'),), + 'managed': False, + 'default_permissions': (), + }, + ), + ] diff --git a/ohmg/core/models/__init__.py b/ohmg/core/models/__init__.py index 8d1d39a4..7d5ba6c7 100644 --- a/ohmg/core/models/__init__.py +++ b/ohmg/core/models/__init__.py @@ -3,4 +3,5 @@ from .layerset import LayerSet, LayerSetCategory # noqa: F401 from .map import Map # noqa: F401 from .mapgroup import MapGroup # noqa: F401 +from .permissions import Permissions # noqa: F401 from .region import Region, RegionCategory # noqa: F401 diff --git a/ohmg/core/models/permissions.py b/ohmg/core/models/permissions.py new file mode 100644 index 00000000..2b68bbb3 --- /dev/null +++ b/ohmg/core/models/permissions.py @@ -0,0 +1,19 @@ +# Source - https://stackoverflow.com/a/37988537 +# Posted by Dmitry, modified by community. See post 'Timeline' for change history +# Retrieved 2026-06-05, License - CC BY-SA 4.0 + +from django.db import models + + +class Permissions(models.Model): + """Permissions model for holding custom Django permissions objects + that are not tied to any specific model.""" + + class Meta: + managed = False # No database table creation or deletion \ + # operations will be performed for this model. + + default_permissions = () # disable "add", "change", "delete" + # and "view" default permissions + + permissions = (("use_helmert", "Can use Helmert transformation"),) From c6d3ac3c50ab9e90cc719b5c03a99cbc636995f9 Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Fri, 5 Jun 2026 10:18:39 -0500 Subject: [PATCH 09/10] restrict helmert availability via permissions --- .../src/components/Georeferencer.svelte | 27 +++++++++++++++---- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/ohmg/frontend/svelte_components/src/components/Georeferencer.svelte b/ohmg/frontend/svelte_components/src/components/Georeferencer.svelte index e9a3ad37..3b6d3340 100644 --- a/ohmg/frontend/svelte_components/src/components/Georeferencer.svelte +++ b/ohmg/frontend/svelte_components/src/components/Georeferencer.svelte @@ -154,9 +154,24 @@ let currentTransformation = 'poly1'; let minGCPs = 3; $: transformations = [ - { id: 'poly1', name: 'Polynomial', enabled: true }, - { id: 'tps', name: 'Thin Plate Spline', enabled: true }, - { id: 'helmert', name: 'Helmert 4-param', enabled: gcpList.length == 2 }, + { + id: 'poly1', + name: 'Polynomial', + enabled: true, + available: true + }, + { + id: 'tps', + name: 'Thin Plate Spline', + enabled: true, + available: true + }, + { + id: 'helmert', + name: 'Helmert 4-param', + enabled: gcpList.length == 2, + available: CONTEXT.user.perms.includes("core.use_helmert") + }, ]; $: { if (gcpList.length == 3 && currentTransformation == "helmert") { @@ -753,7 +768,7 @@ } function submitSession() { - if (gcpList.length < 2) { + if (gcpList.length < minGCPs) { previewMode = 'n/a'; return; } @@ -1046,7 +1061,9 @@ Transformation: From 6a0267aa95d832495a066f9006c07135405511e6 Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Fri, 5 Jun 2026 11:21:18 -0500 Subject: [PATCH 10/10] clean args --- ohmg/georeference/georeferencer.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/ohmg/georeference/georeferencer.py b/ohmg/georeference/georeferencer.py index 1c3a9712..e3ebaa58 100644 --- a/ohmg/georeference/georeferencer.py +++ b/ohmg/georeference/georeferencer.py @@ -443,7 +443,7 @@ def make_warped_vrt( ], coordinateOperation=pipeline, format="VRT", - dstSRS="EPSG:3857", + dstSRS=f"{self.crs_code}", transformerOptions=["SRC_METHOD=NO_GEOTRANSFORM"], dstAlpha=True, resampleAlg="nearest", @@ -472,8 +472,6 @@ def make_warped_vrt( ], format="VRT", dstSRS=f"{self.crs_code}", - # srcNodata=src_nodata, - # srcAlpha=True, dstAlpha=True, resampleAlg="nearest", )