Re: Strange output when using "min" resampling method with WarpedVRT


Sean Gillies
 

Hi Daniel,

There was a bug recently reported at https://github.com/OSGeo/gdal/issues/4098 which might be related. I'll see about reproducing what you see. Could you read the GDAL ticket and see if anything in it is relevant?

On Tue, Aug 24, 2021 at 10:23 AM <daniel.mannarino@...> wrote:
Hi folks!

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 [1]: import rasterio

 

In [2]: res_file_16 = rasterio.open('resampled_2048th_min.tif')

 

In [3]: res_data_16 = res_file_16.read([1, 2, 3], masked=True)

 

In [4]: print(res_data_16)

[[[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- 2202 2202 -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- 2202 -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- 2202 2202]

  [-- -- -- -- -- -- -- -- -- -- -- -- 2330 3258 -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]]

 

 [[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]]

 

 [[-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]

  [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --]]]


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



--
Sean Gillies

Join main@rasterio.groups.io to automatically receive all group messages.