I have a 16x16 TIFF with 3 uint16 bands and no-data values of 0. There are a few non-no-data values in band 1, but bands 2 and 3 are all 0s (no-data). When I explicitly scale by 0.5 (i.e. effectively "zooming out" one zoom level) with rasterio.warp.reproject using either the "nearest" or "min" resampling method, the result is as expected. But when I implicitly scale by reading from a WarpedVRT with shape_out=(3, 8, 8) and resampling method "min", unexpected values of 1 appear in the the resulting bands 2 and 3. After further testing, it seems that they appear even when NOT scaling, reading into a shape_out of (3, 16, 16). This problem does not occur with a WarpedVRT and a resampling method of "nearest". Am I using the WarpedVRT incorrectly, or is this a bug?
The anomalous data was initially noticed in Linux with GDAL 3.2.0 but confirmed on MacOS Catalina with GDAL 3.2.0, 3.2.3 and 3.3.1 using rasterio 1.1.8, 1.2.3 and 1.2.6. Some commands and their output (performed with GDAL 3.3.1, rasterio 1.2.6 and Python 3.8.11):
In [6]: print(res_file_16.profile)
{'driver': 'GTiff', 'dtype': 'uint16', 'nodata': 0.0, 'width': 16, 'height': 16, 'count': 3, 'crs': CRS.from_epsg(3857), 'transform': Affine(39135.75848225398, 0.0, -2504688.542848654,
0.0, -39135.75848225398, 3757032.814319771), 'blockxsize': 16, 'blockysize': 16, 'tiled': True, 'interleave': 'pixel'}
Now let's use reproject with resampling method "min" to scale by 0.5:
In [8]: import numpy as np
In [9]: from affine import Affine
In [10]: from rasterio.enums import Resampling
In [11]: from rasterio.warp import reproject
In [12]: dst = np.zeros((3, 8, 8))
...: dst_transform = res_file_16.transform * Affine.scale(0.5)
...: reproject(
...: res_data_16,
...: dst,
...: src_transform=res_file_16.transform,
...: src_crs=res_file_16.crs,
...: dst_crs=res_file_16.crs,
...: resampling=Resampling.min,
...: src_nodata=0
...: )
Out[12]:
(array([[[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 2202., 0.],
[ 0., 0., 0., 0., 0., 0., 2330., 2202.],
[ 0., 0., 0., 0., 0., 0., 0., 0.]],
[[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.]],
[[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0., 0., 0.]]]),
Affine(78271.51696450796, 0.0, -2504688.542848654,
0.0, -78271.51696450796, 3757032.814319771))
Looks good! Exactly what I would hope for. A cluster of three 2202 values became one 2202, a cluster of two of them became one 2202, and the min of 2330 and 3258 became 2330. And of course everything in bands 2 and 3 is still masked.
Now let's try making a VRT, creating a WarpedVRT on top of it, and reading into an out_shape of (3, 8, 8) to implicitly scale it:
daniel.mannarino$ gdalbuildvrt resampled_2048th_min.vrt resampled_2048th_min.tif
0...10...20...30...40...50...60...70...80...90...100 - done.
In [15]: from rasterio.vrt import WarpedVRT
In [16]: with rasterio.open("resampled_2048th_min.vrt") as src:
...: with WarpedVRT(src, resampling=Resampling.min) as vrt:
...: res_vrt_8 = vrt.read([1, 2, 3], out_shape=(3, 8, 8), masked=True)
...:
In [17]: print(res_vrt_8)
[[[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- 2202 --]
[-- -- -- -- -- -- 3258 --]
[-- -- -- -- -- -- -- --]]
[[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- 1 --]
[-- -- -- -- -- -- 1 --]
[-- -- -- -- -- -- -- --]]
[[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- --]
[-- -- -- -- -- -- 1 --]
[-- -- -- -- -- -- 1 --]
[-- -- -- -- -- -- -- --]]]
Oye! What are those 1s doing in bands 2 and 3? Also, where is the other value of 2202 (relative to the reproject output)?
Let's see what comes out if we read from the VRT WITHOUT scaling, for the heck of it:
In [18]: with rasterio.open("resampled_2048th_min.vrt") as src:
...: with WarpedVRT(src, resampling=Resampling.min) as vrt:
...: res_vrt_16 = vrt.read([1, 2, 3], out_shape=(3, 16, 16), masked=True)
...: print(res_vrt_16)
[[[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- 2202 2202 -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- 2202 -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- 2202 2202]
[-- -- -- -- -- -- -- -- -- -- -- -- 2330 3258 -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]]
[[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- 1 1 -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- 1 -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- 1 1]
[-- -- -- -- -- -- -- -- -- -- -- -- 1 1 -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]]
[[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- 1 1 -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- 1 -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- 1 1]
[-- -- -- -- -- -- -- -- -- -- -- -- 1 1 -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]
[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]]]
The band 1 output looks fine, but bands 2 and 3 have those values of 1! Where did they come from?
Here is some further info for debugging:
daniel.mannarino$ gdalinfo -stats resampled_2048th_min.tif
Driver: GTiff/GeoTIFF
Files: resampled_2048th_min.tif
resampled_2048th_min.tif.aux.xml
Size is 16, 16
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["Popular Visualisation Pseudo-Mercator",
METHOD["Popular Visualisation Pseudo Mercator",
ID["EPSG",1024]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Web mapping and visualisation."],
AREA["World between 85.06°S and 85.06°N."],
BBOX[-85.06,-180,85.06,180]],
ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-2504688.542848654091358,3757032.814319770783186)
Pixel Size = (39135.758482253979309,-39135.758482253979309)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (-2504688.543, 3757032.814) ( 22d30' 0.00"W, 31d57' 7.78"N)
Lower Left (-2504688.543, 3130860.679) ( 22d30' 0.00"W, 27d 3'32.85"N)
Upper Right (-1878516.407, 3757032.814) ( 16d52'30.00"W, 31d57' 7.78"N)
Lower Right (-1878516.407, 3130860.679) ( 16d52'30.00"W, 27d 3'32.85"N)
Center (-2191602.475, 3443946.746) ( 19d41'15.00"W, 29d32' 6.83"N)
Band 1 Block=16x16 Type=UInt16, ColorInterp=Gray
Min=2202.000 Max=3258.000
Minimum=2202.000, Maximum=3258.000, Mean=2371.143, StdDev=364.742
NoData Value=0
Metadata:
STATISTICS_MAXIMUM=3258
STATISTICS_MEAN=2371.1428571429
STATISTICS_MINIMUM=2202
STATISTICS_STDDEV=364.74156352583
STATISTICS_VALID_PERCENT=2.734
Band 2 Block=16x16 Type=UInt16, ColorInterp=Undefined
ERROR 1: resampled_2048th_min.tif, band 2: Failed to compute statistics, no valid pixels found in sampling.
NoData Value=0
Metadata:
STATISTICS_VALID_PERCENT=0
Band 3 Block=16x16 Type=UInt16, ColorInterp=Undefined
ERROR 1: resampled_2048th_min.tif, band 3: Failed to compute statistics, no valid pixels found in sampling.
NoData Value=0
Metadata:
STATISTICS_VALID_PERCENT=0
Any ideas? Happy to file a Github Issue or similar if people think it's a bug.
Thanks,
Daniel Mannarino
