Support for non-UTM coordinate systems

app/admin.py

@@ -44,7 +44,28 @@ def has_add_permission(self, request):
     list_display = ('id', 'name', 'project', 'processing_node', 'created_at', 'status', 'last_error')
     list_filter = ('status', 'project',)
     search_fields = ('id', 'name', 'project__name')
-
+    exclude = ('orthophoto_extent', 'dsm_extent', 'dtm_extent', 'crop', )
+    readonly_fields = ('orthophoto_extent_wkt', 'dsm_extent_wkt', 'dtm_extent_wkt', 'crop_wkt', )
+
+    def orthophoto_extent_wkt(self, obj):
+        if obj.orthophoto_extent:
+            return obj.orthophoto_extent.wkt
+        return None
+
+    def dsm_extent_wkt(self, obj):
+        if obj.dsm_extent:
+            return obj.dsm_extent.wkt
+        return None
+
+    def dtm_extent_wkt(self, obj):
+        if obj.dtm_extent:
+            return obj.dtm_extent.wkt
+        return None
+
+    def crop_wkt(self, obj):
+        if obj.crop:
+            return obj.crop.wkt
+        return None
 
 admin.site.register(Task, TaskAdmin)
 

app/api/tasks.py

@@ -36,7 +36,7 @@
 from app.security import path_traversal_check
 from django.utils.translation import gettext_lazy as _
 from .fields import PolygonGeometryField
-from app.geoutils import geom_transform_wkt_bbox
+from app.geoutils import geom_transform_wkt_bbox, get_srs_name_units_from_epsg
 from webodm import settings
 
 def flatten_files(request_files):
@@ -59,6 +59,7 @@ class TaskSerializer(serializers.ModelSerializer):
     extent = serializers.SerializerMethodField()
     tags = TagsField(required=False)
     crop = PolygonGeometryField(required=False, allow_null=True)
+    srs = serializers.SerializerMethodField()
 
     def get_processing_node_name(self, obj):
         if obj.processing_node is not None:
@@ -90,6 +91,9 @@ def get_can_rerun_from(self, obj):
 
     def get_extent(self, obj):
         return obj.get_extent()
+
+    def get_srs(self, obj):
+        return get_srs_name_units_from_epsg(obj.epsg)
 
     class Meta:
         model = models.Task

app/api/tiler.py

@@ -25,12 +25,13 @@
 from .hillshade import LightSource
 from .formulas import lookup_formula, get_algorithm_list, get_auto_bands
 from .tasks import TaskNestedView
-from app.geoutils import geom_transform_wkt_bbox
+from app.geoutils import geom_transform_wkt_bbox, get_rasterio_to_meters_factor
 from rest_framework import exceptions
 from rest_framework.response import Response
 from worker.tasks import export_raster, export_pointcloud
 from django.utils.translation import gettext as _
 import warnings
+from functools import lru_cache
 
 # Disable: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix be returned.
 warnings.filterwarnings("ignore", category=NotGeoreferencedWarning)
@@ -47,7 +48,19 @@
 for custom_colormap in custom_colormaps:
     colormap = colormap.register(custom_colormap)
 
-
+@lru_cache(maxsize=128)
+def get_colormap_encoded_values(cmap):
+    values = colormap.get(cmap).values()
+    # values = [[R, G, B, A], [R, G, B, A], ...]
+    encoded_values = []
+    for rgba in values:
+        # Pack R, G, B, A (each 0-255) into a 32-bit integer
+        # Format: 0xRRGGBBAA (R in most significant byte)
+        encoded = (rgba[0] << 24) | (rgba[1] << 16) | (rgba[2] << 8) | rgba[3]
+        encoded_values.append(encoded)
+
+    return encoded_values
+
 def get_zoom_safe(src_dst):
     minzoom, maxzoom = src_dst.spatial_info["minzoom"], src_dst.spatial_info["maxzoom"]
     if maxzoom < minzoom:
@@ -171,8 +184,13 @@ def get(self, request, pk=None, project_pk=None, tile_type=""):
         raster_path = get_raster_path(task, tile_type)
         if not os.path.isfile(raster_path):
             raise exceptions.NotFound()
+
+        to_meter = 1.0
         try:
             with COGReader(raster_path) as src:
+                if tile_type in ['dsm', 'dtm']:
+                    to_meter = get_rasterio_to_meters_factor(src.dataset)
+
                 band_count = src.dataset.meta['count']
                 if boundaries_feature is not None:
                     cutline = create_cutline(src.dataset, boundaries_feature, CRS.from_string('EPSG:4326'))
@@ -265,13 +283,22 @@ def get(self, request, pk=None, project_pk=None, tile_type=""):
         info['color_maps'] = []
         info['algorithms'] = algorithms
         info['auto_bands'] = auto_bands
+
+        if to_meter != 1.0:
+            for b in info['statistics']:
+                info['statistics'][b]['min'] *= to_meter
+                info['statistics'][b]['max'] *= to_meter
+                info['statistics'][b]['std'] *= to_meter
+                info['statistics'][b]['percentiles'][0] *= to_meter
+                info['statistics'][b]['percentiles'][1] *= to_meter
+                info['statistics'][b]['histogram'][1] = [n * to_meter for n in info['statistics'][b]['histogram'][1]]
 
         if colormaps:
             for cmap in colormaps:
                 try:
                     info['color_maps'].append({
                         'key': cmap,
-                        'color_map': colormap.get(cmap).values(),
+                        'color_map': get_colormap_encoded_values(cmap),
                         'label': cmap_labels.get(cmap, cmap)
                     })
                 except FileNotFoundError:
@@ -369,6 +396,7 @@ def get(self, request, pk=None, project_pk=None, tile_type="", z="", x="", y="",
         if not os.path.isfile(url):
             raise exceptions.NotFound()
 
+        to_meter = 1.0
         with COGReader(url) as src:
             if not src.tile_exists(z, x, y):
                 raise exceptions.NotFound(_("Outside of bounds"))
@@ -393,6 +421,9 @@ def get(self, request, pk=None, project_pk=None, tile_type="", z="", x="", y="",
             else:
                 vrt_options = None
 
+            if tile_type in ['dsm', 'dtm']:
+                to_meter = get_rasterio_to_meters_factor(src.dataset)
+
             # Handle N-bands datasets for orthophotos (not plant health)
             if tile_type == 'orthophoto' and expr is None:
                 ci = src.dataset.colorinterp
@@ -451,6 +482,8 @@ def get(self, request, pk=None, project_pk=None, tile_type="", z="", x="", y="",
             intensity = None
             try:
                 rescale_arr = list(map(float, rescale.split(",")))
+                if tile_type in ['dsm', 'dtm']:
+                    rescale_arr = [v / to_meter for v in rescale_arr]
             except ValueError:
                 raise exceptions.ValidationError(_("Invalid rescale value"))
 

app/geoutils.py

@@ -1,6 +1,17 @@
 import rasterio.warp
 import numpy as np
 from rasterio.crs import CRS
+from rasterio.warp import transform_bounds
+from osgeo import osr, gdal
+from functools import lru_cache
+osr.DontUseExceptions()
+
+SUPPORTED_UNITS = ["m", "ft", "US survey foot"]
+UNIT_TO_M = {
+    "m": 1.0,
+    "ft": 0.3048,
+    "US survey foot": 1200.0 / 3937.0,
+}
 
 # GEOS has some weird bug where
 # we can't simply call geom.tranform(srid)
@@ -79,3 +90,94 @@ def geom_transform(geom, epsg):
         return list(zip(tx, ty))
     else:
         raise ValueError("Cannot transform complex geometries to WKT")
+
+
+def epsg_from_wkt(wkt):
+    srs = osr.SpatialReference()
+    if srs.ImportFromWkt(wkt) != 0:
+        return None
+
+    epsg = srs.GetAuthorityCode(None)
+    if epsg is not None:
+        return None
+
+    # Try to get the 2D component
+    if srs.IsCompound():
+        if srs.DemoteTo2D() != 0:
+            return None
+
+    epsg = srs.GetAuthorityCode(None)
+    if epsg is not None:
+        return epsg
+
+
+def get_raster_bounds_wkt(raster_path, target_srs="EPSG:4326"):
+    with rasterio.open(raster_path) as src:
+        if src.crs is None:
+            return None
+
+        left, bottom, right, top = src.bounds
+        w, s, e, n = transform_bounds(
+            src.crs, 
+            target_srs, 
+            left, bottom, right, top
+        )
+
+        wkt = f"POLYGON(({w} {s}, {w} {n}, {e} {n}, {e} {s}, {w} {s}))"
+        return wkt
+
+@lru_cache(maxsize=1000)
+def get_srs_name_units_from_epsg(epsg):
+    if epsg is None:
+        return {'name': '', 'units': 'm'}
+
+    srs = osr.SpatialReference()
+    if srs.ImportFromEPSG(epsg) != 0:
+        return {'name': '', 'units': 'm'}
+
+    name = srs.GetAttrValue("PROJCS")
+    if name is None:
+        name = srs.GetAttrValue("GEOGCS")
+
+    if name is None:
+        return {'name': '', 'units': 'm'}
+
+    units = srs.GetAttrValue('UNIT')
+    if units is None:
+        units = 'm'  # Default to meters
+    elif units not in SUPPORTED_UNITS:
+        units = 'm' # Unsupported
+
+    return {'name': name, 'units': units}
+
+def get_rasterio_to_meters_factor(rasterio_ds):
+    units = rasterio_ds.units
+    if len(units) >= 1:
+        unit = units[0]
+        if unit is not None and unit != "" and unit in SUPPORTED_UNITS:
+            return UNIT_TO_M.get(unit, 1.0)
+    return 1.0
+
+
+def get_raster_dem_to_meters_factor(raster_path):
+    unit = get_raster_dem_units(raster_path)
+    return UNIT_TO_M.get(unit, 1.0)
+
+def get_raster_dem_units(raster_path):
+    try:
+        ds = gdal.Open(raster_path, gdal.GA_ReadOnly)
+        if ds is None:
+            raise IOError(f"Cannot open {raster_path}")
+
+        band = ds.GetRasterBand(1)
+        unit = band.GetUnitType()
+        ds = None
+
+        if unit is None or unit == "":
+            return "m"
+        elif unit in SUPPORTED_UNITS:
+            return unit
+        else:
+            return "m"
+    except Exception as e:
+        return "m"

app/models/task.py

@@ -22,8 +22,6 @@
 import requests
 from PIL import Image
 Image.MAX_IMAGE_PIXELS = 4096000000
-from django.contrib.gis.gdal import GDALRaster
-from django.contrib.gis.gdal import OGRGeometry
 from django.contrib.gis.geos import GEOSGeometry
 from django.contrib.postgres import fields
 from django.core.files.uploadedfile import InMemoryUploadedFile
@@ -41,7 +39,7 @@
 from app.pointcloud_utils import is_pointcloud_georeferenced
 from app.testwatch import testWatch
 from app.security import path_traversal_check
-from app.geoutils import geom_transform
+from app.geoutils import geom_transform, epsg_from_wkt, get_raster_bounds_wkt, get_srs_name_units_from_epsg
 from nodeodm import status_codes
 from nodeodm.models import ProcessingNode
 from pyodm.exceptions import NodeResponseError, NodeConnectionError, NodeServerError, OdmError
@@ -1012,32 +1010,16 @@ def extract_assets_and_complete(self):
                 except IOError as e:
                     logger.warning("Cannot create Cloud Optimized GeoTIFF for %s (%s). This will result in degraded visualization performance." % (raster_path, str(e)))
 
-                # Read extent and SRID
-                raster = GDALRaster(raster_path)
-                extent = OGRGeometry.from_bbox(raster.extent)
-
-                # Make sure PostGIS supports it
-                with connection.cursor() as cursor:
-                    cursor.execute("SELECT SRID FROM spatial_ref_sys WHERE SRID = %s", [raster.srid])
-                    if cursor.rowcount == 0:
-                        raise NodeServerError(gettext("Unsupported SRS %(code)s. Please make sure you picked a supported SRS.") % {'code': str(raster.srid)})
-
-                # It will be implicitly transformed into the SRID of the model’s field
-                # self.field = GEOSGeometry(...)
-                setattr(self, field, GEOSGeometry(extent.wkt, srid=raster.srid))
-
-                logger.info("Populated extent field with {} for {}".format(raster_path, self))
+                # Read extent
+                extent_wkt = get_raster_bounds_wkt(raster_path)
+                if extent_wkt is not None:
+                    extent = GEOSGeometry(extent_wkt, srid=4326)
+                    setattr(self, field, extent)
+                    logger.info("Populated extent field with {} for {}".format(raster_path, self))
+                else:
+                    logger.warning("Cannot populate extent field with {} for {}, not georeferenced".format(raster_path, self))
 
         self.check_ept()
-
-        # Flushes the changes to the *_extent fields
-        # and immediately reads them back into Python
-        # This is required because GEOS screws up the X/Y conversion
-        # from the raster CRS to 4326, whereas PostGIS seems to do it correctly :/
-        self.status = status_codes.RUNNING # avoid telling clients that task is completed prematurely
-        self.save()
-        self.refresh_from_db()
-
         self.update_available_assets_field()
         self.update_epsg_field()
         self.update_orthophoto_bands_field()
@@ -1158,6 +1140,7 @@ def get_map_items(self):
                     'camera_shots': camera_shots,
                     'ground_control_points': ground_control_points,
                     'epsg': self.epsg,
+                    'srs': get_srs_name_units_from_epsg(self.epsg),
                     'orthophoto_bands': self.orthophoto_bands,
                     'crop': self.crop is not None,
                     'extent': self.get_extent(),
@@ -1225,8 +1208,18 @@ def update_epsg_field(self, commit=False):
                 try:
                     with rasterio.open(asset_path) as f:
                         if f.crs is not None:
-                            epsg = f.crs.to_epsg()
-                            break # We assume all assets are in the same CRS
+                            code = f.crs.to_epsg()
+                            if code is not None:
+                                epsg = code
+                                break # We assume all assets are in the same CRS
+                            else:
+                                # Try to get code from WKT
+                                wkt = f.crs.to_wkt()
+                                if wkt is not None:
+                                    code = epsg_from_wkt(wkt)
+                                    if code is not None:
+                                        epsg = code
+                                        break
                 except Exception as e:
                     logger.warning(e)