Friday, 10 March 2017

Get a window from a raster in rasterio using coordinates instead of row/column offsets, width and height


I am trying to open a windowed dataset in rasterio, but I have a bounding box to work with. I thought it would be easy to do this but it seems that it involves some complicated process with a "Mixin", something I am having a very hard time understanding, even after reading the thread here.


The documentation says:



A subclass with this mixin MUST provide the following properties: transform, height and width




But I don't know how I can implement these properties if I don't know ahead of time what the height and width are going to be, and I'm not even sure that I am implementing the code properly anyway. This is what I have so far:


import rasterio
from rasterio.windows import WindowMethodsMixin, Window
from rasterio.enums import Resampling

class MyWindow(WindowMethodsMixin, Window):
pass

with rasterio.open("flask/Docker/app/dem/slope_sm.tif") as src:

rst1 = src.read(1, window=MyWindow.window(...))

I know this is hardly complete but I am really confused about how to proceed.


My IDE tells me that the call to MyWindow.window() takes the following parameters:



  1. self

  2. left

  3. bottom

  4. right

  5. top


  6. precision (optional)


But I don't know what to pass to the "self" parameter. The rest of the parameters are understandable - they are the edges of my bounding box from which I am constructing the window.


The documentation doesn't supply any working example of this, which I find odd because this seems to me like a pretty typical operation.


Can somebody please show me how to do this, or at least explain what I am missing here?



Answer



Use the rasterio.windows.from_bounds function. No need for a class or mixin.


import rasterio
from rasterio.windows import from_bounds
from rasterio.enums import Resampling



with rasterio.open(filepath) as src:
rst = src.read(1, window=from_bounds(left, bottom, right, top, src.transform))

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...