main@rasterio.groups.io

Rasterio 1.2.8

Sean Gillies

#869

Hi all,

We had a regression in 1.2.7 affecting the transform module's xy and rowcol functions and, by introducing randomization of test execution order, found and fixed a long-lurking bug in warp's transform function. See https://github.com/mapbox/rasterio/blob/master/CHANGES.txt#L14 for issue links.

--

Sean Gillies

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

Sean Gillies

#868

Hi Daniel,

Yes, please do create an issue. One with some code that will reproduce the problem would be perfect.

Show quoted text

On Wed, Sep 8, 2021 at 11:35 AM <daniel.mannarino@...> wrote:

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
_,_._,

--

Sean Gillies

Re: Help with saving images using lossless jpeg compression

Sean Gillies

#867

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.

Show quoted text

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@...

#866

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

#865

Hi Ghedia,

Try Deflate or LZW compression.

https://gdal.org/drivers/raster/gtiff.html

Show quoted text

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

#864

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

#863

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@...

#862

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@...

#861

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]]]

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]]

[[ 1 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

#860

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

#858

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

Show quoted text

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

#857

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?

Show quoted text

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@...

#856

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., 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.]]]),

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 --]

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

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

#855

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.

Show quoted text

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

#854

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.

Show quoted text

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@...

#853

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

#852

Dear Sean,

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

SDe

Show quoted text

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

#851

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

backtrace.txt

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

Sean Gillies

#850

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

#849

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.