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


daniel.mannarino@...
 

Hi Sean,

Thanks for taking a look! It doesn't sound to ME like that bug would account for the behavior I'm seeing (at least not the anomalous data in the previously all-zero bands), but I'm not sure. Just in case, I decided to build and install GDAL HEAD and rasterio 1.2.6 against it to test, but it doesn't seem to have solved my problem. In fact when I tried to make it easier to observe by writing a script with synthetic data rather than just the interactive session output, the problem surfaced even when just using rasterio.warp.reproject. Further, the min algorithm seems to have selected 10000 from a masked 0, 10000, 2500, and 5000! But perhaps you can spy a mistake I made. Here's the script and the output I obtain when using git GDAL from a few days ago. Note that I decided to use smaller sample data and fewer channels just to cut down on the output. Oh, and this is with Python 3.9.6 on MacOS 10.15.7 (Catalina):

# begin resampling-test.py:
import subprocess
 
import numpy as np
import rasterio
from affine import Affine
from pyproj import CRS
from rasterio.enums import Resampling
from rasterio.profiles import DefaultGTiffProfile
from rasterio.vrt import WarpedVRT
from rasterio.warp import reproject
 
 
synthetic_data = np.zeros((2, 8, 8), dtype="uint16")
synthetic_data[0][0][1] = 2500
synthetic_data[0][1][0] = 5000
synthetic_data[0][1][1] = 10000
 
print("Original synthetic data:")
print(synthetic_data)
 
crs = CRS.from_epsg(3857)
transform = Affine(
    10.0,  # a: width of a pixel
    0.0,   # b: row rotation
    -40.0, # c: x-coord of the upper-left corner of the upper-left pixel
    0.0,   # d: column rotation
    -10.0, # e: height of a pixel (typically negative)
    40.0   # y-coord of the of the upper-left corner of the upper-left pixel
)
 
synthetic_data_profile = {
    **DefaultGTiffProfile(count=synthetic_data.shape[0]),
    "crs": crs,
    "dtype": "uint16",
    "height": synthetic_data.shape[1],
    "width": synthetic_data.shape[2],
    "transform": transform
}
 
print(f"Synthetic file profile:")
print(synthetic_data_profile)
 
with rasterio.open('synthetic_data.tif', 'w', **synthetic_data_profile) as dst:
    dst.write(synthetic_data)
 
with rasterio.open('synthetic_data.tif') as src:
    synthetic_data_saved = src.read([1,2], masked=True)
 
print(
    "Synthetic data written and re-read (should be identical to original "
    "except now 0s are masked):"
)
print(synthetic_data_saved)
 
print("Now let's scale it using rasterio.warp.reproject")
 
reprojection_dst = np.zeros((2, 4, 4), dtype="uint16")
 
reproject(
    synthetic_data_saved,
    reprojection_dst,
    src_transform=transform,
    src_crs=crs.to_dict(),
    dst_crs=crs.to_dict(),
    resampling=Resampling.min,
    src_no_data=0
)
 
print("Resulting data:")
print(reprojection_dst)
 
 
print(
    "Now let's try reading through a VRT. First through a regular VRT, "
    "then through a WarpedVRT without scaling, and finally through a "
    "WarpedVRT into an out_shape that implicitly causes scaling"
)
 
print("Creating the VRT")
subprocess.run(
    ["gdalbuildvrt", "synthetic_data.vrt", "synthetic_data.tif"],
    capture_output=True
)
 
print(f"Data as read through plain VRT:")
with rasterio.open('synthetic_data.vrt') as vrt_src:
    synthetic_data_vrt = vrt_src.read([1,2], masked=True)
print(synthetic_data_vrt)
 
print(
    "Now the data as read through a WarpedVRT with Resampling mode nearest "
    "but no out_shape specified:")
with rasterio.open('synthetic_data.vrt') as vrt_src:
    with WarpedVRT(vrt_src, resampling=Resampling.nearest) as warped_vrt:
        synthetic_data_warped_nearest = warped_vrt.read([1,2], masked=True)
print(synthetic_data_warped_nearest)
 
print(
    "Now the data as read through a WarpedVRT with Resampling mode min "
    "but no out_shape specified:")
with rasterio.open('synthetic_data.vrt') as vrt_src:
    with WarpedVRT(vrt_src, resampling=Resampling.min) as warped_vrt:
        synthetic_data_warped_min = warped_vrt.read([1,2], masked=True)
print(synthetic_data_warped_min)
 
print(
    "Finally the data as read through a WarpedVRT with Resampling mode min "
    "and an out_shape that causes implicit scaling by a factor of 0.5:")
with rasterio.open('synthetic_data.vrt') as vrt_src:
    with WarpedVRT(vrt_src, resampling=Resampling.min) as warped_vrt:
        synthetic_data_warped_min_scaled = warped_vrt.read(
            [1,2], out_shape=(2, 4, 4), masked=True
        )
print(synthetic_data_warped_min_scaled)
# end resampling-test.py:

Output when I run it locally:

daniel.mannarino$ PATH=$PATH:/usr/local/Cellar/gdal/HEAD-a0ba563_3/bin/ python resampling_test.py 

Original synthetic data:

[[[    0  2500     0     0     0     0     0     0]

  [ 5000 10000     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]]

 

 [[    0     0     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]

  [    0     0     0     0     0     0     0     0]]]

Synthetic file profile:

{'driver': 'GTiff', 'interleave': 'band', 'tiled': True, 'blockxsize': 256, 'blockysize': 256, 'compress': 'lzw', 'nodata': 0, 'dtype': 'uint16', 'count': 2, 'crs': <Projected CRS: EPSG:3857>

Name: WGS 84 / Pseudo-Mercator

Axis Info [cartesian]:

- X[east]: Easting (metre)

- Y[north]: Northing (metre)

Area of Use:

- name: World between 85.06°S and 85.06°N.

- bounds: (-180.0, -85.06, 180.0, 85.06)

Coordinate Operation:

- name: Popular Visualisation Pseudo-Mercator

- method: Popular Visualisation Pseudo Mercator

Datum: World Geodetic System 1984 ensemble

- Ellipsoid: WGS 84

- Prime Meridian: Greenwich

, 'height': 8, 'width': 8, 'transform': Affine(10.0, 0.0, -40.0,

       0.0, -10.0, 40.0)}

Synthetic data written and re-read (should be identical to original except now 0s are masked):

[[[-- 2500 -- -- -- -- -- --]

  [5000 10000 -- -- -- -- -- --]

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

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

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

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

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

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

 

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

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

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

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

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

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

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

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

Now let's scale it using rasterio.warp.reproject

/Users/daniel.mannarino/Sandbox/gdal_git/venv/lib/python3.9/site-packages/pyproj/crs/crs.py:1216: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems

  return self._crs.to_proj4(version=version)

Resulting data:

[[[2500    0    0    0]

  [   0    0    0    0]

  [   0    0    0    0]

  [   0    0    0    0]]

 

 [[   1    0    0    0]

  [   0    0    0    0]

  [   0    0    0    0]

  [   0    0    0    0]]]

Now let's try reading through a VRT. First through a regular VRT, then through a WarpedVRT without scaling, and finally through a WarpedVRT into an out_shape that implicitly causes scaling

Creating the VRT

Data as read through plain VRT:

[[[-- 2500 -- -- -- -- -- --]

  [5000 10000 -- -- -- -- -- --]

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

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

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

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

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

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

 

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

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

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

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

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

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

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

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

Now the data as read through a WarpedVRT with Resampling mode nearest but no out_shape specified:

[[[-- 2500 -- -- -- -- -- --]

  [5000 10000 -- -- -- -- -- --]

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

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

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

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

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

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

 

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

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

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

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

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

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

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

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

Now the data as read through a WarpedVRT with Resampling mode min but no out_shape specified:

[[[-- 2500 -- -- -- -- -- --]

  [5000 10000 -- -- -- -- -- --]

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

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

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

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

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

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

 

 [[-- 1 -- -- -- -- -- --]

  [1 1 -- -- -- -- -- --]

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

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

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

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

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

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

Finally the data as read through a WarpedVRT with Resampling mode min and an out_shape that causes implicit scaling by a factor of 0.5:

[[[10000 -- -- --]

  [-- -- -- --]

  [-- -- -- --]

  [-- -- -- --]]

 

 [[1 -- -- --]

  [-- -- -- --]

  [-- -- -- --]

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

Thanks for any help!
Daniel

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