Skip to content

Helmert iiif #397

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mradamcox wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Helmert iiif #397

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
66d74e8
better handling auto switch to poly1; address #393
mradamcox Jul 14, 2026
f770250
refactor transformer for iiif #391; cleanup
mradamcox Jul 14, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 19 additions & 18 deletions ohmg/extensions/iiif.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand All @@ -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]
Expand Down Expand Up @@ -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,)),
Expand Down
2 changes: 1 addition & 1 deletion ohmg/extensions/urls.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 (
Expand All @@ -9,7 +10,6 @@
IIIFResourceView,
IIIFSelectorView,
)
from ohmg.extensions.feeds import PlaceFeed

register_converter(PlaceConverter, "place-slug")

Expand Down
74 changes: 38 additions & 36 deletions ohmg/frontend/svelte_components/src/components/interfaces/Georeferencer.svelte
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 = [
Expand Down Expand Up @@ -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();
}

Expand Down Expand Up @@ -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;
Expand All @@ -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
Expand All @@ -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}`;
},
Expand Down Expand Up @@ -901,17 +897,23 @@
}}
on:unload={cancelSession}
/>
<div style="height:25px;">
Create 3 or more ground control points to georeference this document. <Link
<div>
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 <Link
href="https://docs.oldinsurancemaps.net/guides/georeferencing/"
external={true}>Learn more</Link
>
external={true}>in the docs</Link
>.
</div>

<ModalConfirm id="modal-submit-helmert" yesAction={submitSession}>
<ModalConfirm id="modal-submit-helmert"
yesButtonText="Submit with two GCPs"
yesAction={submitSession}
noButtonText="Back to georeferencing"
>
<p>
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.
</p>
</ModalConfirm>
<ModalConfirm id="modal-cancel"
Expand Down
111 changes: 39 additions & 72 deletions ohmg/georeference/georeferencer.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ class HelmertParams:
offset_x: float
offset_y: float


gdal.SetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS")

TRANSFORMATION_LOOKUP = {
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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)
Expand All @@ -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)

Expand All @@ -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",
Expand Down