Date   

Unexpected result using rasterio windowed writing - Access window out of range in RasterIO()

umbertofilippo@...
 

Hello everybody,

this is my first post in this group. It is copied directly from my question on gis.stackexchange.com in the hope I can receive a feedback to solve (or understand) my problem.
The question is this one (https://gis.stackexchange.com/questions/349183/unexpected-result-using-rasterio-windowed-writing-access-window-out-of-range-i).

I am trying to write (copy and paste) a raster into another raster using raserio's blocks and windows.

I have two rasters (`A` and `B`):

`A` has a greater extent than `B` and it completely contains B.
The two rasters share the same cell size (roughly 30m) and CRS (EPSG:4326).

The lower and right boundaries of `A` and `B` overlap like below:



I am having hard times understanding the problem of my code. I think the issue may be in the order I am iterating through each block when using `block_windows()` method, but I might be wrong.

Looking at the picture in the rasterio's documentation and printing the block indexes at each iteration, it seems the blocks should write east to west (let's say "horizontally"), and once the upper part is all written, the process should jump to the next lower row.



However, when I run my code, it seems the writing process happens "vertically" (north to south), and then jump to the next column). For this problem, when I reach the end of the first column, the code tries to write into a block that is outside `A` extent and I receive the error:

RasterioIOError('Read or write failed. /home/umberto/test/output/mosaic2.tif: Access window out of range in RasterIO().  Requested (42379,56159) of size 256x256 on raster of 85196x56160.')

Here is what I see in QGIS:



The black rectangle represents raster `A` where I am trying to paste `B` block by block. The whitish part in the middle of the raster represents the 256X256 blocks in the first column of `B` written to `A`.

Unfortunately, the writing process cannot jump to the next column as the the indexes of rows makes the code trying to access a part of the raster which is below the raster extent itself.

How am I doing it

The idea is to write each window of raster `B` into raster `A`. However, as the window object column and row offsets are different when considering the two rasters, I am translating for each block the offsets from `B` to `A` (window and new_window in the below code).

Here is my code (simplified but still relevant):

    with rasterio.open(inputfile, 'r') as src:
        with rasterio.open(
            os.path.join(dest_dir, 'mosaic2.tif'),
            'w',
            driver='GTiff',
            compress="DEFLATE",
            height=56160,
            width=85196,
            count=1,
            dtype=np.float64,
            crs='+proj=latlong',
            transform=out_transform,
        ) as mosaic_raster:
            for ji, window in src.block_windows(1):
                r = src.read(1, window=window)
                left, top = src.xy(
                    window.col_off, window.row_off, offset='ul')
                right, bottom = src.xy(
                    window.col_off + window.width, window.row_off + window.height, offset='ul')
                new_window = from_bounds(
                    left, bottom, right, top, mosaic_raster.transform, window.height, window.width)
                mosaic_raster.write(r, 1, window=Window(int(new_window.col_off), int(
                    new_window.row_off), int(new_window.width), int(new_window.height)))


Re: Memory error, rasterio merge

Sean Gillies
 

Hi Simon,

On Wed, Jan 29, 2020 at 7:29 AM Simon via Groups.Io <twinbirds=protonmail.com@groups.io> wrote:
Hello everyone,

I tried to merge/stitch my 1000 geotif images, each of around 4gb using rasterio.

files = glob.glob('*.tif')
src_files = []
for fi in files:
    src = rasterio.open(fi)
    src_files.append(src)
mosaic, out_trans = merge(src_files)
out_meta = src.meta.copy()
out_meta.update({"driver": "GTiff",
                  "height": mosaic.shape[1],
                  "width": mosaic.shape[2],
                  "transform": out_trans,
                  "crs": src.crs})

with rasterio.open(fo, "w", **out_meta) as dest:
     dest.write(mosaic)

Unfortunately, I could not proceed with the Memory error.
I could merge/stitch them using gdal translate, though it was very slow. So, physical memory should not be problem. To speed up, I  tried rasterio, but did not work. Any help to make it possible would be appreciated.

Thank you.
Simon

Rasterio's merge function is not well optimized with respect to memory. An output array is allocated such that it covers all the inputs. If your inputs don't overlap, this could be 4TB or more, which probably exceeds the available memory on your computer.

Best wishes,

--
Sean Gillies


Memory error, rasterio merge

Simon
 

Hello everyone,

I tried to merge/stitch my 1000 geotif images, each of around 4gb using rasterio.

files = glob.glob('*.tif')
src_files = []
for fi in files:
    src = rasterio.open(fi)
    src_files.append(src)
mosaic, out_trans = merge(src_files)
out_meta = src.meta.copy()
out_meta.update({"driver": "GTiff",
                  "height": mosaic.shape[1],
                  "width": mosaic.shape[2],
                  "transform": out_trans,
                  "crs": src.crs})

with rasterio.open(fo, "w", **out_meta) as dest:
     dest.write(mosaic)

Unfortunately, I could not proceed with the Memory error.
I could merge/stitch them using gdal translate, though it was very slow. So, physical memory should not be problem. To speed up, I  tried rasterio, but did not work. Any help to make it possible would be appreciated.

Thank you.
Simon


Re: Silencing NotGeoreferencedWarning

Yann-Sebastien Tremblay-Johnston
 

Just chiming in to say I've gotten busy and haven't had a chance to polish https://github.com/mapbox/rasterio/pull/1845 off yet


On Thu, Jan 23, 2020 at 9:56 AM Guillaume Lostis <g.lostis@...> wrote:
Hi Sean,

Thanks for your answer.
I understand why you would want to keep this as a warning so that a user can act upon it if needed.

If I understand correctly, rasterio currently lacks the ability to georeference an image properly through its RPC/GCP, so you will keep the warning as is for now?
I'd be happy to contribute what I can, at my level, to help move forward on these topics.
On RPCs, @underchemist has a WIP open here: https://github.com/mapbox/rasterio/pull/1845
On GCPs, Vincent contributed a method to estimate a transform matrix here: https://github.com/mapbox/rasterio/pull/1749

On Thu, 23 Jan 2020 at 18:10, Sean Gillies via Groups.Io <sean=mapbox.com@groups.io> wrote:
Hi Guillaume,

Rasterio logs, mostly at DEBUG level, so that you can run programs and see what the state was before a failure. You don't see them unless you turn the warning level down to DEBUG. A program can't easily react to log messages, they go out to the environment and are inspected (or discarded) after the program terminates.

Rasterio warns in a few places instead of using logging because at these places a user may want their program to react to the warning and use different logic, retry, add missing georeferencing, fix the alpha band shadowing, or whatever. An exception would be too strong in the georeferencing case, I think. But we can, with a couple lines of code, turn any NotGeoreferencedWarning into an exception that can be handled appropriately.

Certainly rasterio is lagging in the RPC/GCP realm and we should think about not warning after we add the missing capabilities to the package.

On Thu, Jan 23, 2020 at 7:49 AM Guillaume Lostis <g.lostis@...> wrote:
Hi Vincent,

Option 1 with RPCs and GCPs would be perfect for my use case because all of the datasets I open have one of them.
But it would still leave warnings when you write a png or jpg with no georeference as you pointed out.

I have to admit I didn't know rasterio had logs... What's the rationale behind deciding what goes to logs, and what goes to stderr/warnings?

On Thu, 23 Jan 2020 at 15:40, <vincent.sarago@...> wrote:
Hi Guillaume,
thanks for raising this, it happens to me often too.

I'm thinking about multiple solutions:
- check if the file has other georeferenced information (like rpc or gpc) internally before raising this warning. 
- instead of raising a warning maybe a logging.info would be sufficent ? 

I also encounter this when saving file like png or jpg that don't have georeference which appears to me unnecessary. 

--
Sean Gillies


Re: Silencing NotGeoreferencedWarning

Guillaume Lostis
 

Hi Sean,

Thanks for your answer.
I understand why you would want to keep this as a warning so that a user can act upon it if needed.

If I understand correctly, rasterio currently lacks the ability to georeference an image properly through its RPC/GCP, so you will keep the warning as is for now?
I'd be happy to contribute what I can, at my level, to help move forward on these topics.
On RPCs, @underchemist has a WIP open here: https://github.com/mapbox/rasterio/pull/1845
On GCPs, Vincent contributed a method to estimate a transform matrix here: https://github.com/mapbox/rasterio/pull/1749

On Thu, 23 Jan 2020 at 18:10, Sean Gillies via Groups.Io <sean=mapbox.com@groups.io> wrote:
Hi Guillaume,

Rasterio logs, mostly at DEBUG level, so that you can run programs and see what the state was before a failure. You don't see them unless you turn the warning level down to DEBUG. A program can't easily react to log messages, they go out to the environment and are inspected (or discarded) after the program terminates.

Rasterio warns in a few places instead of using logging because at these places a user may want their program to react to the warning and use different logic, retry, add missing georeferencing, fix the alpha band shadowing, or whatever. An exception would be too strong in the georeferencing case, I think. But we can, with a couple lines of code, turn any NotGeoreferencedWarning into an exception that can be handled appropriately.

Certainly rasterio is lagging in the RPC/GCP realm and we should think about not warning after we add the missing capabilities to the package.

On Thu, Jan 23, 2020 at 7:49 AM Guillaume Lostis <g.lostis@...> wrote:
Hi Vincent,

Option 1 with RPCs and GCPs would be perfect for my use case because all of the datasets I open have one of them.
But it would still leave warnings when you write a png or jpg with no georeference as you pointed out.

I have to admit I didn't know rasterio had logs... What's the rationale behind deciding what goes to logs, and what goes to stderr/warnings?

On Thu, 23 Jan 2020 at 15:40, <vincent.sarago@...> wrote:
Hi Guillaume,
thanks for raising this, it happens to me often too.

I'm thinking about multiple solutions:
- check if the file has other georeferenced information (like rpc or gpc) internally before raising this warning. 
- instead of raising a warning maybe a logging.info would be sufficent ? 

I also encounter this when saving file like png or jpg that don't have georeference which appears to me unnecessary. 

--
Sean Gillies


Re: Silencing NotGeoreferencedWarning

Sean Gillies
 

Hi Guillaume,

Rasterio logs, mostly at DEBUG level, so that you can run programs and see what the state was before a failure. You don't see them unless you turn the warning level down to DEBUG. A program can't easily react to log messages, they go out to the environment and are inspected (or discarded) after the program terminates.

Rasterio warns in a few places instead of using logging because at these places a user may want their program to react to the warning and use different logic, retry, add missing georeferencing, fix the alpha band shadowing, or whatever. An exception would be too strong in the georeferencing case, I think. But we can, with a couple lines of code, turn any NotGeoreferencedWarning into an exception that can be handled appropriately.

Certainly rasterio is lagging in the RPC/GCP realm and we should think about not warning after we add the missing capabilities to the package.


On Thu, Jan 23, 2020 at 7:49 AM Guillaume Lostis <g.lostis@...> wrote:
Hi Vincent,

Option 1 with RPCs and GCPs would be perfect for my use case because all of the datasets I open have one of them.
But it would still leave warnings when you write a png or jpg with no georeference as you pointed out.

I have to admit I didn't know rasterio had logs... What's the rationale behind deciding what goes to logs, and what goes to stderr/warnings?

On Thu, 23 Jan 2020 at 15:40, <vincent.sarago@...> wrote:
Hi Guillaume,
thanks for raising this, it happens to me often too.

I'm thinking about multiple solutions:
- check if the file has other georeferenced information (like rpc or gpc) internally before raising this warning. 
- instead of raising a warning maybe a logging.info would be sufficent ? 

I also encounter this when saving file like png or jpg that don't have georeference which appears to me unnecessary. 

--
Sean Gillies


Re: Silencing NotGeoreferencedWarning

Guillaume Lostis
 

Hi Vincent,

Option 1 with RPCs and GCPs would be perfect for my use case because all of the datasets I open have one of them.
But it would still leave warnings when you write a png or jpg with no georeference as you pointed out.

I have to admit I didn't know rasterio had logs... What's the rationale behind deciding what goes to logs, and what goes to stderr/warnings?

On Thu, 23 Jan 2020 at 15:40, <vincent.sarago@...> wrote:
Hi Guillaume,
thanks for raising this, it happens to me often too.

I'm thinking about multiple solutions:
- check if the file has other georeferenced information (like rpc or gpc) internally before raising this warning. 
- instead of raising a warning maybe a logging.info would be sufficent ? 

I also encounter this when saving file like png or jpg that don't have georeference which appears to me unnecessary. 


Re: Silencing NotGeoreferencedWarning

vincent.sarago@...
 

Hi Guillaume,
thanks for raising this, it happens to me often too.

I'm thinking about multiple solutions:
- check if the file has other georeferenced information (like rpc or gpc) internally before raising this warning. 
- instead of raising a warning maybe a logging.info would be sufficent ? 

I also encounter this when saving file like png or jpg that don't have georeference which appears to me unnecessary. 


Silencing NotGeoreferencedWarning

Guillaume Lostis
 

Hi,

I'm frequently working with GeoTIFF datasets that don't have a geotransform because they are georeferenced through other means (Sentinel-1 SLC images through the satellite state vectors, or non-orthorectified optical imagery through RPCs). Whenever I open such a dataset with rasterio, I get a warning:

NotGeoreferencedWarning: Dataset has no geotransform set. The identity matrix may be returned.

This warning is not useful to me as it is raised on every file I open, so I was looking for a way to silence it systematically. I've found that by doing:

import warnings
import rasterio
warnings.filterwarnings("ignore", category=rasterio.errors.NotGeoreferencedWarning)

it silences most warnings, but sometimes some of them get through the filter in for loops for example, or in the context of pytest. And it also has the disadvantage of having to show up everywhere in my code.

Is there a more radical way to always silence this warning?

Thanks!

Guillaume Lostis


Re: returning subset of raster data via API

Sean Gillies
 

Hi Tom,

On Fri, Jan 17, 2020 at 8:35 AM Tom Kralidis <tom.kralidis@...> wrote:

Hi all: as part of working on OGC API - Coverages support in pygeoapi, I'm writing the equivalent of what WCS GetCoverage would do; that is, clip by bbox and subset by bands if requested.

An initial implementation can be found in this branch.

At the moment, the data is returned as JSON-ified ndarrays. I would like to return the data in its native format. What would be a general approach here? Would one use a MemoryFile, write to it, and then read it back to the caller?

Thanks

When you say that you want to return the data in its native format, do you mean a stream of bytes that contain, for example, a GeoTIFF? The MemoryFile class would serve you well in that case if you wanted to avoid writing to disk, or one of Python's temporary file objects would do just as well if you didn't mind or wanted a temp file on disk. You would be limited to single-file formats if you used MemoryFile, our abstraction doesn't cover multiple files. It seems to me that there would be a related problem for multiple files on the serialization end, whether you used MemoryFile or not: how do you put multiple files in a single stream of bytes?

--
Sean Gillies


Re: Can Rasterio generate a geoTIFF from a matrix?

Luke
 

Is your matrix of values a numpy array? You can easily create a GeoTIFF in rasterio if you have an array of data and the basic metadata (georeferencing, dimensions, band count, data type).

https://rasterio.readthedocs.io/en/latest/quickstart.html#creating-data
https://rasterio.readthedocs.io/en/latest/topics/profiles.html
https://rasterio.readthedocs.io/en/latest/topics/writing.html


returning subset of raster data via API

 

Hi all: as part of working on OGC API - Coverages support in pygeoapi, I'm writing the equivalent of what WCS GetCoverage would do; that is, clip by bbox and subset by bands if requested.

An initial implementation can be found in this branch.

At the moment, the data is returned as JSON-ified ndarrays. I would like to return the data in its native format. What would be a general approach here? Would one use a MemoryFile, write to it, and then read it back to the caller?

Thanks

..Tom


Can Rasterio generate a geoTIFF from a matrix?

Peter Oelslager <poelslager@...>
 

Hi,

I would like to create geoTIFF overlays for RF propagation.  Does Rasterio allow users to generate geoTIFF files from a matrix of values?  If that capability does not currently exist can I try contributing to the project?  I've been able to bit-bang out a couple geoTIFF files but its all sort of hacked together at the moment.  If this kind of feature sounds useful I would like to learn how to use git and contribute code to an actual project.

Thanks.
-Pete


Re: window from_bounds behaviour

Sean Gillies
 

Hi Tom,

Thanks for bringing this up and for joining us here. I've made an issue here: https://github.com/mapbox/rasterio/issues/1857.

Happy New Year!

On Thu, Jan 9, 2020 at 4:53 AM Tom Kralidis <tom.kralidis@...> wrote:

Hi Dion: thanks for the info. Should transform be a required argument (instead of a TypeError traceback)? Happy to file an issue if useful.



--
Sean Gillies


Re: window from_bounds behaviour

 

Hi Dion: thanks for the info. Should transform be a required argument (instead of a TypeError traceback)? Happy to file an issue if useful.


Re: window from_bounds behaviour

Dion Häfner
 

Hey Tom,

rasterio windows are constructed at a specific resolution, so you need to supply a transform to use. It is usually easiest to use the transform of the raster you want to use the window on.

So something like this should work:

import rasterio, rasterio.windows
with rasterio.open('myraster.tif') as src:
... window = rasterio.windows.from_bounds(
... -150, 40, -45, 90, transform=src.transform
... )
... data = src.read(window=window)

Hope that helps.

Dion


On 08/01/2020 23.49, Tom Kralidis via Groups.Io wrote:
Hi all: first off, kudos to the rasterio project, very well done!
Using 1.1.2 I'm getting some weird behaviour when trying to create a window from bounding coordinates.
|>>> from rasterio.windows import from_bounds >>> window = from_bounds(-150, 40, -45, 90) Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/Users/tomkralidis/Dev/pygeoapi/lib/python3.7/site-packages/rasterio/windows.py", line 278, in from_bounds transform, left, top, op=float, precision=precision) File "/Users/tomkralidis/Dev/pygeoapi/lib/python3.7/site-packages/rasterio/transform.py", line 238, in rowcol invtransform = ~transform TypeError: bad operand type for unary ~: 'NoneType' |
Any idea what I am doing wrong here?
Thanks
..Tom
P.S. aside: The pygeoapi <https://pygeoapi.io> project is implementing the OGC API - Coverages specification and will be using rasterio as the default coverage provider.


window from_bounds behaviour

 

Hi all: first off, kudos to the rasterio project, very well done!

Using 1.1.2 I'm getting some weird behaviour when trying to create a window from bounding coordinates.

>>> from rasterio.windows import from_bounds
>>> window = from_bounds(-150, 40, -45, 90)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/tomkralidis/Dev/pygeoapi/lib/python3.7/site-packages/rasterio/windows.py", line 278, in from_bounds
    transform, left, top, op=float, precision=precision)
  File "/Users/tomkralidis/Dev/pygeoapi/lib/python3.7/site-packages/rasterio/transform.py", line 238, in rowcol
    invtransform = ~transform
TypeError: bad operand type for unary ~: 'NoneType'

Any idea what I am doing wrong here?

Thanks

..Tom

P.S. aside: The pygeoapi project is implementing the OGC API - Coverages specification and will be using rasterio as the default coverage provider.


Re: need help installing rasterio with WebP

Howard Butler
 

The reason this doesn't work is Conda/Conda Forge uses an external libtiff build which does not have WebP support enabled. The "WEBP" GDAL driver by itself should work because that is linked differently (but that isn't a TIFF/WebP either). I have initiated an update to Conda Forge's libtiff to add WebP support, but this is likely to take a while for this build to filter back through. See https://github.com/conda-forge/libtiff-feedstock/pull/48 for more details.

Howard


On Thu, Dec 5, 2019 at 1:45 PM <tgertin@...> wrote:
Thank you for looking into this. I am running these commands using the conda-forge channel:

```
conda create --name rasterio_w_webp_test1 python=3.7
 
conda config --add channels conda-forge
 
conda install gdal libgdal
 
conda install rasterio
 
pip install rio-cogeo
```

Then I am using the rio cogeo command that uses WebP:
```
rio cogeo create -p webp input.tif output.tif
```

but I get the following error:
```
File "rasterio/shutil.pyx", line 139, in rasterio.shutil.copy
  File "rasterio/_err.pyx", line 205, in rasterio._err.exc_wrap_pointer
rasterio._err.CPLE_AppDefinedError: Cannot create TIFF file due to missing codec for WEBP.
```

so I think something didn't install correctly. I don't need to use conda, but it seems like a good option because it uses virtual environments.



Re: Reading NetCDF file as inMemoryFile

Even Rouault
 

On vendredi 20 décembre 2019 09:45:45 CET vincent.sarago@... wrote:
Yes it is set to yes,
You were responding to "userfaultfd support: yes" ?

Hum, then I'm not sure. Are you running in a container ? Maybe there are
some restrictions by default. Dunno. But you should see GDAL error messages
if userfaulfd system calls fail at runtime. There are quite a lot of them in
https://github.com/OSGeo/gdal/blob/master/gdal/port/cpl_userfaultfd.cpp

And fallback to the magic solution of open source projects (the reason why
we all use open source, right ;-) ?): take your favorite debugger and break at
https://github.com/OSGeo/gdal/blob/master/gdal/frmts/netcdf/netcdfdataset.cpp#L7274
and follow what happens then...
It should normally go to the call to nc_open_mem()

--
Spatialys - Geospatial professional services
http://www.spatialys.com


Re: Reading NetCDF file as inMemoryFile

vincent.sarago@...
 

Yes it is set to yes,

here the full log https://gist.github.com/vincentsarago/36473e6322336e84cf928ef445db64cc