diff --git a/georeference/components/src/Georeference.svelte b/georeference/components/src/Georeference.svelte index 39edc529..a3e97731 100644 --- a/georeference/components/src/Georeference.svelte +++ b/georeference/components/src/Georeference.svelte @@ -136,6 +136,12 @@ const transformations = [ {id: 'tps', name: 'Thin Plate Spline'}, ]; +let currentTargetProjection = "EPSG:3857" +const availableProjections = [ + {id: 'EPSG:3857', name: 'Pseudo Mercator'}, + {id: 'ESRI:102009', name: 'Lambert North America'}, +]; + // generate a uuid, code from here: // https://www.cloudhadoop.com/2018/10/guide-to-unique-identifiers-uuid-guid function uuid() { @@ -740,6 +746,7 @@ function process(operation){ const data = JSON.stringify({ "gcp_geojson": asGeoJSON(), "transformation": currentTransformation, + "projection": currentTargetProjection, "operation": operation, "sesh_id": session_id, }); @@ -980,6 +987,17 @@ if (DOCUMENT.layer) { + diff --git a/georeference/georeferencer.py b/georeference/georeferencer.py index 7baf69bf..f2972d25 100644 --- a/georeference/georeferencer.py +++ b/georeference/georeferencer.py @@ -1,6 +1,7 @@ import os import time import logging +import requests from osgeo import gdal, osr, ogr from io import StringIO @@ -72,6 +73,22 @@ def get_path_variant(original_path, variant, outdir=None): return os.path.join(outdir, filename) +def retrieve_srs_wkt(code): + + cache_dir = ".srs_cache" + cache_path = os.path.join(cache_dir, f"{code}-wkt.txt") + if os.path.isfile(cache_path): + with open(cache_path, "r") as o: + wkt = o.read() + else: + url = f"https://epsg.io/{code}.wkt" + response = requests.get(url) + wkt = response.content.decode("utf-8") + # with open(cache_path, "w") as o: + # o.write(wkt) + + return wkt + class Georeferencer(object): TRANSFORMATIONS = { @@ -110,11 +127,13 @@ class Georeferencer(object): def __init__(self, epsg_code=None, transformation=None, + crs_code=None, gdal_gcps=[], ): self.gcps = gdal_gcps self.epsg_code = epsg_code + self.crs_code = crs_code self.transformation = None if transformation: self.set_transformation(transformation) @@ -189,7 +208,11 @@ def set_workspace(self, directory): def get_spatial_reference(self): sr = osr.SpatialReference() - sr.ImportFromEPSG(self.epsg_code) + + code = self.crs_code.split(":")[1] + wkt = retrieve_srs_wkt(code) + + sr.ImportFromWkt(wkt) return sr def add_overviews(self, image_path): @@ -233,7 +256,7 @@ def make_warp_options(self, src_path, output_format): f'MAX_GCP_ORDER={self.transformation["gdal_code"]}', ], format=output_format, - dstSRS=f"EPSG:{self.epsg_code}", + dstSRS=f"{self.crs_code}", srcNodata=src_nodata, dstAlpha=True, ) @@ -307,6 +330,30 @@ def make_vrt(self, src_path, output_directory=None): self.run_warp(dst_path, vrt_with_gcps, warp_options) + if self.crs_code != "EPSG:3857": + new_path = dst_path.replace(".vrt", "_3857.vrt") + src_nodata = None + if src_path.endswith(".jpg"): + src_nodata = "255 255 255" + new_options = gdal.WarpOptions( + creationOptions=[ + "TILED=YES", + "COMPRESS=DEFLATE", + ## the following is apparently in the COG spec but doesn't work?? + # "COPY_SOURCE_OVERVIEWS=YES", + ], + transformerOptions = [ + f'DST_SRS={retrieve_srs_wkt(3857)}', + f'MAX_GCP_ORDER={self.transformation["gdal_code"]}', + ], + format="VRT", + dstSRS=f"EPSG:3857", + srcNodata=src_nodata, + dstAlpha=True, + ) + gdal.Warp(new_path,dst_path,options=new_options) + dst_path = new_path + return dst_path diff --git a/georeference/models/sessions.py b/georeference/models/sessions.py index 4ec6f915..93ad15d0 100644 --- a/georeference/models/sessions.py +++ b/georeference/models/sessions.py @@ -506,9 +506,12 @@ def run(self): self.save() try: + # assume EPSG code for now, as making this completely + # flexible is still in-development. see views.py line 277 + crs_code = f"EPSG:{self.data['epsg']}" g = Georeferencer( transformation=self.data['transformation'], - epsg_code=self.data['epsg'], + crs_code=crs_code, ) g.load_gcps_from_geojson(self.data['gcps']) except Exception as e: diff --git a/georeference/utils.py b/georeference/utils.py index 156ae4b7..eaf425a4 100644 --- a/georeference/utils.py +++ b/georeference/utils.py @@ -145,7 +145,7 @@ def add_layer(self, file_path): ' OFFSITE 255 255 255\n', ' TRANSPARENCY 100\n' ' PROJECTION\n', - ' "init=epsg:3857"\n', + ' AUTO\n', ' END\n', ' END # Layer\n', # '\n', diff --git a/georeference/views.py b/georeference/views.py index 16e64732..3da3c65c 100644 --- a/georeference/views.py +++ b/georeference/views.py @@ -228,6 +228,7 @@ def post(self, request, docid): body = json.loads(request.body) gcp_geojson = body.get("gcp_geojson", {}) transformation = body.get("transformation", "poly1") + projection = body.get("projection", "EPSG:3857") operation = body.get("operation", "preview") sesh_id = body.get("sesh_id", None) @@ -243,7 +244,7 @@ def post(self, request, docid): if operation == "preview": # prepare Georeferencer object - g = Georeferencer(epsg_code=3857) + g = Georeferencer(crs_code=projection) g.load_gcps_from_geojson(gcp_geojson) g.set_transformation(transformation) try: @@ -272,7 +273,11 @@ def post(self, request, docid): if operation == "submit": - sesh.data['epsg'] = 3857 + # ultimately, should be putting the whole "EPSG:3857" + # in the session data, but for now stick with just the number + # see models.sessions.py line 510 + epsg_code = int(projection.split(":")[1]) + sesh.data['epsg'] = epsg_code sesh.data['gcps'] = gcp_geojson sesh.data['transformation'] = transformation sesh.save(update_fields=["data"])