From 9cc3d2e6c7367eb3eca4a397599ec53a1fa717ee Mon Sep 17 00:00:00 2001 From: David Manthey Date: Thu, 8 Apr 2021 13:39:19 -0400 Subject: [PATCH 1/2] Handle geospatial files with GCP information. --- .../gdal/large_image_source_gdal/__init__.py | 34 +++++++++++++++---- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/sources/gdal/large_image_source_gdal/__init__.py b/sources/gdal/large_image_source_gdal/__init__.py index 62ea0aa08..2b2383063 100644 --- a/sources/gdal/large_image_source_gdal/__init__.py +++ b/sources/gdal/large_image_source_gdal/__init__.py @@ -355,7 +355,7 @@ def getProj4String(self): :returns: The proj4 string or None. """ with self._getDatasetLock: - wkt = self.dataset.GetProjection() + wkt = self.dataset.GetProjection() or self.dataset.GetGCPProjection() if not wkt: if hasattr(self, '_netcdf') or self._getDriver() in {'NITF'}: return InitPrefix + 'epsg:4326' @@ -398,6 +398,20 @@ 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 (not self.dataset.GetProjection() and + 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 +421,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 +535,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 +580,13 @@ 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['GCPProjection'] = self.dataset.GetGCPProjection() + result['GCPs'] = None if not self.dataset.GetGCPs() else [ + {'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 +920,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 From 523624fb9ff982c0eebc4c163a11bfc51285266e Mon Sep 17 00:00:00 2001 From: David Manthey Date: Mon, 12 Apr 2021 09:43:51 -0400 Subject: [PATCH 2/2] Add tests. --- CHANGELOG.md | 1 + .../gdal/large_image_source_gdal/__init__.py | 18 ++++++++----- test/datastore.py | 3 +++ test/test_source_gdal.py | 26 +++++++++++++++++++ 4 files changed, 41 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4ccad0056..ae43e83e8 100644 --- a/CHANGELOG.md +++ b/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) diff --git a/sources/gdal/large_image_source_gdal/__init__.py b/sources/gdal/large_image_source_gdal/__init__.py index 2b2383063..b85b473ad 100644 --- a/sources/gdal/large_image_source_gdal/__init__.py +++ b/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() or self.dataset.GetGCPProjection() + 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' @@ -407,8 +410,7 @@ def _getGeoTransform(self): """ with self._getDatasetLock: gt = self.dataset.GetGeoTransform() - if (not self.dataset.GetProjection() and - self.dataset.GetGCPProjection() and self.dataset.GetGCPs()): + if (self.dataset.GetGCPProjection() and self.dataset.GetGCPs()): gt = gdal.GCPsToGeoTransform(self.dataset.GetGCPs()) return gt @@ -582,11 +584,13 @@ def getInternalMetadata(self, **kwargs): result['RasterYSize'] = self.dataset.RasterYSize result['GeoTransform'] = self._getGeoTransform() result['Projection'] = self.dataset.GetProjection() + result['proj4Projection'] = self.getProj4String() result['GCPProjection'] = self.dataset.GetGCPProjection() - result['GCPs'] = None if not self.dataset.GetGCPs() else [ - {'id': gcp.Id, 'line': gcp.GCPLine, 'pixel': gcp.GCPPixel, - 'x': gcp.GCPX, 'y': gcp.GCPY, 'z': gcp.GCPZ} - for gcp in self.dataset.GetGCPs()] + 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) diff --git a/test/datastore.py b/test/datastore.py index d82bf9f81..536903727 100644 --- a/test/datastore.py +++ b/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 diff --git a/test/test_source_gdal.py b/test/test_source_gdal.py index d91030a2e..dade5c244 100644 --- a/test/test_source_gdal.py +++ b/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']