From 66d74e8da02a2f8dcb724042de84c541ed76e54f Mon Sep 17 00:00:00 2001 From: Adam Cox Date: Tue, 14 Jul 2026 10:15:12 -0500 Subject: [PATCH 1/2] better handling auto switch to poly1; address #393 --- .../interfaces/Georeferencer.svelte | 74 ++++++++++--------- 1 file changed, 38 insertions(+), 36 deletions(-) diff --git a/ohmg/frontend/svelte_components/src/components/interfaces/Georeferencer.svelte b/ohmg/frontend/svelte_components/src/components/interfaces/Georeferencer.svelte index 73a53141..1785c5fe 100644 --- a/ohmg/frontend/svelte_components/src/components/interfaces/Georeferencer.svelte +++ b/ohmg/frontend/svelte_components/src/components/interfaces/Georeferencer.svelte @@ -173,21 +173,6 @@ available: CONTEXT.user.perms.includes("core.use_helmert") }, ]; - $: { - 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 = [ @@ -590,6 +575,19 @@ note: props.note, }); }); + + // special handling for helmert, do this before the preview is regenerated + 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" + } getPreview(); } @@ -736,6 +734,16 @@ return featureCollection; }; + const preparePayload = function () { + return { + gcp_geojson: asGeoJSON(), + transformation: currentTransformation, + projection: currentTargetProjection, + sesh_id: sessionId, + last_preview_id: currentPreviewId, + } + } + function getPreview() { if (currentTransformation == "helmert") { minGCPs = 2; @@ -750,13 +758,7 @@ `/georeference/${REGION.id}/`, CONTEXT.ohmg_post_headers, 'preview', - { - gcp_geojson: asGeoJSON(), - transformation: currentTransformation, - projection: currentTargetProjection, - sesh_id: sessionId, - last_preview_id: currentPreviewId, - }, + preparePayload(), (result) => { previewMode = previewMode == 'n/a' ? 'transparent' : previewMode; // updating this variable will trigger the preview layer to be @@ -779,13 +781,7 @@ `/georeference/${REGION.id}/`, CONTEXT.ohmg_post_headers, 'submit', - { - gcp_geojson: asGeoJSON(), - transformation: currentTransformation, - projection: currentTargetProjection, - sesh_id: sessionId, - last_preview_id: currentPreviewId, - }, + preparePayload(), () => { window.location.href = `/map/${REGION.map}`; }, @@ -901,17 +897,23 @@ }} on:unload={cancelSession} /> -
- Create 3 or more ground control points to georeference this document. + Create three or more ground control points to georeference this document. First, click on a recognizable spot + in the document on the left to begin the GCP. Then, click on that same location in the web map on the right + to complete it. Learn more about the process Learn more + external={true}>in the docs.
- +

- 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. + You have only created two GCPs, but we recommend having three or more. + If possible, please add more before submitting your work.

Date: Tue, 14 Jul 2026 11:44:58 -0500 Subject: [PATCH 2/2] refactor transformer for iiif #391; cleanup --- ohmg/extensions/iiif.py | 37 +++++----- ohmg/extensions/urls.py | 2 +- ohmg/georeference/georeferencer.py | 111 ++++++++++------------------- 3 files changed, 59 insertions(+), 91 deletions(-) diff --git a/ohmg/extensions/iiif.py b/ohmg/extensions/iiif.py index f1d0c492..6284a1e5 100644 --- a/ohmg/extensions/iiif.py +++ b/ohmg/extensions/iiif.py @@ -34,20 +34,23 @@ def get_target(self): wgs84.ImportFromEPSG(4326) mask_geojson = json.loads(self.layer.mask.geojson) coords = mask_geojson["coordinates"][0] - ct = CoordTransform(SpatialReference("WGS84"), SpatialReference("EPSG:3857")) + + target_crs = f"EPSG:{self.region.gcpgroup.crs_epsg}" + ct = CoordTransform(SpatialReference("WGS84"), SpatialReference(target_crs)) polygon = Polygon(coords) polygon.transform(ct) - g = Georeferencer( - crs=f"EPSG:{self.region.gcpgroup.crs_epsg}", - transformation=self.region.gcpgroup.transformation, - gcps_geojson=self.region.gcpgroup.as_geojson, - ) in_path = ( self.region.file.url if self.region.file.url.startswith("http") else self.region.file.path ) + + g = Georeferencer( + crs=target_crs, + transformation=self.region.gcpgroup.transformation, + gcps_geojson=self.region.gcpgroup.as_geojson, + ) g.make_gcps_vrt(in_path) ds = gdal.Open(g.gcps_vrt.get_vsi_url()) transformer = gdal.Transformer( @@ -57,10 +60,7 @@ def get_target(self): ds, # Transformer options that are ultimately passed to 'GDALCreateGenImgProjTransformer2()' # https://gdal.org/api/gdal_alg.html#_CPPv432GDALCreateGenImgProjTransformer212GDALDatasetH12GDALDatasetHPPc - [ - ## need to update this, different number if thin plate spline - "MAX_GCP_ORDER=1", - ], + g.make_transformer_options(), ) transposed, status = transformer.TransformPoints(False, polygon.coords[0]) coords_str = [f"{i[0]},{i[1]}" for i in transposed] @@ -114,14 +114,15 @@ def get_gcps(self): return features def get_body(self): - if self.region.gcpgroup.transformation == "poly1": - transformation = {"type": "polynomial", "options": {"order": 1}} - elif self.region.gcpgroup.transformation == "tps": - transformation = { - "type": "thinPlateSpline", - } - else: - raise Exception("invalid transformation", self.region.gcpgroup.transformation) + match t := self.region.gcpgroup.transformation: + case "poly1": + transformation = {"type": "polynomial", "options": {"order": 1}} + case "tps": + transformation = {"type": "thinPlateSpline"} + case "helmert": + transformation = {"type": "helmert"} + case _: + raise Exception(f"invalid transformation: {t}") body = { "id": full_reverse("iiif_gcps_view", args=(self.layerid,)), diff --git a/ohmg/extensions/urls.py b/ohmg/extensions/urls.py index 65408ed4..e85ebc56 100644 --- a/ohmg/extensions/urls.py +++ b/ohmg/extensions/urls.py @@ -1,5 +1,6 @@ from django.urls import path, register_converter +from ohmg.extensions.feeds import PlaceFeed from ohmg.places.converters import PlaceConverter from .views import ( @@ -9,7 +10,6 @@ IIIFResourceView, IIIFSelectorView, ) -from ohmg.extensions.feeds import PlaceFeed register_converter(PlaceConverter, "place-slug") diff --git a/ohmg/georeference/georeferencer.py b/ohmg/georeference/georeferencer.py index d621bb98..36ac9f37 100644 --- a/ohmg/georeference/georeferencer.py +++ b/ohmg/georeference/georeferencer.py @@ -30,6 +30,7 @@ class HelmertParams: offset_x: float offset_y: float + gdal.SetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS") TRANSFORMATION_LOOKUP = { @@ -300,6 +301,34 @@ def _get_helmert_params(self) -> HelmertParams: offset_y=offset_y, ) + def _get_helmert_proj_pipeline(self) -> str: + ## fit the four-parameter Helmert model in one shot + params = self._get_helmert_params() + + ## convert the rotation to arcseconds for the proj pipeline + arcseconds = (params.rotation + 90) * 3600 + + pipeline = ( + "+proj=pipeline " + "+step +proj=axisswap +order=2,1 " + f"+step +proj=helmert +x={params.offset_x} +y={params.offset_y} " + f"+theta={arcseconds} +s={params.scale}" + ) + + return pipeline + + def make_transformer_options(self) -> list[str]: + match self.transformation["id"]: + case "helmert": + pipeline = self._get_helmert_proj_pipeline() + logger.debug(f"applying: {pipeline}") + return ["SRC_METHOD=NO_GEOTRANSFORM", f"COORDINATE_OPERATION={pipeline}"] + case _: + return [ + f"DST_SRS={self.crs_wkt}", + f"MAX_GCP_ORDER={self.transformation['gdal_code']}", + ] + def cleanup_files(self): for vrt in [ self.gcps_vrt, @@ -355,61 +384,16 @@ def make_warped_vrt( self.make_gcps_vrt(src_path, out_name) self.warped_vrt = VRTHandler(out_name, as_variant="modified") - if self.transformation["id"] == "helmert": - ## fit the four-parameter Helmert model in one shot - params = self._get_helmert_params() - - ## convert the rotation to arcseconds for the proj pipeline - arcseconds = (params.rotation + 90) * 3600 - - pipeline = ( - "+proj=pipeline " - "+step +proj=axisswap +order=2,1 " - f"+step +proj=helmert +x={params.offset_x} +y={params.offset_y} " - f"+theta={arcseconds} +s={params.scale}" - ) - - logger.debug(f"applying: {pipeline}") - - wo = gdal.WarpOptions( - creationOptions=[ - "BLOCKXSIZE=512", - "BLOCKYSIZE=512", - ], - coordinateOperation=pipeline, - format="VRT", - dstSRS=f"{self.crs_code}", - 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}", - dstAlpha=True, - resampleAlg="nearest", - ) + wo = gdal.WarpOptions( + creationOptions=[ + "BLOCKXSIZE=512", + "BLOCKYSIZE=512", + ], + transformerOptions=self.make_transformer_options(), + format="VRT", + dstAlpha=True, + resampleAlg="nearest", + ) try: gdal.Warp(str(self.warped_vrt.get_path()), self.gcps_vrt.get_vsi_url(), options=wo) @@ -432,15 +416,6 @@ def make_trimmed_vrt(self, src_path, multimask_json_file: Path, layer_name: str) cutlineLayer=multimask_json_file.stem, cutlineWhere=f"layer='{layer_name}'", cropToCutline=True, - # srcAlpha = True, - # dstAlpha = True, - # creationOptions= [ - # 'COMPRESS=JPEG', - # ] - # creationOptions= [ - # 'COMPRESS=DEFLATE', - # 'PREDICTOR=2', - # ] ) gdal.Warp(str(self.trimmed_vrt.get_path()), self.warped_vrt.get_vsi_url(), options=wo) @@ -458,17 +433,9 @@ def make_cog( self.cog = Path(settings.TEMP_DIR, Path(src_path).stem + "-modified.tif") to = gdal.TranslateOptions( - # format="GTiff", format="COG", - # maskBand="mask", creationOptions=[ - # 'COMPRESS=DEFLATE', - # 'PREDICTOR=2', "COMPRESS=JPEG", - # 'TILED=YES', - # 'BLOCKXSIZE=512', - # 'BLOCKYSIZE=512', - # "PHOTOMETRIC=YCBCR", "TILING_SCHEME=GoogleMapsCompatible", ], resampleAlg="nearest",