Date   

Re: Merge resulting in one pixel shift for half the output file

Denis Rykov
 

Do you use the latest version of rasterio?


Merge resulting in one pixel shift for half the output file

melda.salhab.14@...
 

I'm using merge with method = 'max' to merge two rasters with identical metadata. When I open the output file and the input files in QGIS, the top part of the output raster is perfectly merged, but the bottom part of the output raster seems to have the merged values of the pixel to it's left. Very bizarre. Has anyone experienced anything like this? Code below.


    combinedFloodList = [FluvialCropped_outputfile] + [PluvialCropped_outputfile] 
 
    # Iterate over raster files and add them to list in 'read mode'            
    allFiles = []
    for fp in combinedFloodList:
        src = rs.open(fp)
        allFiles.append(src)
 
    # Merge flood files. Using max pixel by pixel method
    floodArray, out_trans = merge(allFiles, method='max')
 
    # Save results as raster
    ## Update the metadata
    out_meta = allFiles[0].meta.copy()
    out_meta.update({"height": floodArray.shape[1], 
                      "width": floodArray.shape[2],
                      "transform": out_trans})
 
    # Output file name
    floodName = f"Flood_{countryCode}_{adm1Index}.tif"
 
    ## export as a new geotiff 
    FloodMerged_outputfile = Flood_outputFolderPath / floodName
    with rs.open(FloodMerged_outputfile, 'w', **out_meta, compress = 'LZW') as dest:
        dest.write(floodArray)


Different values when I use a window

adrianocorbelinoii@...
 

Hello,  I've been having a problem using windows, probably a misunderstanding on my part.
My goal is to dump a specific data from a band into a csv file, first I tried to use a window but the results differ from the results that I see on Qgis(the results appear to be shifted 1 position).
So I decided to try to read the whole band, instead to use a window, and the result match with Qgis.
I think that I am doing something wrong when I create/use the window, what I tried to do was:
  1. create a window from the bounds of the desired area: 
    bandWin = rasterio.windows.from_bounds(*boundingBox,band.transform)
  2. get the transform matrix of that window:
    win_transform = band.window_transform(bandWin)
  3. Obtain the window data:
    win_data = band.read(1, window = bandWin)
  4. Get the desired window indexes:
    indexes = rasterio.transform.rowcol(win_transform, [Xs], [Ys])
  5. Get the values that I want:
    csvData.append(win_data[indexes[0],indexes[1]][0])
The only differences between when I don't use a window is that instead use win_transform I use band.transform and I read the whole band.

All help is appreciated.


List of numpy arrays to raster

ashnair0007@...
 

I originally had a geotiff (size= 4490 x 7341) which I tiled into tiles of size 256 x 256 to perform an operation. The result is a list of 522 numpy array. From the original geotiff I have the geo transform as well.  What would be the recommended way to merge them to get the raster? 


Re: Window from bounds is shifted by 1 Pixel

thomaswanderer@...
 

Thanks for letting me know! Thumbs up for your work!


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.