Date   
Re: Window from bounds is shifted by 1 Pixel

Sean Gillies
 

Hi Thomas,

I think Denis found the root of the problem and we have a fix at https://github.com/mapbox/rasterio/pull/1938 that will be a part of the 1.1.5 release. Thank you for the test cases and for your patience!


On Tue, May 26, 2020 at 4:39 PM <thomaswanderer@...> wrote:

I am not 100 sure, but I now simply rounded all windows (origin + offsets) and this seems to have fixed it.



--
Sean Gillies

Re: Create raster aligned with one created by gdal_translate

Denis Rykov
 

Re: can resampling not populate cell if source has nodata?

Nick Forfinski-Sarkozi
 

FYI...I was able to achieve the desired behavior by "lying" about the nodata value.  The actual nodata was nan, but when I specified something else (e.g., nodata=-9999) in the profile, the resampled raster didn't contain values for cells corresponding to source nans.


On Wed, May 27, 2020 at 4:20 PM Sean Gillies via groups.io <sean=mapbox.com@groups.io> wrote:
I see what you mean. Rasterio doesn't have a built in option for this, unless there is a GDAL resampling trick that I am not aware of.

On Wed, May 27, 2020 at 9:55 AM <nick.forfinski@...> wrote:
Hi Sean.  Thanks for the message.  I understand your example, but I may have been confusing in my email.  I think I'm asking something different.  Is there a way to force the resampling process to not populate a resampled value where there are corresponding nodata values in the original raster?  For example, in the image below, the lighter 5-m purple raster is the resampled version of the darker 1-m blue raster.  Many of the resampled raster's pixels correspond to gaps in the original raster.  How could I resample so that all 25 source cells are required to produce a resampled value?  
 
Thanks again for your help and a great module,
Nick
 



--
Sean Gillies

Re: can resampling not populate cell if source has nodata?

Alan Snow
 

Here is an example using rioxarray: https://corteva.github.io/rioxarray/stable/examples/interpolate_na.html

It uses scipy griddata.

Re: can resampling not populate cell if source has nodata?

Sean Gillies
 

I see what you mean. Rasterio doesn't have a built in option for this, unless there is a GDAL resampling trick that I am not aware of.

On Wed, May 27, 2020 at 9:55 AM <nick.forfinski@...> wrote:
Hi Sean.  Thanks for the message.  I understand your example, but I may have been confusing in my email.  I think I'm asking something different.  Is there a way to force the resampling process to not populate a resampled value where there are corresponding nodata values in the original raster?  For example, in the image below, the lighter 5-m purple raster is the resampled version of the darker 1-m blue raster.  Many of the resampled raster's pixels correspond to gaps in the original raster.  How could I resample so that all 25 source cells are required to produce a resampled value?  
 
Thanks again for your help and a great module,
Nick
 



--
Sean Gillies

Re: can resampling not populate cell if source has nodata?

Nick Forfinski-Sarkozi
 

Hi Sean.  Thanks for the message.  I understand your example, but I may have been confusing in my email.  I think I'm asking something different.  Is there a way to force the resampling process to not populate a resampled value where there are corresponding nodata values in the original raster?  For example, in the image below, the lighter 5-m purple raster is the resampled version of the darker 1-m blue raster.  Many of the resampled raster's pixels correspond to gaps in the original raster.  How could I resample so that all 25 source cells are required to produce a resampled value?  
 
Thanks again for your help and a great module,
Nick
 

Re: can resampling not populate cell if source has nodata?

Sean Gillies
 

Hi Nick,

Ignoring nodata values is the default behavior for rasterio, as it is for the GDAL library that rasterio uses. Here's an example using one of rasterio's test files.

$ rio insp tests/data/RGB.byte.tif
Rasterio 1.1.5dev Interactive Inspector (Python 3.6.6)
Type "src.meta", "src.read(1)", or "help(src)" for more information.
>>> show(src.read())

That command uses matplotlib to display the raw data and I've zoomed in to the bottom to show the hard edge along the image collar. Valid data on the inside, nodata on the outside.

Figure_1.png
Now, I'll dilate the image by a factor of 3 on reading and use bilinear resampling to fill in the extra pixels.

>>> from rasterio.enums import Resampling
>>> resampled = src.read(out_shape=(3, src.height * 3, src.width * 3), resampling=Resampling.bilinear)
>>> show(resampled)

Zooming in again to the collar, you can see that the valid data has been smoothed, but the hard nodata edge remains. Nodata values were ignored by the resampling system.
Figure_2.png

On Tue, May 26, 2020 at 11:30 AM <nick.forfinski@...> wrote:
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



--
Sean Gillies

Re: Get bounds after rasterio merge

Luke
 

On Wed, 27 May 2020 at 04:18, <ashnair0007@...> wrote:
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?

Re: Window from bounds is shifted by 1 Pixel

thomaswanderer@...
 

I am not 100 sure, but I now simply rounded all windows (origin + offsets) and this seems to have fixed it.

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?