Date   

Re: MemoryFile workflow - should closing a dataset close the memfile?

Dion Häfner <dion.haefner@...>
 

You can use [contextlib.ExitStack](https://docs.python.org/3/library/contextlib.html#contextlib.ExitStack) in situations like these to save indentation levels (or enter multiple contexts in a comprehension).

On 29/07/2019 09.45, ronipay via Groups.Io wrote:

Thanks for the reply Luke!

A context manager would work, but it's ergonomically annoying, no? I have a several functions I would like to pipe, each returning a raster. In that case I will have a lot of nested with statements:

with do_something1(ds1) as ds2:

    with do_something2(ds2) as ds3:

        with ...


Currently my solution was to change MemoryFile to NamedTemporaryFile, which solves the memory leak issue at the expense of some speed (but it's not terrible when you have an SSD).


Re: MemoryFile workflow - should closing a dataset close the memfile?

ronipay@...
 

Thanks for the reply Luke!

A context manager would work, but it's ergonomically annoying, no? I have a several functions I would like to pipe, each returning a raster. In that case I will have a lot of nested with statements:

with do_something1(ds1) as ds2:

    with do_something2(ds2) as ds3:

        with ...


Currently my solution was to change MemoryFile to NamedTemporaryFile, which solves the memory leak issue at the expense of some speed (but it's not terrible when you have an SSD).


Re: MemoryFile workflow - should closing a dataset close the memfile?

Luke
 

You could use a context manager to clean it up automatically:

from contextlib import contextmanager
 
import rasterio
from rasterio import MemoryFile
 
@contextmanager
def mem_raster(data, **profile):
    with MemoryFile() as memfile:
        with memfile.open(**profile) as dataset_writer:
            dataset_writer.write(data)
 
        with memfile.open() as dataset_reader:
            yield dataset_reader
 
 
#setup, get data from somewhere, copy or create profile etc...
 
with mem_raster(data, **profile) as ds:
    do_something_with(ds)
 
# the memfile is cleaned up after exiting the with context


MemoryFile workflow - should closing a dataset close the memfile?

ronipay@...
 

Hey all,

I have a general workflow which causes a memory leak. I have several functions which receive 1 or more rasters, perform some operations, and return a MemoryFile based raster, on which I perform other operations, which may also return a memory based raster, etc.

Closing this raster (or deleting the object) does not free the MemoryFile memory. So for now I have to propagate the memfile to all consumer functions which is annoying, not ergonomic, and prone to mistakes.

I'm familiar with rioxarray which might be useful, but we have a lot of code based on pure rasterio.

What is the best way to deal with this? Is it possible to create a specialized DatasetReader/DatasetWriter which closes the memory file on __del__?

Thanks


Re: Is it possible to open a binary file with rasterio.open?

Sean Gillies
 

Ryan,

It looks like you're trying to pass an opened file object into Python's open() function. That won't work. It only takes strings and paths.

On Mon, Jul 22, 2019 at 2:18 PM Ryan Avery <ravery@...> wrote:
Thanks for the suggestion Vincent, it looks like I still get the error 

TypeError: expected str, bytes or os.PathLike object, not _io.BytesIO\

Because this happens in the python open(image_bytes, 'rb') call not the  MemoryFile call. 

When I tried following this post to convert the BytesIO object to a bytes object like so but I don't think it actually converts it to a bytes object

with open(image_bytes.read(), 'rb') as f, MemoryFile(f) as memfile:

I get this error: ValueError: embedded null byte within the python open() call

Any other suggestions are super appreciated!




--
Sean Gillies


Re: Is it possible to open a binary file with rasterio.open?

Ryan Avery
 

Thanks for the suggestion Vincent, it looks like I still get the error 

TypeError: expected str, bytes or os.PathLike object, not _io.BytesIO\

Because this happens in the python open(image_bytes, 'rb') call not the  MemoryFile call. 

When I tried following this post to convert the BytesIO object to a bytes object like so but I don't think it actually converts it to a bytes object

with open(image_bytes.read(), 'rb') as f, MemoryFile(f) as memfile:

I get this error: ValueError: embedded null byte within the python open() call

Any other suggestions are super appreciated!



Re: Is it possible to open a binary file with rasterio.open?

vincent.sarago@...
 

Hi Ryan, 
This is just an idea but can you try with

```
with open(image_bytes, 'rb') as f, MemoryFile(f.read()) as memfile:
    with memfile.open() as src:
        arr = reshape_as_image(src.read())
```


Re: Is it possible to open a binary file with rasterio.open?

Ryan Avery
 

Hey all,

I'm trying something similar, trying to read a bytes object from a POST request, which originally comes from a tif file. Here's my code:

with open(image_bytes, 'rb') as f, MemoryFile(f) as memfile:
with memfile.open() as src:
arr = reshape_as_image(src.read())


but I get this error
```

TypeError: expected str, bytes or os.PathLike object, not _io.BytesIO

```

```
{

    "TaskId": 1828,

    "Status": "failed: <class 'TypeError'>\nTraceback (most recent call last):\n  File \"/app/keras_iNat_api/runserver.py\", line 60, in detect\n    arr_for_detection, image_for_drawing = keras_detector.open_image(image_bytes)\n  File \"./keras_detector.py\", line 32, in open_image\n    with open(image_bytes, 'rb') as f, MemoryFile(f) as memfile:\nTypeError: expected str, bytes or os.PathLike object, not _io.BytesIO\n",

    "Timestamp": "2019-07-22 18:33:06",

    "Endpoint": "uri"

 

}
```

I can' seem to find how to convert a _io.BytesIO object into a bytes object, any suggestions?

Cheers,


Re: build_overviews

James David Smith
 

A belated thanks for this Vincent !

On Fri, 12 Jul 2019 at 14:56, <vincent.sarago@...> wrote:

Hi James,
The level of overviews is usually linked to the size of the raster and its internal tiles (block).
here is an example in rio-cogeo: https://github.com/cogeotiff/rio-cogeo/blob/master/rio_cogeo/utils.py#L95-L118

Let say you have a 10 000 x 10 000 px raster with 256x256 block/tile size (to allow fast reading). The overview level is then 6 which results in `[2, 4, 8, 16, 32, 64]` levels (https://github.com/cogeotiff/rio-cogeo/blob/master/rio_cogeo/cogeo.py#L223).


Rendering a raster in a Cartopy GeoAxes object

john.o.woods@...
 

I tried to send this question once but I'm not sure it went through.

I'm trying to render a rasterio raster onto a Cartopy GeoAxes. This used to be possible, as evidenced here: https://github.com/mapbox/rasterio/issues/1723

I've also tried the method that Alan Snow kindly provided, in his response to this issue: https://github.com/mapbox/rasterio/issues/1724

In both cases, I can render Cartopy features. I cannot get the raster to show up on the same set of axes as the features, however. If I plot only the raster in such a way that it occupies the full plot, I can get the raster, but the Cartopy features won't show up.

Is it still possible to plot Cartopy features together with a raster? Can anyone point me toward a non-deprecated example?

Thank you!
John


Re: Writing raster to GeoPackage layer

Ryan
 

Thanks Sean! Here is a full example, using your suggestion, but still getting the same file format error:

import rasterio
import fiona
from shapely.geometry import Polygon, mapping
raster_fn = 'rasterio/tests/data/RGB.byte.tif'
out_fn = '/tmp/vector-and-raster.gpkg'
with rasterio.open(raster_fn) as r_ds:
    array = r_ds.read()
    l,b,r,t = r_ds.bounds
    crs = r_ds.crs
poly = Polygon([(l,b),(r,b),(r,t),(l,t)])
fiona_config = {
    'schema':{
        'geometry':'Polygon',
        'properties': {('name','str')}
    },
    'driver':'GPKG',
    'crs':crs,
}
with fiona.open(out_fn, 'w', **fiona_config) as ds:
    rec = {
        'properties': {'name':raster_fn},
        'geometry': mapping(poly)
    }
    ds.write(rec)
with rasterio.open(out_fn, "r+", raster_table="new_table", append_subdataset=True) as dst_dataset:
    dst_dataset.write(array)

The GPKG is created and valid. Just fails on the rasterio.open call.


Re: Writing raster to GeoPackage layer

Sean Gillies
 

Ryan,

I don't use geopackage rasters and haven't tested the code I'm about to suggest, but I think this is worth a try:

with rasterio.open("testing.gpkg", "r+", raster_table="new_table", append_subdataset=True) as dst_dataset:
    dst_dataset.write(array)

The format driver has a big set of options and it's possible that the ones above are not entirely sufficient for your purposes. See

On Sat, Jun 29, 2019 at 1:11 PM Sean Gillies via Groups.Io <sean.gillies=gmail.com@groups.io> wrote:
It may be that rasterio doesn't support multilayer raster datasets. I'll make a note to open an issue in the tracker about this. 

On Sat, Jun 29, 2019, 10:19 AM Ryan <code@...> wrote:
Installed rasterio through conda:

conda create -n gpkg-test -c conda-forge python=3.7 rasterio

--
Sean Gillies


Re: build_overviews

vincent.sarago@...
 

Hi James, 
The level of overviews is usually linked to the size of the raster and its internal tiles (block).
here is an example in rio-cogeo: https://github.com/cogeotiff/rio-cogeo/blob/master/rio_cogeo/utils.py#L95-L118

Let say you have a 10 000 x 10 000 px raster with 256x256 block/tile size (to allow fast reading). The overview level is then 6 which results in `[2, 4, 8, 16, 32, 64]` levels (https://github.com/cogeotiff/rio-cogeo/blob/master/rio_cogeo/cogeo.py#L223).


build_overviews

James David Smith
 

Hello,

Would anyone happen to know the answer to this question?

https://stackoverflow.com/questions/57005499/rasterio-version-of-arcgis-raster-pyramids

Essentially I am trying to copy the ArcGIS build overviews dialog box.
But I am unsure how many overviews to use or at what levels they
should be.

I guess this is more of a general GDAL question that rasterio specific.

Thanks, James


Re: Combining rasters using rasterio

James David Smith
 

Oh right. Well that sounds promising then. So for a set of single band
geotiffs, would I be doing something like this .... ?

file_one = rasterio.open('raster_one.tif', 'r')
file_two = rasterio.open('raster_two.tif', 'r')
file_three = rasterio.open('raster_three.tif', 'r')

new_raster = rasterio.merge.merge([file_one, file_two, file_three], indexes=1)

Thanks again, James

On Fri, 5 Jul 2019 at 17:33, <vincent.sarago@...> wrote:

It seems the docs it wrong, looking at the code it take the max extent of all the dataset https://github.com/mapbox/rasterio/blob/b9f34ee559039239c7c0c97bd911b466701a39cd/rasterio/merge.py#L77-L92


Re: Combining rasters using rasterio

vincent.sarago@...
 

It seems the docs it wrong, looking at the code it take the max extent of all the dataset https://github.com/mapbox/rasterio/blob/b9f34ee559039239c7c0c97bd911b466701a39cd/rasterio/merge.py#L77-L92


Re: Combining rasters using rasterio

James David Smith
 

Hi Vincent,

Thanks for the reply. I looked at merge (
https://rasterio.readthedocs.io/en/stable/api/rasterio.merge.html ),
but was a bit confused about what to do for the bounds. The
documentation says:

"Geospatial bounds and resolution of a new output file in the units of
the input file coordinate reference system may be provided and are
otherwise taken from the first input file."

So in my example, the bounds would be set from raster_one . What I
actually need to do is set the bounds as the min_x of all 3 rasters,
the min_y of all 3 rasters, the max_x of all 3 rasters, and the max_y
of all 3 rasters.

I couldn't figure out how to do that.

Thanks, James

On Fri, 5 Jul 2019 at 17:06, <vincent.sarago@...> wrote:

On Wed, Jul 3, 2019 at 06:52 AM, James David Smith wrote:

g. If anyone can point me in the right direction or provide some examples I'd be very grateful. A previous colleague completed this work using ArcPY, and I'm doing my best to re-write the code in open-source. Thanks!

Hi James,
Rasterio has a cli command called `merge` (`rio merge`) which is the equivalent of `gdal_merge.py` with rasterio. If you are looking for a way to do it within python you can check https://github.com/mapbox/rasterio/blob/master/rasterio/rio/merge.py

Vincent


Re: Combining rasters using rasterio

vincent.sarago@...
 

On Wed, Jul 3, 2019 at 06:52 AM, James David Smith wrote:
g. If anyone can point me in the right direction or provide some examples I'd be very grateful. A previous colleague completed this work using ArcPY, and I'm doing my best to re-write the code in open-source. Thanks!
Hi James,
Rasterio has a cli command called `merge` (`rio merge`) which is the equivalent of `gdal_merge.py` with rasterio. If you are looking for a way to do it within python you can check https://github.com/mapbox/rasterio/blob/master/rasterio/rio/merge.py

Vincent 


Combining rasters using rasterio

James David Smith
 

Hello,

Thanks for the help last week with the OVR files. I think I've solved that now.

My next question is about combining rasters of shared space, but different extents. Each raster has one band. I want to combine them into one raster, keeping the values of the highest cell value where they overlap. Hopefully the attached image explains better than I can in writing. If anyone can point me in the right direction or provide some examples I'd be very grateful. A previous colleague completed this work using ArcPY, and I'm doing my best to re-write the code in open-source. Thanks!

Capture.PNG


Re: Understanding `+wktext` more

Sean Gillies
 

Hi,

On Mon, Jul 1, 2019 at 9:54 AM David Hoese <dhoese@...> wrote:
I posted this originally on github, but have moved it here at Sean's request. Additional details of this question/issue are discussed in the `pyproj` repository [here](https://github.com/pyproj4/pyproj/issues/357).
 
I discovered that `nsper` and `ob_tran` projections required the `+wktext` parameter to be accepted by rasterio's `CRS` objects. It wasn't a big deal to add this until recently where I started using pyproj's `CRS` objects as my library's container for CRS information. Using pyproj 2.0+ and its `CRS` objects, the `+wktext` of my PROJ dictionary is being ignored/removed. This means that doing `rasterio.crs.CRS.from_dict(CRS.to_dict())` fails because rasterio needs the `+wktext` parameter.
 
I'm not sure if this should be considered a bug in pyproj or in rasterio or not a bug at all. I found the related issues in this repository that mentioned `+wktext` but I had trouble understanding what it was for and why it was needed. Any help and information would be appreciated.
 
## Examples
 
from pyproj import CRS                                                                                                                                                            
from rasterio.crs import CRS as rioCRS                                                                                                                                            
 
proj_dict = {'proj': 'ob_tran', 'o_proj': 'eqc', 'o_lat_p': 30.0, 'o_lon_p': 10.0, 'lon_0': -10.0, 'a': 6371000.0, 'b': 6371000.0, 'wktext': True}
 
CRS.from_dict(proj_dict).to_dict()                                                                                                                                                
# {'proj': 'ob_tran', 'o_proj': 'eqc', 'o_lat_p': 30.0, 'o_lon_p': 10.0, 'lon_0': -10.0, 'a': 6371000.0, 'b': 6371000.0, 'type': 'crs'}
 
rioCRS.from_dict(proj_dict)                                                                                                                                                       
# CRS.from_wkt('PROJCS["unnamed",GEOGCS["unnamed ellipse",DATUM["unknown",SPHEROID["unnamed",6371000,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["custom_proj4"],EXTENSION["PROJ4","+proj=ob_tran +o_proj=eqc +o_lat_p=30.0 +o_lon_p=10.0 +lon_0=-10.0 +a=6371000.0 +b=6371000.0 +wktext=True"]]')
 
rioCRS.from_dict(CRS.from_dict(proj_dict).to_dict()) 
# CRSError: The PROJ4 dict could not be understood. OGR Error code 5
The latest version of pyproj uses PROJ version 6, which is based on WKT (version 2) and not on the old PROJ.4 key/value representation of projections. I think +wktext may be gone from PROJ 6 as its no longer needed to carry extra WKT information into the PROJ.4 realm. If you really want to use the pyproj.CRS class like this, you'll need to seek help from that project.

Furthermore dicts are no longer the highest fidelity interface between pyproj and rasterio. To make a rasterio CRS object from a pyproj CRS object, export WKT (version 1) from the latter [1] and send that to rasterio.crs.CRS.from_wkt(). I haven't tested this, but am confident that it's the best approach.


--
Sean Gillies

681 - 700 of 956