Date   

Re: Help with saving images using lossless jpeg compression

Sean Gillies
 

Hi,

Lossless JPEG compression is only possible with the JPEG2000 format drivers. It's not an option for GeoTIFF.

See https://gdal.org/drivers/raster/jp2openjpeg.html#lossless-compression for the details. There are 4 other JPEG2000 format drivers for GDAL and they may have different options.


On Wed, Sep 8, 2021 at 6:47 AM Varisht Ghedia <varisht@...> wrote:
Hi,

While saving transformed and cut geotiff files to local filesystem, I wanted to save them in a lossless format. Currently, I use the "GTiff" driver to save the rasters in a tif format. While this does preserve the data, the file size is quite large (60MB~100MB). On saving them in a compressed format using JPEG a lot of quality is lost and the images lack the original details. Is there a provision to save these rasters in a different lossless format which is lighter and offers good compression? I explored lossless jpeg compression and wanted to use it. What is the right way to implement this?

Thanks and regards,
Varisht Ghedia


--
Sean Gillies


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

daniel.mannarino@...
 

Hi folks,

I'm wondering what my next steps should be for this. I tried the script I posted a week ago after installing GDAL 3.3.2 and Python 3.9.7 via Homebrew (GDAL 3.3.2 is not yet in the tree for Homebrew but I edited the recipe to point to 3.3.2 and it built without difficulty), then rasterio 1.2.7, pyproj, and affine via pip. The results are still incorrect (there are values of 1 in channels which should all be masked 0s and the value selected for the minimum is not actually the smallest). Should I create a Github issue at this point?
Thanks,
Daniel Mannarino


Re: Help with saving images using lossless jpeg compression

Idan Miara
 

Hi Ghedia, 

Try Deflate or LZW compression. 


On Wed, 8 Sep 2021, 14:47 Varisht Ghedia, <varisht@...> wrote:
Hi,

While saving transformed and cut geotiff files to local filesystem, I wanted to save them in a lossless format. Currently, I use the "GTiff" driver to save the rasters in a tif format. While this does preserve the data, the file size is quite large (60MB~100MB). On saving them in a compressed format using JPEG a lot of quality is lost and the images lack the original details. Is there a provision to save these rasters in a different lossless format which is lighter and offers good compression? I explored lossless jpeg compression and wanted to use it. What is the right way to implement this?

Thanks and regards,
Varisht Ghedia


Help with saving images using lossless jpeg compression

Varisht Ghedia
 

Hi,

While saving transformed and cut geotiff files to local filesystem, I wanted to save them in a lossless format. Currently, I use the "GTiff" driver to save the rasters in a tif format. While this does preserve the data, the file size is quite large (60MB~100MB). On saving them in a compressed format using JPEG a lot of quality is lost and the images lack the original details. Is there a provision to save these rasters in a different lossless format which is lighter and offers good compression? I explored lossless jpeg compression and wanted to use it. What is the right way to implement this?

Thanks and regards,
Varisht Ghedia


Rasterio 1.2.7

Sean Gillies
 

Hi all,

More EPSG coordinate axis order fallout means that rasterio versions 1.2.6 and older will not be fully compatible with GDAL 3.3.2. Rasterio 1.2.7 is on PyPI now and is compatible with all GDAL releases back to 2.3.0. It also includes several bug fixes. See https://github.com/mapbox/rasterio/blob/master/CHANGES.txt#L14 for the details.

The rasterio 1.2.7 wheels for AMD64 linux and macos platforms include GDAL 3.3.2 and PROJ 7.2.1.

--
Sean Gillies


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

daniel.mannarino@...
 

Oops, I mangled my message a bit during editing. To clarify, "the problem [with 1s appearing in the other channel] surfaced even when just using rasterio.warp.reproject. Further, [when reading from a WarpedVRT] the min algorithm seems to have selected 10000 from a masked 0, 10000, 2500, and 5000"


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


mask.mask() / Is prohibiting 'crop and invert' necessary?

Shinji Suzuki
 

Hello.
Specifically https://github.com/mapbox/rasterio/blob/8cb216ca83e57284f8a56bafbe8eda4334a34db6/rasterio/mask.py#L71

I don't see that they are exclusive or redundant to each other. And a
sample code I wrote works though both flags are specified to be True
once the check is removed.
Was there confusion between masking outside/inside of feature and
masking inside/outside of window when the relevant code was inserted
perhaps?

-shinji


Re: Tracking down a segmentation fault

Angus Dickey
 

Sean,

I created an issue as suggested.

https://github.com/OSGeo/gdal/issues/4318

Thanks,

Angus

Angus Dickey

Technical Lead

angus@...  |  LinkedIn


MyHEAT.ca  HEAT Maps  |  SOLAR Maps



On Fri, Aug 20, 2021 at 8:05 AM Sean Gillies via groups.io <sean=mapbox.com@groups.io> wrote:
Hi Angus,

Would you be willing to make an issue at https://github.com/OSGeo/gdal/issues. It's possible that we're using GDAL incorrectly, but it should guard against a crash like this.

On Wed, Aug 18, 2021 at 12:12 PM Angus Dickey <angus@...> wrote:

Thanks to both Sean and Even for their responses.

I ended up compiling GDAL 3.3.0 (to include the debug symbols), then building rasterio on top of that. The backtrace (full bt is attached) now provides some more information:

#0  GDALDatasetPool::_CloseDatasetIfZeroRefCount (this=0x7fffb03c4940, pszFileName=pszFileName@entry=0x7fffa85946f0 "/vsis3/bucket/path/01.cog.tiff", pszOwner=pszOwner@entry=0x7fffa83c6eb0 "0x7fffa84f38c0") at gdalproxypool.cpp:378
#1  0x00007fffcb5034a2 in GDALDatasetPool::CloseDatasetIfZeroRefCount (pszFileName=0x7fffa85946f0 "/vsis3/bucket/path/01.cog.tiff", eAccess=eAccess@entry=GA_ReadOnly, pszOwner=pszOwner@entry=0x7fffa83c6eb0 "0x7fffa84f38c0")at gdalproxypool.cpp:520

Seems the issue is here in GDAL. I am not sure if this is a bug or am am doing something "off label". For context I am reading a VRT stored in S3 that is made up of a series of COGS.

Should I be moving this out of the rasterio list to GDAL? It might not be easy to replicate with GDAL alone though.

Thanks again for any thoughts,

Angus


--
Sean Gillies


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


Strange output when using "min" resampling method with WarpedVRT

daniel.mannarino@...
 

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


Re: Tracking down a segmentation fault

Sean Gillies
 

Hi Angus,

Would you be willing to make an issue at https://github.com/OSGeo/gdal/issues. It's possible that we're using GDAL incorrectly, but it should guard against a crash like this.

On Wed, Aug 18, 2021 at 12:12 PM Angus Dickey <angus@...> wrote:

Thanks to both Sean and Even for their responses.

I ended up compiling GDAL 3.3.0 (to include the debug symbols), then building rasterio on top of that. The backtrace (full bt is attached) now provides some more information:

#0  GDALDatasetPool::_CloseDatasetIfZeroRefCount (this=0x7fffb03c4940, pszFileName=pszFileName@entry=0x7fffa85946f0 "/vsis3/bucket/path/01.cog.tiff", pszOwner=pszOwner@entry=0x7fffa83c6eb0 "0x7fffa84f38c0") at gdalproxypool.cpp:378
#1  0x00007fffcb5034a2 in GDALDatasetPool::CloseDatasetIfZeroRefCount (pszFileName=0x7fffa85946f0 "/vsis3/bucket/path/01.cog.tiff", eAccess=eAccess@entry=GA_ReadOnly, pszOwner=pszOwner@entry=0x7fffa83c6eb0 "0x7fffa84f38c0")at gdalproxypool.cpp:520

Seems the issue is here in GDAL. I am not sure if this is a bug or am am doing something "off label". For context I am reading a VRT stored in S3 that is made up of a series of COGS.

Should I be moving this out of the rasterio list to GDAL? It might not be easy to replicate with GDAL alone though.

Thanks again for any thoughts,

Angus


--
Sean Gillies


Re: Create a Dataset from an array

Sean Gillies
 

I feel it's too early to make a new function for this. First, I'd like to try to use this to eliminate the copying of data inside InMemoryRaster.


On Thu, Aug 19, 2021 at 3:39 PM <ayr035@...> wrote:
Hi Sean,

Thank you for the example code. It seems to work reliably enough for my use case. It was way easier than I imagined it being as I thought it might involve hacking the InMemoryRaster class. Would it be a desirable addition to make it easy for a user to do this in rasterio (at least a helper function generate the MEM string given an array)?

Wrapping the array in a Dataset does not seem to have much of a noticeable performance penalty if any at all. The performance of reading directly from the array (making a copy of any slices) and reading from the Dataset were almost identical.

Thanks again.


--
Sean Gillies


Re: Create a Dataset from an array

ayr035@...
 

Hi Sean,

Thank you for the example code. It seems to work reliably enough for my use case. It was way easier than I imagined it being as I thought it might involve hacking the InMemoryRaster class. Would it be a desirable addition to make it easy for a user to do this in rasterio (at least a helper function generate the MEM string given an array)?

Wrapping the array in a Dataset does not seem to have much of a noticeable performance penalty if any at all. The performance of reading directly from the array (making a copy of any slices) and reading from the Dataset were almost identical.

Thanks again.


Re: Read data from a sequence of TIF files

Shaunak De
 

Dear Sean,

Thank you for the recommendation and quick response. Was able to successfully create the VRT and read and reproject in rasterio.


SDe

On Mon, Aug 16, 2021, 6:06 PM Sean Gillies <sean.gillies@...> wrote:
Hi Shaunak,

A virtual mosaic using GDAL's VRT format https://gdal.org/drivers/raster/vrt.html is the best-tested approach for doing this. Rasterio can open and operate on such VRT XML files.

On Mon, Aug 16, 2021 at 4:58 PM <shaunakde@...> wrote:

Hello,

I was trying to use the Global Surface Water dataset found here: https://global-surface-water.appspot.com/download 

The data is delivered as a sequence of TIF files tiles as 10x10 degree granules. For example the granule starting at 10N and 10E is here: https://storage.googleapis.com/global-surface-water/downloads2020/occurrence/occurrence_120E_40Sv1_3_2020.tif 

I was wondering what the best way to read all the data intersecting an AOI was (without having to manually merge all the data)

Shaunak


--
Sean Gillies


Re: Tracking down a segmentation fault

Angus Dickey
 

Thanks to both Sean and Even for their responses.

I ended up compiling GDAL 3.3.0 (to include the debug symbols), then building rasterio on top of that. The backtrace (full bt is attached) now provides some more information:

#0  GDALDatasetPool::_CloseDatasetIfZeroRefCount (this=0x7fffb03c4940, pszFileName=pszFileName@entry=0x7fffa85946f0 "/vsis3/bucket/path/01.cog.tiff", pszOwner=pszOwner@entry=0x7fffa83c6eb0 "0x7fffa84f38c0") at gdalproxypool.cpp:378
#1  0x00007fffcb5034a2 in GDALDatasetPool::CloseDatasetIfZeroRefCount (pszFileName=0x7fffa85946f0 "/vsis3/bucket/path/01.cog.tiff", eAccess=eAccess@entry=GA_ReadOnly, pszOwner=pszOwner@entry=0x7fffa83c6eb0 "0x7fffa84f38c0")at gdalproxypool.cpp:520

Seems the issue is here in GDAL. I am not sure if this is a bug or am am doing something "off label". For context I am reading a VRT stored in S3 that is made up of a series of COGS.

Should I be moving this out of the rasterio list to GDAL? It might not be easy to replicate with GDAL alone though.

Thanks again for any thoughts,

Angus

 


Re: Question about code fragment in features.geometry_window()

Sean Gillies
 

Hi,

On Tue, Aug 17, 2021 at 11:35 AM Shinji Suzuki <shinji.suzuki@...> wrote:
Oops, I overlooked the fact that min/max was computed after (x,y) gets
transformed.
Please disregard point 2.


On Wed, Aug 18, 2021 at 2:22 AM Shinji Suzuki via groups.io
<shinji.suzuki=gmail.com@groups.io> wrote:
>
> Hi.
> Starting this line, bounds extended by (pad_x, pad_y) are computed and
> later get transformed.
> https://github.com/mapbox/rasterio/blob/master/rasterio/features.py#L451
>
> What  I don't understand are;
> 1. 'pad_x' is subtracted from 'bottom'. Shouldn't 'pad_y' be the one
> that should be related to 'bottom'?
> 2. Same value is computed twice.(e.g. left - pad_x). Since the
> computed values are used later only for obtaining max and min. There
> are no point in yielding a value multiple times.
>
> If my understanding is correct that this code is for computing the
> bounding box after transformation, should not this be changed as
> follows?
>
> diff --git a/rasterio/features.py b/rasterio/features.py
> index 768f6429..e0764e63 100644
> --- a/rasterio/features.py
> +++ b/rasterio/features.py
> @@ -451,12 +451,12 @@ def geometry_window(
>      xs = [
>          x
>          for (left, bottom, right, top) in all_bounds
> -        for x in (left - pad_x, right + pad_x, right + pad_x, left - pad_x)
> +        for x in (left - pad_x, right + pad_x, right - pad_x, left + pad_x)
>      ]
>      ys = [
>          y
>          for (left, bottom, right, top) in all_bounds
> -        for y in (top + pad_y, top + pad_y, bottom - pad_x, bottom - pad_x)
> +        for y in (top + pad_y, top - pad_y, bottom - pad_y, bottom + pad_y)
>      ]
>
>      rows1, cols1 = rowcol(
>
> Thank you for reading.

Indeed I believe you have found a bug! We're tracking it at https://github.com/mapbox/rasterio/issues/2266.

Thanks,

--
Sean Gillies


Re: Tracking down a segmentation fault

Even Rouault
 



I did a shallow search of the GDAL issue tracker and found https://github.com/OSGeo/gdal/issues/3189.

Likely unrelated: the above issue was specific to VRT pansharpening.

-- 
http://www.spatialys.com
My software is free, but my time generally not.


Re: Tracking down a segmentation fault

Sean Gillies
 

Hi Angus,

On Wed, Aug 11, 2021 at 9:41 PM Angus Dickey <angus@...> wrote:
We have a Chalice app that uses Rasterio to read COGs from an AWS S3 bucket, do some processing, and deliver the result as XYZ tiles.

Recently we upgraded the rasterio version and started seeing segmentation faults killing the chalice development server. Rasterio 1.2.3 and lower work fine, but 1.2.4 and up segfaults after a small number of requests.

The app uses a python 3.6 virtual environment and the rasterio manylinux wheels (which are awesome BTW). The backtrace (from gdb) seems to show the problem is in libgdal but doesn't provide much to go on:

Thread 13 "python" received signal SIGSEGV, Segmentation fault.
[Switching to Thread 0x7fff9e5d8700 (LWP 3887)]
0x00007fffad92e0b3 in ?? () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
(gdb) bt
#0  0x00007fffad92e0b3 in ?? () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#1  0x00007fffad92e552 in ?? () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#2  0x00007fffad92e5b1 in GDALProxyPoolDataset::~GDALProxyPoolDataset() ()
   from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#3  0x00007fffad92e679 in GDALProxyPoolDataset::~GDALProxyPoolDataset() ()
   from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#4  0x00007fffad8d0fcf in GDALDataset::ReleaseRef() () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#5  0x00007fffad86448c in VRTSimpleSource::~VRTSimpleSource() () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#6  0x00007fffad864569 in VRTComplexSource::~VRTComplexSource() () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#7  0x00007fffad863e8a in VRTSourcedRasterBand::CloseDependentDatasets() ()
   from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#8  0x00007fffad863eeb in VRTSourcedRasterBand::~VRTSourcedRasterBand() ()
   from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#9  0x00007fffad863f59 in VRTSourcedRasterBand::~VRTSourcedRasterBand() ()
   from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#10 0x00007fffad8d67ef in GDALDataset::~GDALDataset() () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#11 0x00007fffad8409f2 in VRTDataset::~VRTDataset() () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#12 0x00007fffad840a99 in VRTDataset::~VRTDataset() () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#13 0x00007fffad92dd11 in ?? () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#14 0x00007fffad92e311 in ?? () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#15 0x00007fffad92e62c in GDALProxyPoolDataset::~GDALProxyPoolDataset() ()
   from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#16 0x00007fffad92e679 in GDALProxyPoolDataset::~GDALProxyPoolDataset() ()
   from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#17 0x00007fffad8d0fcf in GDALDataset::ReleaseRef() () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#18 0x00007fffad86448c in VRTSimpleSource::~VRTSimpleSource() () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#19 0x00007fffad864569 in VRTComplexSource::~VRTComplexSource() () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#20 0x00007fffad863e8a in VRTSourcedRasterBand::CloseDependentDatasets() ()
   from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#21 0x00007fffad863eeb in VRTSourcedRasterBand::~VRTSourcedRasterBand() ()
   from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#22 0x00007fffad863f59 in VRTSourcedRasterBand::~VRTSourcedRasterBand() ()
   from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#23 0x00007fffad8d67ef in GDALDataset::~GDALDataset() () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#24 0x00007fffad8409f2 in VRTDataset::~VRTDataset() () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#25 0x00007fffad840a99 in VRTDataset::~VRTDataset() () from /home/user/project/env/lib/python3.6/site-packages/rasterio/../rasterio.libs/libgdal-41ff680c.so.29.0.0
#26 0x00007fffaec0ac99 in ?? () from /home/user/project/env/lib/python3.6/site-packages/rasterio/_base.cpython-36m-x86_64-linux-gnu.so
#27 0x00007fffaec21216 in ?? () from /home/user/project/env/lib/python3.6/site-packages/rasterio/_base.cpython-36m-x86_64-linux-gnu.so
#28 0x00007fffaec2d5d0 in ?? () from /home/user/project/env/lib/python3.6/site-packages/rasterio/_base.cpython-36m-x86_64-linux-gnu.so
#29 0x0000000000566bbc in _PyCFunction_FastCallDict ()
#30 0x00000000005a4cd1 in _PyObject_FastCallDict ()
#31 0x00000000005a4fb8 in PyObject_CallFunctionObjArgs ()
#32 0x000000000050d698 in _PyEval_EvalFrameDefault ()
#33 0x00000000005095c8 in ?? ()
#34 0x000000000050a2fd in ?? ()
#35 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#36 0x00000000005095c8 in ?? ()
#37 0x000000000050a2fd in ?? ()
#38 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#39 0x0000000000507be4 in ?? ()
#40 0x0000000000588e5c in ?? ()
#41 0x000000000059fd0e in PyObject_Call ()
#42 0x000000000050d256 in _PyEval_EvalFrameDefault ()
#43 0x00000000005095c8 in ?? ()
#44 0x000000000050a2fd in ?? ()
#45 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#46 0x0000000000507be4 in ?? ()
#47 0x0000000000509900 in ?? ()
#48 0x000000000050a2fd in ?? ()
#49 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#50 0x0000000000507be4 in ?? ()
#51 0x0000000000509900 in ?? ()
#52 0x000000000050a2fd in ?? ()
#53 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#54 0x00000000005095c8 in ?? ()
#55 0x000000000050a2fd in ?? ()
#56 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#57 0x0000000000508cd5 in _PyFunction_FastCallDict ()
#58 0x0000000000594a01 in ?? ()
#59 0x000000000054a971 in ?? ()
#60 0x00000000005a9dac in _PyObject_FastCallKeywords ()
#61 0x000000000050a433 in ?? ()
#62 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#63 0x0000000000507be4 in ?? ()
#64 0x0000000000508ec2 in _PyFunction_FastCallDict ()
#65 0x0000000000594a01 in ?? ()
#66 0x000000000054a971 in ?? ()
#67 0x00000000005a9dac in _PyObject_FastCallKeywords ()
#68 0x000000000050a433 in ?? ()
#69 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#70 0x0000000000508cd5 in _PyFunction_FastCallDict ()
#71 0x0000000000594a01 in ?? ()
#72 0x000000000054a971 in ?? ()
#73 0x00000000005a9dac in _PyObject_FastCallKeywords ()
#74 0x000000000050a433 in ?? ()
#75 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#76 0x0000000000507be4 in ?? ()
#77 0x0000000000509900 in ?? ()
#78 0x000000000050a2fd in ?? ()
#79 0x000000000050cc96 in _PyEval_EvalFrameDefault ()
#80 0x00000000005095c8 in ?? ()
#81 0x000000000050a2fd in ?? ()
#82 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#83 0x00000000005095c8 in ?? ()
#84 0x000000000050a2fd in ?? ()
#85 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#86 0x00000000005095c8 in ?? ()
#87 0x000000000050a2fd in ?? ()
#88 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#89 0x00000000005095c8 in ?? ()
#90 0x000000000050a2fd in ?? ()
#91 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#92 0x0000000000507be4 in ?? ()
#93 0x0000000000508ec2 in _PyFunction_FastCallDict ()
#94 0x0000000000594a01 in ?? ()
#95 0x0000000000549e8f in ?? ()
#96 0x00000000005515c1 in ?? ()
#97 0x00000000005a48ec in _PyObject_FastCallDict ()
#98 0x00000000005f03bc in ?? ()
#99 0x00000000005a9dac in _PyObject_FastCallKeywords ()
#100 0x000000000050a433 in ?? ()
#101 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#102 0x00000000005095c8 in ?? ()
#103 0x000000000050a2fd in ?? ()
#104 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#105 0x0000000000508cd5 in _PyFunction_FastCallDict ()
#106 0x0000000000594a01 in ?? ()
#107 0x000000000059fd0e in PyObject_Call ()
#108 0x000000000050d256 in _PyEval_EvalFrameDefault ()
#109 0x00000000005095c8 in ?? ()
#110 0x000000000050a2fd in ?? ()
#111 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#112 0x00000000005095c8 in ?? ()
#113 0x000000000050a2fd in ?? ()
#114 0x000000000050beb4 in _PyEval_EvalFrameDefault ()
#115 0x0000000000508cd5 in _PyFunction_FastCallDict ()
#116 0x0000000000594a01 in ?? ()
#117 0x000000000059fd0e in PyObject_Call ()
#118 0x00000000005e12f2 in ?? ()
#119 0x0000000000631b94 in ?? ()
#120 0x00007ffff77ca6db in start_thread (arg=0x7fff9e5d8700) at pthread_create.c:463
#121 0x00007ffff7b0371f in clone () at ../sysdeps/unix/sysv/linux/x86_64/clone.S:95


Is there a location where the libgdal version for each resterio release is documented? I presume 1.2.4 uses a newer GDAl version.

Any tips for tracking down what the problem might be in rasterio/GDAL?

Thanks for any help,

Angus


Rasterio's 1.2.3 wheels include GDAL 3.2.2. The 1.2.4 wheels include GDAL 3.3.0. We really should write these down in the readme.

I did a shallow search of the GDAL issue tracker and found https://github.com/OSGeo/gdal/issues/3189. There may be a lead in there or in another issue I didn't find.

--
Sean Gillies


Re: WarpedVRT and mixed data types

Even Rouault
 


In the second case I am surprised that the same exception is not raised. It may be because GDALCreateWarpedVRT (called by rasterio's WarpedVRT constructor) makes a mistake in determining the data types of the source bands and assigns "int32" to all of them.

The warping kernel only supports one single working data type for all bands. The  GDALWarpResolveWorkingDataType() function

https://github.com/OSGeo/gdal/blob/c6f0ea7d8a4da08e0d633eb381677c23f3934be3/gdal/alg/gdalwarper.cpp#L1439

will select the "largest" data type from input bands.

Even

--

http://www.spatialys.com
My software is free, but my time generally not.

21 - 40 of 880