Merge pull request #588 from girder/handle-gcp-projections

Handle geospatial files with GCP information

CHANGELOG.md

@@ -6,6 +6,7 @@
 - Allow converting a single frame of a multiframe image (#579)
 - Add a convert endpoint to the Girder plugin (#578)
 - Added support for creating Aperio svs files (#580)
+- Added support for geospatial files with GCP (#588)
 
 ### Improvements
 - More untiled tiff files are handles by the bioformats reader (#569)

sources/gdal/large_image_source_gdal/__init__.py

@@ -355,7 +355,10 @@ def getProj4String(self):
         :returns: The proj4 string or None.
         """
         with self._getDatasetLock:
-            wkt = self.dataset.GetProjection()
+            if self.dataset.GetGCPs() and self.dataset.GetGCPProjection():
+                wkt = self.dataset.GetGCPProjection()
+            else:
+                wkt = self.dataset.GetProjection()
         if not wkt:
             if hasattr(self, '_netcdf') or self._getDriver() in {'NITF'}:
                 return InitPrefix + 'epsg:4326'
@@ -398,6 +401,19 @@ def getNativeMagnification(self):
             'mm_y': scale * 100 if scale else None,
         }
 
+    def _getGeoTransform(self):
+        """
+        Get the GeoTransform.  If GCPs are used, get the appropriate transform
+        for those.
+
+        :returns: a six-component array with the transform
+        """
+        with self._getDatasetLock:
+            gt = self.dataset.GetGeoTransform()
+            if (self.dataset.GetGCPProjection() and self.dataset.GetGCPs()):
+                gt = gdal.GCPsToGeoTransform(self.dataset.GetGCPs())
+        return gt
+
     def getBounds(self, srs=None):
         """
         Returns bounds of the image.
@@ -407,8 +423,7 @@ def getBounds(self, srs=None):
             used.  None if we don't know the original projection.
         """
         if srs not in self._bounds:
-            with self._getDatasetLock:
-                gt = self.dataset.GetGeoTransform()
+            gt = self._getGeoTransform()
             nativeSrs = self.getProj4String()
             if not nativeSrs:
                 self._bounds[srs] = None
@@ -522,7 +537,10 @@ def getOneBandInformation(self, band):
     def getMetadata(self):
         with self._getDatasetLock:
             metadata = {
-                'geospatial': bool(self.dataset.GetProjection() or hasattr(self, '_netcdf')),
+                'geospatial': bool(
+                    self.dataset.GetProjection() or
+                    (self.dataset.GetGCPProjection() and self.dataset.GetGCPs()) or
+                    hasattr(self, '_netcdf')),
                 'levels': self.levels,
                 'sizeX': self.sizeX,
                 'sizeY': self.sizeY,
@@ -564,8 +582,15 @@ def getInternalMetadata(self, **kwargs):
             result['fileList'] = self.dataset.GetFileList()
             result['RasterXSize'] = self.dataset.RasterXSize
             result['RasterYSize'] = self.dataset.RasterYSize
-            result['GeoTransform'] = self.dataset.GetGeoTransform()
+            result['GeoTransform'] = self._getGeoTransform()
+            result['Projection'] = self.dataset.GetProjection()
+            result['proj4Projection'] = self.getProj4String()
             result['GCPProjection'] = self.dataset.GetGCPProjection()
+            if self.dataset.GetGCPs():
+                result['GCPs'] = [{
+                    'id': gcp.Id, 'line': gcp.GCPLine, 'pixel': gcp.GCPPixel,
+                    'x': gcp.GCPX, 'y': gcp.GCPY, 'z': gcp.GCPZ}
+                    for gcp in self.dataset.GetGCPs()]
             result['Metadata'] = self.dataset.GetMetadata_List()
             for key in ['IMAGE_STRUCTURE', 'SUBDATASETS', 'GEOLOCATION', 'RPC']:
                 metadatalist = self.dataset.GetMetadata_List(key)
@@ -899,8 +924,7 @@ def toNativePixelCoordinates(self, x, y, proj=None, roundResults=True):
         outProj = self._proj4Proj(self.getProj4String())
         px, py = pyproj.Transformer.from_proj(inProj, outProj, always_xy=True).transform(x, y)
         # convert to native pixel coordinates
-        with self._getDatasetLock:
-            gt = self.dataset.GetGeoTransform()
+        gt = self._getGeoTransform()
         d = gt[2] * gt[4] - gt[1] * gt[5]
         x = (gt[0] * gt[5] - gt[2] * gt[3] - gt[5] * px + gt[2] * py) / d
         y = (gt[1] * gt[3] - gt[0] * gt[4] + gt[4] * px - gt[1] * py) / d

test/datastore.py

@@ -36,6 +36,9 @@
     # Thematic landcover sample
     # Source: landcover_sample_1000.tif
     'landcover_sample_1000.tif': 'sha512:0a1e8c4cf29174b19ddece9a2deb7230a31e819ee78b5ec264feda6356abf60d63763106f1ddad9dd04106d383fd0867bf2db55be0552c30f38ffb530bf72dec',  # noqa
+    # Geospatial image defined by GCP
+    # Source: extracted from a public landsat image
+    'region_gcp.tiff': 'sha512:b49d753c0a04888882da60917a6578707ad2d6f5a8f1de3b21b2ae498ad6140ee8a1779a4690cc24fa7b1cbae46e86c3529d2fbd88f392a344fbe549b4447f23',  # noqa
     # Multiframe OME tif
     # Source: sample.ome.tif
     'sample.ome.tif': 'sha512:a65ad3d67bbc8f56eb3aa7c7d6a43326e50fdef7c0a0c7ace2b477c2dfda2e810c1acc148be2ee2da9a9aa7b0195032938f36da4117a7c3f46302d4fbb1e8173',  # noqa

test/test_source_gdal.py

@@ -348,3 +348,29 @@ def testGetRegionWithProjection():
     region, _ = ts.getRegion(output=dict(maxWidth=1024, maxHeight=1024),
                              format=constants.TILE_FORMAT_NUMPY)
     assert region.shape == (1024, 1024, 4)
+
+
+def testGCPProjection():
+    imagePath = datastore.fetch('region_gcp.tiff')
+    source = large_image_source_gdal.open(imagePath)
+    tileMetadata = source.getMetadata()
+    assert tileMetadata['tileWidth'] == 256
+    assert tileMetadata['tileHeight'] == 256
+    assert tileMetadata['sizeX'] == 1204
+    assert tileMetadata['sizeY'] == 512
+    assert tileMetadata['levels'] == 4
+    assert tileMetadata['geospatial']
+
+    source = large_image_source_gdal.open(imagePath, projection='EPSG:3857')
+    tileMetadata = source.getMetadata()
+    assert tileMetadata['tileWidth'] == 256
+    assert tileMetadata['tileHeight'] == 256
+    assert tileMetadata['sizeX'] == 524288
+    assert tileMetadata['sizeY'] == 524288
+    assert tileMetadata['levels'] == 12
+    assert tileMetadata['bounds']['xmax'] == pytest.approx(-10753925, 1)
+    assert tileMetadata['bounds']['xmin'] == pytest.approx(-10871650, 1)
+    assert tileMetadata['bounds']['ymax'] == pytest.approx(3949393, 1)
+    assert tileMetadata['bounds']['ymin'] == pytest.approx(3899358, 1)
+    assert tileMetadata['bounds']['srs'] == '+init=epsg:3857'
+    assert tileMetadata['geospatial']