Date   

Rasterio 1.2.6

Sean Gillies
 

Hi all,

If you've been unable to install rasterio in --non-binary mode today, please see https://github.com/mapbox/rasterio/blob/master/CHANGES.txt#L14-L18 for details on the fix in version 1.2.6.

--
Sean Gillies


Re: Rasterio 1.2.5

Sean Gillies
 

Hi all,

I messed up the 1.2.5 source distribution (the .tar.gz) and there is a breakpoint set in it. I don't see any better fix than to release a 1.2.6 version later today.

The 1.2.5 wheels are fine. The mistake only affects the source distribution.


On Tue, Jun 22, 2021 at 10:52 AM Sean Gillies via groups.io <sean=mapbox.com@groups.io> wrote:
Hi all,

Rasterio 1.2.5 distributions were uploaded to PyPI yesterday evening. Here are the changes: https://github.com/mapbox/rasterio/blob/master/CHANGES.txt#L17-L21.

Thank you for the bug reports and contributions, everyone!

--
Sean Gillies



--
Sean Gillies


Rasterio 1.2.5

Sean Gillies
 

Hi all,

Rasterio 1.2.5 distributions were uploaded to PyPI yesterday evening. Here are the changes: https://github.com/mapbox/rasterio/blob/master/CHANGES.txt#L17-L21.

Thank you for the bug reports and contributions, everyone!

--
Sean Gillies


Re: Shape of resulting merged image is not what I expect

Sean Gillies
 

Hi,

I'm happy report that you won't have this problem with version 1.2.5, which is coming soon to PyPI.

As a part of the work to fix the bug noted at https://github.com/mapbox/rasterio/blob/maint-1.2/CHANGES.txt#L10, we replaced a call to math.ceil() to a call of round() in a key place and now the output shape of merge is what you expect.


On Thu, Jun 10, 2021 at 9:02 AM <juanjo.gomeznavarro@...> wrote:

[Edited Message Follows]
[Reason: typos]

I'm not sure if this is a bug or a subtle unexpected behaviour of the function rasterio.merge.merge.

I'm creating a mosaic with several images for the whole globe. The resolution is such that the image should be 16000x8001 (width x height). The calculation of the resolution is done through these variables:
pixels = 16000
res = 360 / pixels
xmin = -180 - res / 2
xmax = 180 - res / 2
ymin = -90 - res / 2
ymax = 90 + res / 2

Which result in xmin, xmax, ymin, ymax = (-180.01125, 179.98875, -90.01125, 90.01125). OK, the problem is that when I use merge as:  
mosaic, out_trans = merge(files, bounds=(xmin, ymin, xmax, ymax), res=res,  method='max')
The shape os mosaic is NOT (4, 16000, 8001), but (4, 16000, 8002) (note the extra row).

I think the problem comes from a rounding issue. In my python terminal (ymax-ymin)/res = 8001.000000000001 (whereas (xmax-xmin)/res = 16000.0). This rounding error of 1E-12 seems to be forcing the number of rows to be 8002 instead of 8001. This might seem like a tiny approximation error, but the problem is that controlling the shape of the image is very important for further calculations downstream of the algorithm.

Is this what it is supposed to happen? Is there a way to "force" merge to produce an image of the expected shape 16000x8001?

Thanks.
  
--
Sean Gillies


Re: Proposal: Allow other cloud object store providers

Sean Gillies
 

Hi,

Thanks for bringing this up! To be sure that I understand, are you proposing that potentially any cloud platform would be made first class in rasterio, concretely, as in GDAL and as with rasterio/S3 today? AWS is special in rasterio because it's what I use at work and is where most of the important public raster datasets were hosted 2-3 years ago. GDAL makes all cloud platforms first class because people or organizations paid the maintainer to do it and because it's a less complex approach than making GDAL extendable.  

What would you think if rasterio were to take the opposite approach and require users to write ~5 lines of code themselves to adapt output of, say, keystonev3 and swiftclient to a standard interface in rasterio?

On Sun, Jun 13, 2021 at 10:03 PM <ashley.sommer@...> wrote:

Hi Everyone,

I understand that AWS S3 is by far the most common cloud object store provider, especially in the geospatial community. Unfortunately some users don't get to choose which object store provider they're given to use.

Rasterio already has Session support for five different Cloud sessions, including GSSession, OSSSession, AzureSession, SwiftSession and AWSSession. (In this case, a sessions are "classes that configure access to secured resources".)

However it seems that only AWS is a first-class citizen in rasterIO for the following reasons:

  • In the docs, on the topic of the use of objects on the cloud as virtual files, only support for AWS S3 is shown: https://rasterio.readthedocs.io/en/latest/topics/vsi.html
    • Other cloud providers are not even mentioned. I didn't know other sessions existed until I looked into the code.
  • When installing rasterio, you can add support for credentailizing AWS S3 (with boto3) using the package "extras" syntax: `pip install rasterio[s3]`.
    • No other cloud providers have their own addon available at install-time.
  • When configuring a rasterio context-manger with `rasterio.Env()` you have the ability to pass in an AWSSession or aws credentials to credentialize your context.
    • If you don't pass in credentails, `Env` will make a session using `Session.aws_or_dummy()`, but doesn't attempt to check the if other provider sessions should be used.

I propose the change to move the current object-store feature implementation to be less AWS focused, with
  • Ability to install requirement packages for other cloud providers, using additional "extras" options:
    • eg, `rasterio[gs]`, `rasterio[azure]`, `rasterio[swift]`
  • Modify the `session.Env()` context manger to be cloud-platform agnostic, ie, no "aws-or-dummy" behaviour, based on which supporting-packages are installed and which environment variables are present.
  • Document all of the different cloud providers that Rasterio supports, and how to configure and use them.

I understand this is a large piece of work, and an tall proposal for an opensource project. I'm usually one to scratch my own itch, my specific requirement out of this is to be able to use rasterio with my existing OpenStack Swift object store. I want to be able to use swift by passing pre-configured application credentials and using openstack `keystonev3` and `swiftclient` libraries to configure and credentialize the context used for GDAL. Of course I could simply open a PR with some code to patch in that support in the existing SwiftSession, but I think a cleaner solution would involve looking at the issues highlighted above, to make a more generalized fix overall, on top of the extra Swift features that I require.

Thanks for reading.



--
Sean Gillies


Re: How are calculate_default_transform and aligned_target different?

Sean Gillies
 


Hi Yon,

The intent is that rasterio.warp.reproject, given the values output from calculate_default_transform, will produce a result close to that of GDAL's gdalwarp utility *without* its "-tap" option. And that if you wrap the output of calculate_default_transform with aligned_target, adding a resolution, you will get a result close to that of gdalwarp *with* its "-tap" option.

You can certainly do your own warping in chunks, using Python concurrency features, but you may also rely on GDAL's warper for this. rasterio.warp.reproject will take rasterio datasets as inputs and output and can turn over all the details to GDAL. You can set the number of threads and a limit for memory utilization and GDAL's warper will handle the chunking.

Yours,

On Thu, Jun 10, 2021 at 3:12 PM <yon.davies@...> wrote:
Greetings everyone,

I am trying to do reprojection in chunks and have been using rasterio.warp.calculate_default_transform for getting the projected extent and transform of the chunk. However, I noticed rasterio.warp.aligned_target also returns a transform and bounds, but the transform is a bit different from the one returned by calculate_default_transform. The documentation is a little vague – calculate_default_transform: Output dimensions and transform for a reprojection; aligned_target: Aligns target to specified resolution.

Can anyone illuminate this for me? And do you have any insight for which might be better for reprojecting a huge raster in chunks that will eventually be merged later? Thank you.


--
Sean Gillies


Re: Rasters with mixed dtypes

Sean Gillies
 

Hi Amaury,

Sorry for the late reply. Rasterio depends on GDAL and support for mixed data types in that library varies across the different format drivers. GeoTIFF is constrained, but others are not. I'm not aware of any comprehensive report of which formats do allow mixed types and which don't. HDF5 does, at least.

Yours,

On Thu, Jun 10, 2021 at 8:03 AM <amaury.dehecq@...> wrote:

Hello everyone,

I was wondering if there are any good examples/reasons for having rasters that contain bands with mixed dtypes e.g. ('uint8', 'float32').

This seem to be mostly supported in rasterio, but functions like rasterio.open accept only one dtype. Furthermore, mixed dtypes are not supported in GTiff apparently: https://lists.osgeo.org/pipermail/gdal-dev/2010-August/025657.html

I'm just wondering if I should make the effort to accept mixed dtypes in a library I'm developing, which is based upon rasterio.

Thanks!

Amaury


--
Sean Gillies


Proposal: Allow other cloud object store providers

ashley.sommer@...
 

Hi Everyone,

I understand that AWS S3 is by far the most common cloud object store provider, especially in the geospatial community. Unfortunately some users don't get to choose which object store provider they're given to use.

Rasterio already has Session support for five different Cloud sessions, including GSSession, OSSSession, AzureSession, SwiftSession and AWSSession. (In this case, a sessions are "classes that configure access to secured resources".)

However it seems that only AWS is a first-class citizen in rasterIO for the following reasons:

  • In the docs, on the topic of the use of objects on the cloud as virtual files, only support for AWS S3 is shown: https://rasterio.readthedocs.io/en/latest/topics/vsi.html
    • Other cloud providers are not even mentioned. I didn't know other sessions existed until I looked into the code.
  • When installing rasterio, you can add support for credentailizing AWS S3 (with boto3) using the package "extras" syntax: `pip install rasterio[s3]`.
    • No other cloud providers have their own addon available at install-time.
  • When configuring a rasterio context-manger with `rasterio.Env()` you have the ability to pass in an AWSSession or aws credentials to credentialize your context.
    • If you don't pass in credentails, `Env` will make a session using `Session.aws_or_dummy()`, but doesn't attempt to check the if other provider sessions should be used.

I propose the change to move the current object-store feature implementation to be less AWS focused, with
  • Ability to install requirement packages for other cloud providers, using additional "extras" options:
    • eg, `rasterio[gs]`, `rasterio[azure]`, `rasterio[swift]`
  • Modify the `session.Env()` context manger to be cloud-platform agnostic, ie, no "aws-or-dummy" behaviour, based on which supporting-packages are installed and which environment variables are present.
  • Document all of the different cloud providers that Rasterio supports, and how to configure and use them.

I understand this is a large piece of work, and an tall proposal for an opensource project. I'm usually one to scratch my own itch, my specific requirement out of this is to be able to use rasterio with my existing OpenStack Swift object store. I want to be able to use swift by passing pre-configured application credentials and using openstack `keystonev3` and `swiftclient` libraries to configure and credentialize the context used for GDAL. Of course I could simply open a PR with some code to patch in that support in the existing SwiftSession, but I think a cleaner solution would involve looking at the issues highlighted above, to make a more generalized fix overall, on top of the extra Swift features that I require.

Thanks for reading.


How are calculate_default_transform and aligned_target different?

yon.davies@...
 

Greetings everyone,

I am trying to do reprojection in chunks and have been using rasterio.warp.calculate_default_transform for getting the projected extent and transform of the chunk. However, I noticed rasterio.warp.aligned_target also returns a transform and bounds, but the transform is a bit different from the one returned by calculate_default_transform. The documentation is a little vague – calculate_default_transform: Output dimensions and transform for a reprojection; aligned_target: Aligns target to specified resolution.

Can anyone illuminate this for me? And do you have any insight for which might be better for reprojecting a huge raster in chunks that will eventually be merged later? Thank you.


Shape of resulting merged image is not what I expect

juanjo.gomeznavarro@...
 
Edited

I'm not sure if this is a bug or a subtle unexpected behaviour of the function rasterio.merge.merge.

I'm creating a mosaic with several images for the whole globe. The resolution is such that the image should be 16000x8001 (width x height). The calculation of the resolution is done through these variables:
pixels = 16000
res = 360 / pixels
xmin = -180 - res / 2
xmax = 180 - res / 2
ymin = -90 - res / 2
ymax = 90 + res / 2

Which result in xmin, xmax, ymin, ymax = (-180.01125, 179.98875, -90.01125, 90.01125). OK, the problem is that when I use merge as:  
mosaic, out_trans = merge(files, bounds=(xmin, ymin, xmax, ymax), res=res,  method='max')
The shape os mosaic is NOT (4, 16000, 8001), but (4, 16000, 8002) (note the extra row).

I think the problem comes from a rounding issue. In my python terminal (ymax-ymin)/res = 8001.000000000001 (whereas (xmax-xmin)/res = 16000.0). This rounding error of 1E-12 seems to be forcing the number of rows to be 8002 instead of 8001. This might seem like a tiny approximation error, but the problem is that controlling the shape of the image is very important for further calculations downstream of the algorithm.

Is this what it is supposed to happen? Is there a way to "force" merge to produce an image of the expected shape 16000x8001?

Thanks.



 


Rasters with mixed dtypes

amaury.dehecq@...
 

Hello everyone,

I was wondering if there are any good examples/reasons for having rasters that contain bands with mixed dtypes e.g. ('uint8', 'float32').

This seem to be mostly supported in rasterio, but functions like rasterio.open accept only one dtype. Furthermore, mixed dtypes are not supported in GTiff apparently: https://lists.osgeo.org/pipermail/gdal-dev/2010-August/025657.html

I'm just wondering if I should make the effort to accept mixed dtypes in a library I'm developing, which is based upon rasterio.

Thanks!

Amaury


Help with inconsistent errors using WarpedVRT/rio_tiler

fgassert@...
 

Hello,

I have been running into some puzzling errors with windowed reads of COGs in Mollweide crs.

The read raises a CPLE_AppDefinedError on transforming the window to the dataset crs the first five (5) times it is called, then it succeeds. This implies to me that there might be some behavior that is affected by caching, or with an error suppression logic in rio core. I believe the errors are raised in the first place because the COG bounds become (inf, inf, inf, inf) when transformed from Mollweide to WGS84. However, this does not affect the eventual success of the reads. The error occurs in `rasterio._base()._transform()`.

Here is a simple script to reproduce the errors.
```
import rasterio as rio
from rio_tiler.io import COGReader

arr = None
while arr is None:
    try:
        with COGReader('https://storage.googleapis.com/francis_indicators/backup/HFP/HFP_source/2013_meris_lc.tif') as cog:
            reader = cog.part((10,10,20,20))
            arr = reader.data
    except Exception as e:
        print(e)

print(arr.shape)
```

Output:
```
acos/asin: |arg| >1.+1e-14
acos/asin: |arg| >1.+1e-14
acos/asin: |arg| >1.+1e-14
acos/asin: |arg| >1.+1e-14
Reprojection failed, err = -19, further errors will be suppressed on the transform object.
(1, 1013, 1013)
```

Versions
```
rasterio: 1.2.4
rio_tiler: 2.1.0
gdal: 3.3.0
```

Thanks,
Francis


Re: Rasterio 1.2.4

Even Rouault
 

Hi Sean

regarding the change in web mercator, I guess this is due to a subtle change where we don't check anymore if latitude is close enough of PI/2 to decide the point doesn't project. Now we only test lat < PI/2. And as 90° converted in radian is slightly < PI/2, it now projects to a finite (although very big) northing value. So this changes the aspect of how a raster with world wide bounds in EPSG:4326 reprojects in EPSG:3857. Neither the previous (which truncated a few pixels) nor the new behavior (which results in an image with a big height/width ratio) are ideal, but I don't think there isn't an ideal solution given the characteristics of the mercator projection. The commit that modified this is https://github.com/OSGeo/PROJ/commit/5d22d137b1be282bdc14c1d12ab8f669f58d41a6

Regarding discovery of proj.db, I can only think to https://github.com/OSGeo/PROJ/commit/9488b26ca5b7ea763db4e210b311bcf3f04233ca which might have effect on this. It should mostly be a no-op for default builds, although I see that line https://github.com/OSGeo/PROJ/commit/9488b26ca5b7ea763db4e210b311bcf3f04233ca#diff-8a6a421339db369cb2b34097f1a6fede57afb13e7e858e255ed7d6895f814540R1478 might have subtle effects if the proj_lib directory hardcoded at build time doesn't exist at runtime. Previously it will have been used even if not existing, whereas now it is not considered (and thus a proj.db in the current working directory might be used). The new behavior makes more sense to me.

Even

Le 01/06/2021 à 18:37, Sean Gillies via groups.io a écrit :
Hi all,

Here is the list of bug fixes in rasterio 1.2.4: 


The wheels on PyPI for this version now include GDAL 3.3.0 and GEOS 3.9.1. The PROJ version remains the same, 7.2.1 until we get more experience with changes to web mercator projection (I believe I see expanded limits) and discovery of the proj.db file.

Thanks for the bug reports and discussions, everyone! 

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


Rasterio 1.2.4

Sean Gillies
 

Hi all,

Here is the list of bug fixes in rasterio 1.2.4: 


The wheels on PyPI for this version now include GDAL 3.3.0 and GEOS 3.9.1. The PROJ version remains the same, 7.2.1 until we get more experience with changes to web mercator projection (I believe I see expanded limits) and discovery of the proj.db file.

Thanks for the bug reports and discussions, everyone! 

--
Sean Gillies


Re: What do I need to do to visualize Sentinel-1 data with rasterio?

Sean Gillies
 

Thanks, Sebastian!

rasterio.plot.show cannot deal with complex arrays. It looks like a library called mpmath can, might be worth looking in https://mpmath.org/gallery/.


On Wed, May 26, 2021 at 1:46 PM Yann-Sebastien Tremblay-Johnston <yanns.tremblay@...> wrote:
Looks like your data is Single Look Complex (SLC), you'll need to convert it to float. You can do this yourself with numpy, or with https://gdal.org/drivers/raster/derived.html GDAL (which rasterio uses under the hood) will compute the intensity for you. 

e.g.

src = rasterio.open('DERIVED_SUBDATASET:INTENSITY:../x.SAFE/measurement/y-001.tif')

Sebastien
   

On Wed, May 26, 2021 at 12:38 PM <earthdata@...> wrote:
I'm a first time user. I downloaded a Sentinel-1 granule. I'm now trying to visualize it via the following script:
import rasterio
from rasterio.plot import show
path = '.../x.SAFE/measurement/y-001.tiff'
src = rasterio.open(path)
show(src)
However this gives me the following error:
TypeError: Image data of dtype complex128 cannot be converted to float
What else do I need to do?

._,

--
Sean Gillies


Re: What do I need to do to visualize Sentinel-1 data with rasterio?

Yann-Sebastien Tremblay-Johnston
 

Looks like your data is Single Look Complex (SLC), you'll need to convert it to float. You can do this yourself with numpy, or with https://gdal.org/drivers/raster/derived.html GDAL (which rasterio uses under the hood) will compute the intensity for you. 

e.g.

src = rasterio.open('DERIVED_SUBDATASET:INTENSITY:../x.SAFE/measurement/y-001.tif')

Sebastien
   

On Wed, May 26, 2021 at 12:38 PM <earthdata@...> wrote:
I'm a first time user. I downloaded a Sentinel-1 granule. I'm now trying to visualize it via the following script:
import rasterio
from rasterio.plot import show
path = '.../x.SAFE/measurement/y-001.tiff'
src = rasterio.open(path)
show(src)
However this gives me the following error:
TypeError: Image data of dtype complex128 cannot be converted to float
What else do I need to do?


What do I need to do to visualize Sentinel-1 data with rasterio?

earthdata@...
 

I'm a first time user. I downloaded a Sentinel-1 granule. I'm now trying to visualize it via the following script:
import rasterio
from rasterio.plot import show
path = '.../x.SAFE/measurement/y-001.tiff'
src = rasterio.open(path)
show(src)
However this gives me the following error:
TypeError: Image data of dtype complex128 cannot be converted to float
What else do I need to do?


Re: concurrent.futures.ThreadPoolExecutor for s3 COG reads fails - CPLReleaseMutex: Error

Sean Gillies
 

Hi Darren,

I've seen this error message before, and opened an issue at https://github.com/OSGeo/gdal/issues/2278. I haven't tried to reproduce the error with GDAL 3.3 yet.

On Fri, May 7, 2021 at 2:51 PM <dweber.consulting@...> wrote:
This could be similar to https://github.com/mapbox/rasterio/issues/1686

In an AWS-Lambda python 3.7 runtime, there is an error like:
 
    CPLReleaseMutex: Error = 1 (Operation not permitted)
    [WARNING] 2021-05-07T17:38:33.315Z 3e8b4b55-9066-426c-9319-540e39e17709
    CPLE_AppDefined in HTTP response code on https://bucket-xxx-us-west-2.s3.amazonaws.com/prefix/geotiff.tif: 0
    2021-05-07 17:38:33,366 | ERROR | '/vsis3/bucket-xxx-us-west-2/prefix/geotiff.tif' not recognized as a supported file format.

The python concurrency that wraps the rasterio/gdal reader is a common pattern like:

    dataset = []
    with concurrent.futures.ThreadPoolExecutor() as executor:
        future_to_task = {}
        for task in peril_queries:
            future = executor.submit(extract_geotiff_data, task)
            future_to_task[future] = task
    for future in concurrent.futures.as_completed(future_to_task):
        try:
            data = future.result()
            dataset.extend(data)
        except Exception as exc:
            task = future_to_task[future]
            LOGGER.error("Failed Task: %s", task)
            raise
 

Where the 'extract_geotiff_data' function is basically given a task to read an s3-COG,
and that function encapsulates everything to setup a rasterio.Env with gdal env-var
settings, opening the s3 object and reading data from it.

What is causing the `CPLReleaseMutex: Error` and how is it debugged and avoided?

TIA,
Darren

--
Sean Gillies


Re: Reconciling transforms by from_origin and from_bounds

Luke
 
Edited

rasterio.transform.from_origin(westnorthxsizeysize)

Return an Affine transformation given upper left and pixel sizes.

https://rasterio.readthedocs.io/en/latest/api/rasterio.transform.html#rasterio.transform.from_origin

After changing your height calculation from bottom - top to top - bottom, for your coords I get similar transforms:

import math
from rasterio.transform import (from_origin, from_bounds, array_bounds)
from rasterio.coords import BoundingBox

bounds = BoundingBox(
left=1679712.0, bottom=5860848.0, right=1862208.0, top=6029312.0)
width = math.ceil((bounds.right - bounds.left) /
100) # resolution of 100
height = math.ceil((bounds.top - bounds.bottom) / 100) # resolution of 100

T_o = from_origin(bounds.left, bounds.top, 100, 100)

T_b = from_bounds(*bounds
, width, height)

print(T_o)
print(T_b)

print(array_bounds(height, width, T_o))
print(array_bounds(height, width, T_b))

| 100.00, 0.00, 1679712.00| | 0.00,-100.00, 6029312.00| | 0.00, 0.00, 1.00| | 100.00, 0.00, 1679712.00| | 0.00,-99.98, 6029312.00| | 0.00, 0.00, 1.00| (1679712.0, 5860812.0, 1862212.0, 6029312.0) (1679712.0, 5860848.0, 1862208.0, 6029312.0)


Reconciling transforms by from_origin and from_bounds

ayr035@...
 

I am trying to reconcile the transforms that I get by using transform.from_origin and transform.from_bounds

bounds = BoundingBox(left=1679712.0, bottom=5860848.0, right=1862208.0, top=6029312.0)
width = math.ceil((bounds.right - bounds.left) / 100) # resolution of 100
height = math.ceil((bounds.bottom - bounds.top) / 100)
# resolution of 100

T_o = transform.from_origin(bounds.left, bounds.top, width, height)
Affine(1825.0, 0.0, 1679712.0,
      0.0, -1685.0, 6029312.0)

T_b = transform.from_bounds(*bounds, width, height)
Affine(99.99780821917808, 0.0, 1679712.0,
      0.0, -99.9786350148368, 6029312.0)

The transforms seem to differ quite significantly. If we compute the array bounds, I expect to recover the original bounds
transform.array_bounds(height, width, T_o) ->
(1679712.0, 3190087.0, 5010337.0, 6029312.0)   # bottom and right bounds are wrong?
transform.array_bounds(height, width, T_b) -> (
1679712.0, 5860848.0, 1862208.0, 6029312.0)   # These are my original bounds

Am I misusing transform.from_origin?

161 - 180 of 945