Date   

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

thomaswanderer@...
 

Thank you for the info. I saw shorter values tan in QGIS and those came from Python's float limitations. When using the Decimal class I also get the same over-precise values.


Get bounds after rasterio merge

ashnair0007@...
 

I merged two tiffs via the rasterio.merge function which returned the merged numpy array and the geo transform. Is there a way to get the bounds (in latitude and longitude) of this new merged image?


can resampling not populate cell if source has nodata?

Nick Forfinski-Sarkozi
 

Hi rasterio.  New user here.  Thanks for a fantastic module.  During resampling, is there a way to force the cells of the resampled raster to nodata if any of the corresponding original cells are nodata?   (If by chance anyone is familiar with it, I'm looking for functionality analogous to arcpy's Aggregate function with ignore_data=True https://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/aggregate.htm )

Thanks,
Nick


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

Sean Gillies
 

Hi,

On Mon, May 25, 2020 at 3:30 AM <thomaswanderer@...> wrote:

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.

QGIS and rasterio both use GDAL to read raster data and all of these use double precision floats to store values of the georeferencing matrix. I believe you're seeing different string representations of the same numbers. Since the hex representations of the two numbers are the same, I suspect the values shown by QGIS are overly precise and that the last digits are meaningless.

>>> 0.0008333333299974727193.hex()
'0x1.b4e81b312a043p-11'
>>> 0.0008333333299974727.hex()
'0x1.b4e81b312a043p-11'

--
Sean Gillies


Re: Create raster aligned with one created by gdal_translate

thomaswanderer@...
 

I also still fight with the shifted by one pixel problem for 1 extent out of 5 extents. 4 return perfectly aligned clipped subsets of the source raster, 1 does not.


Re: Create raster aligned with one created by gdal_translate

Denis Rykov
 

Thanks Thomas. I've tried your approach but it still gives me a shifted output compared to gdal_translate.


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.

161 - 180 of 698