initial plumbing to support non-mercator projections during warping

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) {
           </select>
         </label>
       </div>
+      <!--
+      <div>
+        <label style="margin-top:5px;" title="Set georeferencing transformation">
+          Projection:
+          <select class="trans-select" style="width:151px;" bind:value={currentTargetProjection} on:change={() => { process("preview"); }}>
+            {#each availableProjections as proj}
+              <option value={proj.id}>{proj.name}</option>
+            {/each}
+          </select>
+        </label>
+      </div>-->
     </div>
   </nav>
 </div>

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
 
 

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:

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',

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"])