Use gdal warp when not applying a style on geospatial sources.

CHANGELOG.md

@@ -2,6 +2,9 @@
 
 ## Unreleased
 
+### Features
+- Regions can be output as tiled tiffs with scale or geospatial metadata (#594)
+
 ### Improvements
 - Cache histogram requests (#598)
 
@@ -12,7 +15,6 @@
 - 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)
-- Regions can be output as tiled tiffs with scale or geospatial metadata (#594)
 
 ### Improvements
 - More untiled tiff files are handles by the bioformats reader (#569)

sources/gdal/large_image_source_gdal/__init__.py

@@ -35,7 +35,8 @@
 from large_image import config
 from large_image.cache_util import LruCacheMetaclass, methodcache, CacheProperties
 from large_image.constants import (
-    SourcePriority, TileInputUnits, TileOutputMimeTypes, TILE_FORMAT_NUMPY, TILE_FORMAT_PIL)
+    SourcePriority, TileInputUnits, TileOutputMimeTypes,
+    TILE_FORMAT_NUMPY, TILE_FORMAT_PIL, TILE_FORMAT_IMAGE)
 from large_image.exceptions import TileSourceException
 from large_image.tilesource import FileTileSource
 
@@ -1064,6 +1065,64 @@ def _encodeTiledImageFromVips(self, vimg, iterInfo, image, **kwargs):
             raise exc
         return pathlib.Path(outputPath), TileOutputMimeTypes['TILED']
 
+    def getRegion(self, format=(TILE_FORMAT_IMAGE, ), **kwargs):
+        """
+        Get a rectangular region from the current tile source.  Aspect ratio is
+        preserved.  If neither width nor height is given, the original size of
+        the highest resolution level is used.  If both are given, the returned
+        image will be no larger than either size.
+
+        :param format: the desired format or a tuple of allowed formats.
+            Formats are members of (TILE_FORMAT_PIL, TILE_FORMAT_NUMPY,
+            TILE_FORMAT_IMAGE).  If TILE_FORMAT_IMAGE, encoding may be
+            specified.
+        :param kwargs: optional arguments.  Some options are region, output,
+            encoding, jpegQuality, jpegSubsampling, tiffCompression, fill.  See
+            tileIterator.
+        :returns: regionData, formatOrRegionMime: the image data and either the
+            mime type, if the format is TILE_FORMAT_IMAGE, or the format.
+        """
+        if not isinstance(format, (tuple, set, list)):
+            format = (format, )
+        # The tile iterator handles determining the output region
+        iterInfo = self._tileIteratorInfo(**kwargs)
+        # Only use gdal.Warp of the original image if the region has not been
+        # styled.
+        useGDALWarp = (
+            iterInfo and
+            not self._jsonstyle and
+            TILE_FORMAT_IMAGE in format and
+            kwargs.get('encoding') == 'TILED')
+        if not useGDALWarp:
+            return super().getRegion(format, **kwargs)
+        srs = self.projection or self.getProj4String()
+        tl = self.pixelToProjection(
+            iterInfo['region']['left'], iterInfo['region']['top'], iterInfo['level'])
+        br = self.pixelToProjection(
+            iterInfo['region']['right'], iterInfo['region']['bottom'], iterInfo['level'])
+        outWidth = iterInfo['output']['width']
+        outHeight = iterInfo['output']['height']
+        gdalParams = large_image.tilesource.base._gdalParameters(
+            defaultCompression='lzw', **kwargs)
+        gdalParams += [
+            '-t_srs', srs,
+            '-te', str(tl[0]), str(br[1]), str(br[0]), str(tl[1]),
+            '-ts', str(int(math.floor(outWidth))), str(int(math.floor(outHeight))),
+        ]
+        fd, outputPath = tempfile.mkstemp('.tiff', 'tiledGeoRegion_')
+        os.close(fd)
+        try:
+            self._logger.info('Using gdal warp %r' % gdalParams)
+            ds = gdal.Open(self._path, gdalconst.GA_ReadOnly)
+            gdal.Warp(outputPath, ds, options=gdalParams)
+        except Exception as exc:
+            try:
+                os.unlink(outputPath)
+            except Exception:
+                pass
+            raise exc
+        return pathlib.Path(outputPath), TileOutputMimeTypes['TILED']
+
 
 def open(*args, **kwargs):
     """

test/test_source_gdal.py

@@ -435,3 +435,64 @@ def testGetTiledRegion16Bit():
     assert tileMetadata['bounds']['ymin'] == pytest.approx(3899358, 1)
     assert '+proj=merc' in tileMetadata['bounds']['srs']
     region.unlink()
+
+
+def testGetTiledRegionWithStyle():
+    imagePath = datastore.fetch('landcover_sample_1000.tif')
+    ts = large_image_source_gdal.open(imagePath, style='{"bands":[]}')
+    region, _ = ts.getRegion(output=dict(maxWidth=1024, maxHeight=1024),
+                             encoding='TILED')
+    result = large_image_source_gdal.open(str(region))
+    tileMetadata = result.getMetadata()
+    assert tileMetadata['bounds']['xmax'] == pytest.approx(2006547, 1)
+    assert tileMetadata['bounds']['xmin'] == pytest.approx(1319547, 1)
+    assert tileMetadata['bounds']['ymax'] == pytest.approx(2658548, 1)
+    assert tileMetadata['bounds']['ymin'] == pytest.approx(2149548, 1)
+    assert '+proj=aea' in tileMetadata['bounds']['srs']
+    region.unlink()
+
+
+def testGetTiledRegionWithProjectionAndStyle():
+    imagePath = datastore.fetch('landcover_sample_1000.tif')
+    ts = large_image_source_gdal.open(imagePath, projection='EPSG:3857', style='{"bands":[]}')
+    # This gets the whole world
+    region, _ = ts.getRegion(output=dict(maxWidth=1024, maxHeight=1024),
+                             encoding='TILED')
+    result = large_image_source_gdal.open(str(region))
+    tileMetadata = result.getMetadata()
+    assert tileMetadata['bounds']['xmax'] == pytest.approx(20037508, 1)
+    assert tileMetadata['bounds']['xmin'] == pytest.approx(-20037508, 1)
+    assert tileMetadata['bounds']['ymax'] == pytest.approx(20037508, 1)
+    assert tileMetadata['bounds']['ymin'] == pytest.approx(-20037508, 1)
+    assert '+proj=merc' in tileMetadata['bounds']['srs']
+    region.unlink()
+
+    # Ask for a smaller part
+    region, _ = ts.getRegion(
+        output=dict(maxWidth=1024, maxHeight=1024),
+        region=dict(left=-8622811, right=-8192317, bottom=5294998,
+                    top=5477835, units='projection'),
+        encoding='TILED')
+    result = large_image_source_gdal.open(str(region))
+    tileMetadata = result.getMetadata()
+    assert tileMetadata['bounds']['xmax'] == pytest.approx(-8192215, 1)
+    assert tileMetadata['bounds']['xmin'] == pytest.approx(-8622708, 1)
+    assert tileMetadata['bounds']['ymax'] == pytest.approx(5477783, 1)
+    assert tileMetadata['bounds']['ymin'] == pytest.approx(5294946, 1)
+    assert '+proj=merc' in tileMetadata['bounds']['srs']
+    region.unlink()
+
+
+def testGetTiledRegion16BitWithStyle():
+    imagePath = datastore.fetch('region_gcp.tiff')
+    ts = large_image_source_gdal.open(imagePath, style='{"bands":[]}')
+    region, _ = ts.getRegion(output=dict(maxWidth=1024, maxHeight=1024),
+                             encoding='TILED')
+    result = large_image_source_gdal.open(str(region))
+    tileMetadata = result.getMetadata()
+    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 '+proj=merc' in tileMetadata['bounds']['srs']
+    region.unlink()