diff --git a/sources/gdal/large_image_source_gdal/__init__.py b/sources/gdal/large_image_source_gdal/__init__.py index 09a18c0d0..51a729bcb 100644 --- a/sources/gdal/large_image_source_gdal/__init__.py +++ b/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): """ diff --git a/test/test_source_gdal.py b/test/test_source_gdal.py index 64281a6e5..c11337a60 100644 --- a/test/test_source_gdal.py +++ b/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()