Date   

A question about the rasterio precision & python floats (Pixel width & height)

thomaswanderer@...
 

When opening my raster with QGIS, I can see the Pixel Size of my raster as:

0.0008333333299974727193,-0.000833333330020215703

When I read he same raster with rasterio, I get pixel size values with a lower precision and I guess this causes some problems when calculating new affines for another extent

0.0008333333299974727 0.0008333333300202157

I have tries to pass the decimal_precision=X when opening my raster with rasterio, but this has no effect. I guess this is a general float problem on my system. My suggestion therefore would be to use python Decimals or any other type representing real floats.


Re: Create raster aligned with one created by gdal_translate

thomaswanderer@...
 
Edited

I had a very similar approach to calculate the new affine & width/height manually as you. But I have now replaced it with just rasterio functions and when comparing the results of the new affine, they were identical :-)

I just wonder, why you add 0.5 when calculating the height & width? My code looked like this (windows is the extent I am matching my raster to, profile is the target raster profile, copied from my source raster and then modified)
if window is not None and window.is_polygon:
                    # Make sure we operate with the same coordinates
                    window = window.transform(profile['crs'].to_epsg())

                    # Get the coordinates from raster and window bounds
                    s_xmin, s_ymin, s_xmax, s_ymax = raster_bounds.bounds
                    w_xmin, w_ymin, w_xmax, w_ymax = window.bounds

                    # Get the source pixel width and height
                    x_pixel = profile['transform'].a
                    y_pixel = profile['transform'].e

                    # Calculate the shift vector from source bounds to target bounds
                    shift_xmin = w_xmin - s_xmin
                    shift_ymin = w_ymin - s_ymin
                    shift_xmax = w_xmax - s_xmax
                    shift_ymax = w_ymax - s_ymax

                    # Normalize the shift vector by the source pixel width and height
                    shift_xmin_norm = ceil(shift_xmin / x_pixel) * x_pixel
                    shift_ymin_norm = floor(shift_ymin / y_pixel) * y_pixel
                    shift_xmax_norm = floor(shift_xmax / x_pixel) * x_pixel
                    shift_ymax_norm = ceil(shift_ymax / y_pixel) * y_pixel

                    # Create the new boundaries of the target raster
                    new_xmin = s_xmin + shift_xmin_norm
                    new_ymin = s_ymin + shift_ymin_norm
                    new_xmax = s_xmax + shift_xmax_norm
                    new_ymax = s_ymax + shift_ymax_norm

                    # Calculate width & height (Pixel) for the target raster
                    new_width = abs(int((new_xmax - new_xmin) / x_pixel))
                    new_height = abs(int((new_ymax - new_ymin) / y_pixel))
                    if new_width == 0 or new_height == 0:
                        raise Exception(f'New raster has invalid measure (width or height is 0 pixels)')

                    # Update profile with windows dimensions
                    profile.update(width=new_width)
                    profile.update(height=new_height)

                    # Create new affine and update profile with it
                    new_affine = rasterio.transform.from_origin(
                        new_xmin,
                        new_ymax,
                        new_width,
                        new_height
                    )
                    profile.update(transform=new_affine)

But i have replaced all this now by:


# Make sure we operate with the same coordinates
# Make sure we operate with the same coordinates window = window.transform(profile['crs'].to_epsg()) # Create a new window and affine from the window extent raster_window = source.window(*window.bounds) raster_window = raster_window.round_offsets(op='ceil').round_lengths(op='floor') # Update target raster with the new affine new_affine = source.window_transform(raster_window) profile.update(transform=new_affine) # Update target raster with new windows shape (Make sure raster stays within extent) c_xmin, c_ymin, c_xmax, c_ymax = source.window_bounds(raster_window) w_xmin, w_ymin, w_xmax, w_ymax = window.bounds x_size, y_size = source.res corrected_width = floor(abs(c_xmin - w_xmax) / x_size) corrected_height = floor(abs(w_ymin - c_ymax) / y_size) profile.update(width=corrected_width) profile.update(height=corrected_height)

Both methods actually created and affine with the same values, that is why I replaced my own code with the rasterio functions.


Re: Window from bounds is shifted by 1 Pixel

thomaswanderer@...
 
Edited

Good to know I am not the only one.
I guess it could be rounding problems, but the extent I use for creating the window aligns perfectly with the input rater (as seen above). But the output window get's shifted and only if my extent overlaps on the western side of the original source raster. All other windows (created from the intersection extents) are just fine for both input and output! Very strange. Ah by the way, I use rasterio 1.13.


Re: Window from bounds is shifted by 1 Pixel

Denis Rykov
 

Hi! I had faced the similar shifting issue yesterday: https://github.com/mapbox/rasterio/issues/1932


Window from bounds is shifted by 1 Pixel

thomaswanderer@...
 

Hi. this is my first post. I am new to raster.io and really like it so far. Thumbs up for all the work :-)

But I have already a first problem: I created a function which creates a new raster for a given geographical extent and copies the values of the source raster into it. The new raster has the same characteristics as the source raster. The geographical extents can now
  • include the source raster (within)
  • lay within the source raster bounds (contained)
  • or partially overlap with the source raster bounds


For all of these cases, I create a new affine for the new raster and then fill it with the values of the source raster (or any other optionally specified value).
I still have a problem when copying over the values from the source to the target raster.
In one case the two rasters (source -> copy) are not aligned but shifted by 1 pixel. (See yellow extent in above images - All other extents create perfectly aligned new rasters)

My strategy for copying data from the source to the target is the following:

  1. I have calculated the exact pixel aligned bounds for the target raster
  2. I have the exact pixel aligned bounds from the source raster
  3. I create an intersection polygon from both
    intersection = INTERSECT(source.bounds, target.bounds)
  4. I use the intersection to create a read-window for the source
    with rasterio.open('source.tif', 'r') as source:
        read_window = source.window(*intersection.bounds))
  5. I use the intersection to create a write-window for the target
    with rasterio.open('target.tif', 'w', **profile) as target:
        write_window = target.window(*intersection.bounds))
  6. Then I read from the source and write to the target
    target.write(source.read(window=read_window), window=write_window)
For all rasters which overlap the western side of the source, I receive values which are shifted by 1 pixel west wards.



My guess is that the write_window = target.window(*intersection.bounds)) is the part where the error happens. While the read_window still reads the correct data, the write window places the copied data shifted by 1 pixel. I use the same intersection extent for the creation of both windows, and the intersection is also pixel aligned as visible in the above image. So I wonder what I do wrong. Is this maybe a rounding or precision error?


Create raster aligned with one created by gdal_translate

Denis Rykov
 

Hi everyone!

I tried to convert gdal_translate command to rasterio and faced some problem I described here. The original command is:

gdal_translate -oo OVERVIEW_LEVEL=1 -projwin 10816494.68930465 3653956.75283452 10819454.39626612 3650628.7941839 /vsicurl/https://raw.githubusercontent.com/satellogic/telluric/master/tests/data/google.xml gdalized.tif

Eventually I came up with the following solution that produces raster aligned with one produced by gdal_translate:

import rasterio
from rasterio.windows import Window
from math import floor
with rasterio.open(
    "/vsicurl/https://raw.githubusercontent.com/satellogic/telluric/master/tests/data/google.xml",
    OVERVIEW_LEVEL=1,
) as src:
    bounds = (10816494.68930465, 3650628.7941839, 10819454.39626612, 3653956.75283452)
    xoff = floor((bounds[0] - src.transform.c) / src.transform.a) - 1
    yoff = floor((bounds[3] - src.transform.f) / src.transform.e) - 1
    xsize = int((bounds[2] - bounds[0]) / src.transform.a + 0.5)
    ysize = int((bounds[1] - bounds[3]) / src.transform.e + 0.5)
    out_window = Window(xoff, yoff, xsize, ysize)
    height = int(out_window.height)
    width = int(out_window.width)
    out_kwargs = src.profile
    out_kwargs.update(
        {
            "driver": "GTiff",
            "height": height,
            "width": width,
            "transform": src.window_transform(out_window),
        }
    )
    with rasterio.open("aligned.tif", "w", **out_kwargs) as out:
        out.write(src.read(window=out_window, out_shape=(src.count, height, width),))

It works fine but I still have couple of questions:

  • Why I have to subtract 1 to get right values of xoff and yoff?
  • How can I adjust spatial extent of raster to match extent of raster produced by gdal_translate? Now there is 1 px difference between them (but it is definitely not shift because all overlapped pixels have the same value).

Any advice would be appreciated.


Re: Create sparse raster with Rasterio. Is this possible?

Brendan Ward
 

Ben - I'm not sure I follow your intent here.

Is the idea to make the smallest possible geotiff after rasterizing all shapes in test.shp to pixels, where SPARSE_OK omits non-rasterized pixels where possible from the geotiff?

In rasterio, the best way to approach that would be to first do a windowed read of your data from raster that overlap the bounding box of all features in test.shp (see geometry_window() https://github.com/mapbox/rasterio/blob/master/rasterio/features.py#L387), then rasterize within that smaller window, then write it out to a new geotiff (note: it will have a different geotransform / width / height than the original).

If your shapes in test.shape are scattered around the raster with large areas of nodata in between, there's not really a way to handle that case aside from windowing out and writing each rasterized shape to a separate geotiff, then maybe making a VRT of the individual geotiffs it in GDAL to mosaic them back together.


Create sparse raster with Rasterio. Is this possible?

BDHudson@...
 

Hi All -

I have benefited greatly from Rasterio - and first want to say thanks.

Second, I am trying to move bash GDAL code into rasterio and python. Is it possible to use rasterio rasterize to create a sparse raster?

In gdal_rasterize setting the SPARSE_OK=TRUE flag does this. Because rasterio rasterize creates a numpy array, (not a geotiff) I am assuming this may not be possible but wanted to ask. I also might imagine a way to use scipy sparse matrices to do something like this. I am doing this to save on memory requirements. 

here is the full call - any advice on how to replicate appreciated. 

gdal_rasterize \
test.shp \
-te $xmin $ymin $xmax $ymax \
-burn 1.0 -tr $res $res -ot Int16 \
-of GTiff -co COMPRESS=ZIP -co BIGTIFF=YES -co SPARSE_OK=TRUE -co NUM_THREADS=ALL_CPUS \
$raster_test

Regards,
Ben


Re: Does Rasterio support tiling of large geotiffs?

Sean Gillies
 

Hi,

On Sun, May 17, 2020 at 6:35 AM <ashnair0007@...> wrote:
I would like to crop my large geotiff into 800 x 800 tiles. While it's possible to do so otherwise, I was wondering if rasterio has an internal functionality to crop large geotiffs to tiles of a specific size.

The rasterio package includes a clip command: https://rasterio.readthedocs.io/en/latest/cli.html#clip. I don't remember how we arrived at calling this "clip" instead of "crop", but it's the same idea. Internally, the clip command uses rasterio's windowed reading feature: https://rasterio.readthedocs.io/en/latest/topics/windowed-rw.html#reading.

-- 
Sean Gillies


Does Rasterio support tiling of large geotiffs?

ashnair0007@...
 

I would like to crop my large geotiff into 800 x 800 tiles. While it's possible to do so otherwise, I was wondering if rasterio has an internal functionality to crop large geotiffs to tiles of a specific size.


Re: rasterio 1.1.4

vincent.sarago@...
 

This is really cool, thanks for the hard work Sean! 


rasterio 1.1.4

Sean Gillies
 

Hi all,

Rasterio 1.1.4 has been released. Here are the changes:


Here are the distributions on PyPI


Please note that there is no wheel for Python 2.7 and OS X. We've lost the ability to build them and I'm not sure why. We may need help from the Travis CI team to get that particular build job unstuck. But as we're moving on from supporting Python 2.7, I'm not super concerned about this right now.

The biggest changes in the wheels on PyPI aren't mentioned in rasterio because they are outside the project scope. See https://github.com/rasterio/rasterio-wheels/blob/master/CHANGES.md#2020-05-06 for details. In a nutshell, we're upgrading NetCDF to 4.6.2, patching GDAL 2.4.4 to fix two different bugs, and are using a patched version of auditwheel to add a rasterio-specific tag, "rasterio", to the SONAME of shared libraries in the linux wheels. In rasterio.libs you will see shared libraries with names like libcurl-rasterio-ea538880.so.4.4.0 instead of libcurl-ea538880.so.4.4.0.

These wheels include 2.4.4. I felt like upgrading to 3.0 or 3.1 while also trying to fix the SONAME collision problem was too much for me, so we'll tackle the GDAL/PROJ upgrades in a future release.

Thanks for your help, everyone!

--
Sean Gillies


Re: Rasterio.features.geometry_window error when geojson has GeometryCollection

filthyhands@...
 

Thanks again Brendan, very helpful. 


Re: Rasterio.features.geometry_window error when geojson has GeometryCollection

Brendan Ward
 

The target release date that will include this fix is 5/1/2020 (https://github.com/mapbox/rasterio/issues/1911)

If you need changes sooner, you have a few options:
1) clone the repo and checkout the maint-1.1 branch, and rebuild the Cython files.  This can be a bit involved depending on your OS / env; you'll need a compiler, GDAL, Cython, etc.  (https://github.com/mapbox/rasterio/blob/master/CONTRIBUTING.rst#development-environment)

2) coerce each GeometryCollection to FeatureCollection

3) use shapely to union each GeometryCollection into a single geometry object before passing into bounds()

4) calculate bounds on each individual geometry within a GeometryCollection, then union the bounds.  Here's how it was done in the Cython file (https://github.com/mapbox/rasterio/pull/1915/files#diff-83bc1acea1f8d435fd6c4798a58fc072R411-R426); you can do the same thing on your end.


Re: MemoryFile loses Profile information

ciaran.evans@...
 

I think the problem is that I'm trying to do the conversion in memory, not on disk. If I've read your latest reply correctly, I'd convert to JPEG2000 on disk and then copy to S3?

Can I do the conversion via rasterio.shutil.copy in memory?


Re: Rasterio.features.geometry_window error when geojson has GeometryCollection

filthyhands@...
 

How can I patch my install of Rasterio to incorporate your changes? 


Re: MemoryFile loses Profile information

Sean Gillies
 

Hi,

On Wed, Apr 29, 2020 at 9:40 AM Guillaume Lostis <g.lostis@...> wrote:

Like Ciaran I'm not very familiar with GDAL's python bindings, but I've tried running a little snippet which I believe does the same thing as what was done with rasterio: it writes a JP2 to a vsimem in-memory file, and then writes the bytes to a file on disk.

from osgeo import gdal

mem, out = "/vsimem/B01.jp2", "new_B01.jp2"
gdal.Translate(mem, "B01.jp2")

# Taken from https://lists.osgeo.org/pipermail/gdal-dev/2016-August/045030.html
f = gdal.VSIFOpenL(mem, "rb")
gdal.VSIFSeekL(f, 0, 2)  # seek to end
size = gdal.VSIFTellL(f)
gdal.VSIFSeekL(f, 0, 0)  # seek to beginning
data = gdal.VSIFReadL(1, size, f)
gdal.Unlink(mem)

with open(out, "wb") as ff:
    ff.write(data)

When I look at the gdalinfo of new_B01.jp2, it has a CRS and a transform, so I am not able to reproduce the behavior we have with rasterio.

(Side note: Interestingly enough, the pixel values of new_B01.jp2 have slightly changed with respect to the original B01.jp2 file, so there is some data loss somewhere in the process. But maybe that is expected of the JP2 format and could be avoided by passing extra arguments to gdal.Translate? I mainly have experience with writing GeoTIFFs, writing JP2s is new to me)

Guillaume

I've only now remembered (not being a regular jp2 user) that JPEG2000 is a create-copy format (see https://gdal.org/drivers/raster/index.html) and as such is not suited for uses cases like

    with rasterio.open("file.jp2", "w", driver="JP2OpenJPEG", ...) as dataset:
        dataset.write(data)

Rasterio tries to abstract over the differences between "create" and "create-copy" formats, but might falter in some cases. We're testing JPEG and PNG in the test suite, but not JPEG2000. If you're constructing a dataset from the bands of multiple other datasets, you should probably be using GeoTIFF as a format, it's well suited for this. And then convert to JPEG2000 before uploading to S3 using rasterio.shutil.copy, which calls on GDAL's GDALCreateCopy and is mostly equivalent to gdal_translate.

--
Sean Gillies


Re: Rasterio.features.geometry_window error when geojson has GeometryCollection

filthyhands@...
 

That's great Brendan, much appreciated. 


Re: MemoryFile loses Profile information

Even Rouault
 

On mercredi 29 avril 2020 08:40:07 CEST Guillaume Lostis wrote:

> Like Ciaran I'm not very familiar with GDAL's python bindings, but I've

> tried running a little snippet which I believe does the same thing as what

> was done with `rasterio`: it writes a JP2 to a `vsimem` in-memory file, and

> then writes the bytes to a file on disk.

>

> ```python

> from osgeo import gdal

>

> mem, out = "/vsimem/B01.jp2", "new_B01.jp2"

> gdal.Translate(mem, "B01.jp2")

>

> # Taken from

> https://lists.osgeo.org/pipermail/gdal-dev/2016-August/045030.html f =

> gdal.VSIFOpenL(mem, "rb")

> gdal.VSIFSeekL(f, 0, 2) # seek to end

> size = gdal.VSIFTellL(f)

> gdal.VSIFSeekL(f, 0, 0) # seek to beginning

> data = gdal.VSIFReadL(1, size, f)

> gdal.Unlink(mem)

>

> with open(out, "wb") as ff:

> ff.write(data)

> ```

>

> When I look at the `gdalinfo` of `new_B01.jp2`, it has a CRS and a

> transform, so I am not able to reproduce the behavior we have with

> `rasterio`.

>

> (Side note: Interestingly enough, the pixel values of `new_B01.jp2` have

> slightly changed with respect to the original `B01.jp2` file, so there is

> some data loss somewhere in the process. But maybe that is expected of the

> JP2 format and could be avoided by passing extra arguments to

> `gdal.Translate`?

 

By default, the JP2OpenJPEG driver uses a lossy compression (my personal take on JPEG2000 is that it has very limited interested when used in its lossless profile. Better use more conventional algorithms that are way faster).

See https://gdal.org/drivers/raster/jp2openjpeg.html#lossless-compression for options to pass to get a lossless JPEG2000 as output.

 

So all in all that should be something like

 

gdal.Translate(mem, "B01.jp2",

format = "JP2OpenJPEG", # in case you have several JP2K drivers available

options = ["REVERSIBLE=YES", "QUALITY=100"])

 

 

--

Spatialys - Geospatial professional services

http://www.spatialys.com


Re: Rasterio.features.geometry_window error when geojson has GeometryCollection

Brendan Ward
 

This looks like an oversight when we implemented bounds().

Rasterio currently supports FeatureCollection or individual geometry objects in this operation.

Logged here: https://github.com/mapbox/rasterio/issues/1914

Working on a quick fix now...