From d03f8e919f10b5932f279ddc5bc803f4306119d1 Mon Sep 17 00:00:00 2001 From: Denis Rykov Date: Mon, 18 Oct 2021 22:58:38 +0200 Subject: [PATCH 1/2] Add UseMaskBand option support --- telluric/base_vrt.py | 40 +++--- telluric/gdalvrt.xsd | 182 ++++++++++++++++++++++++++-- telluric/vrt.py | 54 +++++---- tests/data/raster/google_israel.vrt | 2 +- tests/data/raster/overlap2.vrt | 2 +- tests/test_georaster.py | 8 ++ tests/test_vrt.py | 1 - 7 files changed, 237 insertions(+), 52 deletions(-) diff --git a/telluric/base_vrt.py b/telluric/base_vrt.py index 0492f99b..b7d72d8b 100644 --- a/telluric/base_vrt.py +++ b/telluric/base_vrt.py @@ -88,29 +88,35 @@ def add_band(self, dtype, band_idx, color_interp, return vrtrasterband - def add_band_simplesource(self, vrtrasterband, band_idx, dtype, relative_to_vrt, - file_name, rasterxsize, rasterysize, blockxsize=None, blockysize=None, - src_rect=None, dst_rect=None, nodata=None - ): - simplesource = ET.SubElement(vrtrasterband, 'SimpleSource') - self._setup_band_simplesource(simplesource, band_idx, dtype, relative_to_vrt, file_name, - rasterxsize, rasterysize, blockxsize, blockysize, nodata) - srcrect_element = ET.SubElement(simplesource, 'SrcRect') + def add_band_complexsource( + self, vrtrasterband, band_idx, dtype, relative_to_vrt, + file_name, rasterxsize, rasterysize, blockxsize=None, blockysize=None, + src_rect=None, dst_rect=None, nodata=None + ): + complexsource = ET.SubElement(vrtrasterband, 'ComplexSource') + self._setup_band_complexsource( + complexsource, band_idx, dtype, relative_to_vrt, file_name, + rasterxsize, rasterysize, blockxsize, blockysize, nodata) + srcrect_element = ET.SubElement(complexsource, 'SrcRect') self._setup_rect(srcrect_element, src_rect.col_off, src_rect.row_off, src_rect.width, src_rect.height) - dstrect_element = ET.SubElement(simplesource, 'DstRect') + dstrect_element = ET.SubElement(complexsource, 'DstRect') self._setup_rect(dstrect_element, dst_rect.col_off, dst_rect.row_off, dst_rect.width, dst_rect.height) - return simplesource, srcrect_element, dstrect_element - - def _setup_band_simplesource(self, simplesource, band_idx, dtype, relative_to_vrt, file_name, - rasterxsize, rasterysize, blockxsize, blockysize, nodata): - sourcefilename = ET.SubElement(simplesource, 'SourceFilename') + usemaskband = ET.SubElement(complexsource, 'UseMaskBand') + usemaskband.text = 'true' + return complexsource, srcrect_element, dstrect_element + + def _setup_band_complexsource( + self, complexsource, band_idx, dtype, relative_to_vrt, file_name, + rasterxsize, rasterysize, blockxsize, blockysize, nodata + ): + sourcefilename = ET.SubElement(complexsource, 'SourceFilename') sourcefilename.attrib['relativeToVRT'] = "1" if relative_to_vrt else "0" sourcefilename.text = vsi_path(parse_path(file_name)) - sourceband = ET.SubElement(simplesource, 'SourceBand') + sourceband = ET.SubElement(complexsource, 'SourceBand') sourceband.text = str(band_idx) - sourceproperties = ET.SubElement(simplesource, 'SourceProperties') + sourceproperties = ET.SubElement(complexsource, 'SourceProperties') sourceproperties.attrib['RasterXSize'] = str(rasterxsize) sourceproperties.attrib['RasterYSize'] = str(rasterysize) if blockxsize is not None and blockysize is not None: @@ -123,7 +129,7 @@ def _setup_band_simplesource(self, simplesource, band_idx, dtype, relative_to_vr # so till we figure it out I leave it commented out # if nodata is not None: - # nodata_elem = ET.SubElement(simplesource, 'NODATA') + # nodata_elem = ET.SubElement(complexsource, 'NODATA') # nodata_elem.text = str(nodata) def _setup_rect(self, sub_element, xoff, yoff, xsize, ysize): diff --git a/telluric/gdalvrt.xsd b/telluric/gdalvrt.xsd index 6c030b76..d8c2e5e9 100644 --- a/telluric/gdalvrt.xsd +++ b/telluric/gdalvrt.xsd @@ -34,7 +34,7 @@ - + @@ -44,14 +44,36 @@ + + - - + + + + + + + + + + + + + + + + + + + + + + @@ -63,6 +85,7 @@ + @@ -149,6 +172,7 @@ + @@ -186,10 +210,9 @@ + + - - - @@ -245,6 +268,29 @@ + + + + + + + + + + + + + + + + + + + + + + + @@ -343,7 +389,8 @@ - + + @@ -436,4 +483,123 @@ - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/telluric/vrt.py b/telluric/vrt.py index 1bc698a2..bf638232 100644 --- a/telluric/vrt.py +++ b/telluric/vrt.py @@ -90,9 +90,10 @@ def wms_vrt(wms_file, bounds=None, resolution=None): band_element = vrt.add_band(data_type, bidx, band) dst_window = Window(0, 0, dst_width, dst_height) - vrt.add_band_simplesource(band_element, bidx, data_type, False, os.path.abspath(wms_file), - orig_width, orig_height, blockx, blocky, - src_window, dst_window) + vrt.add_band_complexsource( + band_element, bidx, data_type, False, os.path.abspath(wms_file), + orig_width, orig_height, blockx, blocky, + src_window, dst_window) return vrt @@ -132,9 +133,10 @@ def boundless_vrt_doc( if background is not None: src_window = Window(0, 0, background.width, background.height) dst_window = Window(0, 0, width, height) - vrt.add_band_simplesource(band_element, bidx, dtype, False, background.name, - width, height, block_shape[1], block_shape[0], - src_window, dst_window) + vrt.add_band_complexsource( + band_element, bidx, dtype, False, background.name, + width, height, block_shape[1], block_shape[0], + src_window, dst_window) src_window = Window(0, 0, src_dataset.width, src_dataset.height) xoff = (src_dataset.transform.xoff - transform.xoff) / transform.a @@ -142,9 +144,10 @@ def boundless_vrt_doc( xsize = src_dataset.width * src_dataset.transform.a / transform.a ysize = src_dataset.height * src_dataset.transform.e / transform.e dst_window = Window(xoff, yoff, xsize, ysize) - vrt.add_band_simplesource(band_element, bidx, dtype, False, src_dataset.name, - width, height, block_shape[1], block_shape[0], - src_window, dst_window, nodata=src_dataset.nodata) + vrt.add_band_complexsource( + band_element, bidx, dtype, False, src_dataset.name, + width, height, block_shape[1], block_shape[0], + src_window, dst_window, nodata=src_dataset.nodata) if all(MaskFlags.per_dataset in flags for flags in src_dataset.mask_flag_enums): mask_band = vrt.add_mask_band('Byte') @@ -154,8 +157,9 @@ def boundless_vrt_doc( xsize = src_dataset.width ysize = src_dataset.height dst_window = Window(xoff, yoff, xsize, ysize) - vrt.add_band_simplesource(mask_band, 'mask,1', 'Byte', False, src_dataset.name, - width, height, block_shape[1], block_shape[0], src_window, dst_window) + vrt.add_band_complexsource( + mask_band, 'mask,1', 'Byte', False, src_dataset.name, + width, height, block_shape[1], block_shape[0], src_window, dst_window) return vrt @@ -173,7 +177,7 @@ def raster_list_vrt(rasters, relative_to_vrt=True, nodata=None, mask_band=None): rasters : list The list of GeoRasters. relative_to_vrt : bool, optional - If True the bands simple source url will be related to the VRT file + If True the bands complex source url will be related to the VRT file nodata : int, optional If supplied is the note data value to be used mask_band: int, optional @@ -196,7 +200,7 @@ def raster_collection_vrt(fc, relative_to_vrt=True, nodata=None, mask_band=None) rasters : FeatureCollection The FeatureCollection of GeoRasters. relative_to_vrt : bool, optional - If True the bands simple source url will be related to the VRT file + If True the bands complex source url will be related to the VRT file nodata : int, optional If supplied is the note data value to be used mask_band: int, optional @@ -253,18 +257,20 @@ def max_resolution(): else: file_name = os.path.join(os.getcwd(), raster.source_file) - vrt.add_band_simplesource(band_element, band_idx, raster.dtype, relative_to_vrt, file_name, - raster.width, raster.height, - raster.block_shape(i)[1], raster.block_shape(i)[0], - src_window, dst_window) + vrt.add_band_complexsource( + band_element, band_idx, raster.dtype, relative_to_vrt, file_name, + raster.width, raster.height, + raster.block_shape(i)[1], raster.block_shape(i)[0], + src_window, dst_window) if i == mask_band: - vrt.add_band_simplesource(mask_band_elem, "mask,%s" % (mask_band + 1), - "Byte", - relative_to_vrt, - file_name, - raster.width, raster.height, - raster.block_shape(i)[1], raster.block_shape(i)[0], - src_window, dst_window) + vrt.add_band_complexsource( + mask_band_elem, "mask,%s" % (mask_band + 1), + "Byte", + relative_to_vrt, + file_name, + raster.width, raster.height, + raster.block_shape(i)[1], raster.block_shape(i)[0], + src_window, dst_window) vrt.add_metadata(items={band_names_tag: json.dumps(band_names)}) return vrt diff --git a/tests/data/raster/google_israel.vrt b/tests/data/raster/google_israel.vrt index 09c88001..67bf0465 100644 --- a/tests/data/raster/google_israel.vrt +++ b/tests/data/raster/google_israel.vrt @@ -1 +1 @@ -b'PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]3820628.4218062493,1.0,0.0,3879332.059529266,0.0,-1.0PIXELRed--root-folder--/tests/data/google.xml1Green--root-folder--/tests/data/google.xml2Blue--root-folder--/tests/data/google.xml3' \ No newline at end of file +b'PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]3820628.4218062493,1.0,0.0,3879332.059529266,0.0,-1.0PIXELRed--root-folder--/tests/data/google.xml1trueGreen--root-folder--/tests/data/google.xml2trueBlue--root-folder--/tests/data/google.xml3true' \ No newline at end of file diff --git a/tests/data/raster/overlap2.vrt b/tests/data/raster/overlap2.vrt index c17f6e36..e3526c71 100644 --- a/tests/data/raster/overlap2.vrt +++ b/tests/data/raster/overlap2.vrt @@ -1 +1 @@ -PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]1375843.2008425537,9.999898988594913,0.0,5176043.981091773,0.0,-10.000101012425429Redtests/data/raster/overlap2.tif1Greentests/data/raster/overlap2.tif2Bluetests/data/raster/overlap2.tif3tests/data/raster/overlap2.tifmask,1 \ No newline at end of file +PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]1375843.2008425537,9.999898988594913,0.0,5176043.981091773,0.0,-10.000101012425429Redtests/data/raster/overlap2.tif1trueGreentests/data/raster/overlap2.tif2trueBluetests/data/raster/overlap2.tif3truetests/data/raster/overlap2.tifmask,1true \ No newline at end of file diff --git a/tests/test_georaster.py b/tests/test_georaster.py index dd2ee559..3648400e 100644 --- a/tests/test_georaster.py +++ b/tests/test_georaster.py @@ -908,6 +908,14 @@ def rasters_for_testing_chunks(): return rasters +def test_join_use_mask_band(): + # https://github.com/OSGeo/gdal/issues/1148 + r1 = GeoRaster2.open("tests/data/raster/overlap1.tif") + r2 = GeoRaster2.open("tests/data/raster/overlap2.tif") + joined = join([r1, r2]) + assert not joined.get(r2.corners()["bl"]).mask.all() + + @pytest.mark.parametrize("raster", rasters_for_testing_chunks()) def test_chunks_with_pad(raster): shape = (256, 256) diff --git a/tests/test_vrt.py b/tests/test_vrt.py index c7fd3ed9..26068f5d 100644 --- a/tests/test_vrt.py +++ b/tests/test_vrt.py @@ -7,7 +7,6 @@ from tempfile import NamedTemporaryFile, TemporaryDirectory from telluric.util.raster_utils import build_vrt from telluric import GeoFeature, GeoRaster2, constants, FeatureCollection -from telluric.georaster import join from telluric.vrt import wms_vrt, boundless_vrt_doc, raster_list_vrt, raster_collection_vrt From 23eac72683ee33c9221eeaa5cd4f8792213e5601 Mon Sep 17 00:00:00 2001 From: Denis Rykov Date: Tue, 19 Oct 2021 11:20:15 +0200 Subject: [PATCH 2/2] Add a small comment --- tests/test_georaster.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_georaster.py b/tests/test_georaster.py index 3648400e..1287efe7 100644 --- a/tests/test_georaster.py +++ b/tests/test_georaster.py @@ -910,6 +910,7 @@ def rasters_for_testing_chunks(): def test_join_use_mask_band(): # https://github.com/OSGeo/gdal/issues/1148 + # The test checks the pixel at the coordinate of bottom left corner of the r2 raster is unmasked. r1 = GeoRaster2.open("tests/data/raster/overlap1.tif") r2 = GeoRaster2.open("tests/data/raster/overlap2.tif") joined = join([r1, r2])