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

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

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

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

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

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

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

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

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

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

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

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

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: Unexpected result using rasterio windowed writing - Access window out of range in RasterIO()

umbertofilippo@...
 

Replying to myself to show my last update on the situation...

I have made some interesting improvements:

I am now able to write the full extent of raster `A` correctly over raster `B`, having the same pixel values in the same locations for both rasters.

The only thing left before having a complete solution is that the output raster has some (apparently regularly spaced) stripes, sign that there might be some adjustments to do in the window coordinates translation.

Here is a picture of the updated situation:



zoomed to see the stripes:



Finally, the updated code:

    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(
                    row=window.row_off,
                    col=window.col_off)
                right, bottom = src.xy(
                    row=window.row_off + window.height,
                    col=window.col_off + window.width)
                new_window = from_bounds(
                    left=left,
                    bottom=bottom,
                    right=right,
                    top=top,
                    transform=mosaic_raster.transform,
                    height=window.height,
                    width=window.width)
                mosaic_raster.write(r, 1, window=new_window)

Trying to make custom WKT for a non-Earth projection to spatially reference images taken with microscope

Ryan Avery
 

I'm working with non spatially-referenced images produced from a laser ablation system. Along with the image, the software produces an "align" file of the same name - an xml file that contains the XY center point of each image (coordinates are in microns) and the XY size (in microns) along with any rotation.

I would like to display the images (we typically generate about 1000 images / day) in their correct relative positions and sizes.

The laser has a stage with a cartesian positive XY coordinate system (in microns), and the origin in the SW corner.

So far I have gotten to the step of generating an affine transform for each scanned image based on it's size and center coordinate in microns and rotation metadata. What I don't have is a definition of a custom CRS in WKT. I don't have any experience with the WKT format but it seems like the way to represent non-Earth CRSs. Once I have this I think I reproject each image and save as a geotiff for later viewing in QGIS.

Are there any resources for defining a custom CRS in WKT or has anyone done something like this before? I know the units are in microns, the center point of my CRS, and the height and width of the coordinate reference system (and that it is planar). Any help is much appreciated, my progress so far is located here: https://github.com/rbavery/custom_spatial_reference/blob/master/spatial_ref_imgs.ipynb

Best,
Ryan

Using Rasterio with GDAL 2.4.x

Ratcliff, Christina (A&F, Waite Campus)
 

Hi,

 

I’m trying to determine if it is possible to install Rasterio against GDAL 2.4.x on windows using the Gohlke wheels.

 

When installing rasterio versions 1.0.25 and above, it detects that the GDAL requirement is not met and attempts to build and install GDAL 3.0.4 and fails (log attached).
The environment I’m installing into is QGIS3 has GDAL 2.4.1 which shouldn’t be upgraded.

 

Is it possible to install the current release of Rasterio against GDAL 2.4.x ? The Rasterio readme suggests that this is supported.

 

On the assumption that this behavior is in error, I have been looking at source files and have a couple of queries.

 

1. In the .travis.yml are the following lines

- GDALINST=$HOME/gdalinstall
- GDALBUILD=$HOME/gdalbuild
- PROJINST=$HOME/gdalinstall
- PROJBUILD=$HOME/projbuild

 

      Should the 3rd line read $home/projinstall?

 

2. I have when installing using the wheel file, that GDAL is being matched to ~=3.0.1. (see attached log)

Collecting gdal~=3.0.1

  Using cached GDAL-3.0.4.tar.gz (577 kB)


    Where is this specified in the source files? 
If Rasterio indeed supports GDAL 2.4.1 should this in fact be >= 2.x.x ?


Many thanks in advance,

Christina

Re: Trying to make custom WKT for a non-Earth projection to spatially reference images taken with microscope

Sean Gillies
 

Hi Ryan,

On Sun, Feb 9, 2020 at 1:55 PM Ryan Avery <ravery@...> wrote:

I'm working with non spatially-referenced images produced from a laser ablation system. Along with the image, the software produces an "align" file of the same name - an xml file that contains the XY center point of each image (coordinates are in microns) and the XY size (in microns) along with any rotation.

I would like to display the images (we typically generate about 1000 images / day) in their correct relative positions and sizes.

The laser has a stage with a cartesian positive XY coordinate system (in microns), and the origin in the SW corner.

So far I have gotten to the step of generating an affine transform for each scanned image based on it's size and center coordinate in microns and rotation metadata. What I don't have is a definition of a custom CRS in WKT. I don't have any experience with the WKT format but it seems like the way to represent non-Earth CRSs. Once I have this I think I reproject each image and save as a geotiff for later viewing in QGIS.

Are there any resources for defining a custom CRS in WKT or has anyone done something like this before? I know the units are in microns, the center point of my CRS, and the height and width of the coordinate reference system (and that it is planar). Any help is much appreciated, my progress so far is located here: https://github.com/rbavery/custom_spatial_reference/blob/master/spatial_ref_imgs.ipynb

Best,
Ryan


This sounds like a job for WKT2 EngineeringCRS.


I do not know how to use it with rasterio, but I would love to learn how.

--
Sean Gillies