mdbartos / pysheds

:earth_americas: Simple and fast watershed delineation in python.
GNU General Public License v3.0
708 stars 191 forks source link

loading data from AWS #191

Open cmacritchie opened 2 years ago

cmacritchie commented 2 years ago

I have a DEM file stored in an S3 bucket. I'm reading the file using rasterio and the data file is type: numpy.ndarray. I tried following the suggestion from "Problem with loading DEM as numpy array #99" but was unable add the grid data, is add_gridded_data deprecated? Is there a way to grab the data from an S3 bucket (without downloading it locally) and use it in pysheds? Here is what I currently have:

`src = rasterio.open('/vsis3/link-to-aws-DEM.tif') data = src.read(1) # type: numpy.ndarray dem = gaussian_filter(data,sigma=8)

affine = Affine(1,0,0,0,1,0) grid = Grid() grid.add_gridded_data(dem,'dem',affine=affine, crs={'init':'epsg:32614'}) dem_view = grid.view('dem',apply_mask=False) ` The error returns: 'sGrid' object has no attribute 'add_gridded_data'

Thank you!

mdbartos commented 2 years ago

Yes, add_gridded_data is deprecated as of v0.3. Also, datasets as named parameters of grid have been removed. All grid methods now operate on Raster objects, which are essentially ndarrays with added spatial referencing metadata.

This example shows how you can construct a Raster from an ndarray: https://mattbartos.com/pysheds/raster.html#instantiating-rasters

cmacritchie commented 2 years ago

This is great, thanks Matt! this ended up working, but I couldn't seem to maintain the latitude/ longitude in the raster which made it difficult to delineate a watershed when selecting the outlet based on lat/lng coordinates. Do you have any suggestions for maintaining latlng after converting to a raster from numpy.ndarray?

mdbartos commented 2 years ago

To get the spatial referencing to work out, you will need to know the lat/lon coordinates (in degrees) of the upper left corner of the upper left pixel of your array, and the width and height of each grid cell (in lat/lon degrees).

These coefficients then go into the affine transformation matrix, where:

a: width of pixel c: x-coordinate of upper-left corner of upper-left pixel e: height of pixel f: y-coordinate of upper-left corner of upper-left pixel

(See: https://www.perrygeo.com/python-affine-transforms.html)

So let's say your pixels are 0.5 deg wide and 0.5 deg tall, with an upper left-hand corner of lat: 32.2, lon=-97.3. Then you would do:

affine = Affine(0.5, 0, -97.3, 0, -0.5, 32.2)

Let me know if that helps.