Get the dem.coords in dem.shape / add a generated raster to existing grid and save it as tif

[ShervanGharari]

Hello,

This might look a very simple question, but I did not really find an answer for.

I have a dem with 6,000 rows and columns. However, the dem.coords return an array of 36,000,000 by 2 which is a pair of the coordinate of each grid on the DEM.

My question: is there any way to get the dem.coords in exact same shape of dem.shape for lats and longs. If a solution doesn't exits within the package, how can I reshape the coordinates to get two array with size of dem.shape for lats and longs. It seems the dem.coords start from top left and are arranges row by row.

Thank you in advance.

PS. I would like to thank you for this very straightforward and nice package. I was working with HAND for hydrological modelling during my master and PhD (for example, https://doi.org/10.5194/hess-15-3275-2011). It is beautiful to see day(s) of coding 10 years ago is now readily accommodated in this package.

[mdbartos]

Greetings,

Glad this has been helpful for you.

I believe you can just do:

dem_coords = dem.coords
dem_coords_reshaped = [dem_coords[:,0].reshape(dem.shape),
                       dem_coords[:,1].reshape(dem.shape)]

For verification, see:

import numpy as np

n = 10
x = np.arange(n)
y = np.arange(n)
shape = (n, n)

xy = np.meshgrid(x, y)
xy_flat = np.vstack(np.dstack(xy))
xy_reshaped = [xy_flat[:,0].reshape(shape), xy_flat[:,1].reshape(shape)]
assert((i == j).all() for i, j in zip(xy, xy_reshaped))

[ShervanGharari]

Awesome, thank you for the clarification. Maybe I can ask an additional question which is somehow related to the coordinates of the points. How can I write a Numpy Array with the same shape of the input raster as tif file with pysheds? In this example, how can I save z1 a tiff file with the same information as input dem. I know how to do this with gdal, basically reading the DEM, passing all the information to a new file replacing the values with z1 but was just wondering if there is quick solution for this with pysheds. Thank you again.

[mdbartos]

Greetings. Because z1 has the same coordinates as the current grid, you can just do:

# Add z1 to grid
grid.add_gridded_data(z1, data_name='z1', affine=grid.affine,
                      shape=z1.shape, crs=grid.crs, nodata=np.nan)

# Write to raster
grid.to_raster('z1', 'z1.tif', view=False)

# Read new raster
grid.read_raster('z1.tif', data_name='z1_copy')

# Check that they are the same
assert np.where(~np.isnan(grid.z1) & ~np.isnan(grid.z1_copy),
         grid.z1 == grid.z1_copy, True).all()

If the affine transform, crs, or shape are different, it gets a bit more complicated.

[ShervanGharari]

Awesome! thank you very much. I would suggest you can maybe add this to the example on snap points. would be a good conclusion for that exercise to save the sub-basins. I have also encounter an strange issue which seems to be reported here. This result in the final sub-basin to have small holes especially on very flat land or large lakes. Is there any way to solve this manually? Sorry if I am asking very broad topics under the same issue, but hopefully it helps others.

[ShervanGharari]

Greetings again! How can I save a raster which is int similar to flow direction for example. grid.to_raster('dir', '/Users/shg096/Desktop/DEM_practice/dir.tif', view=False)

[ShervanGharari]

Another question, if I may ask. For the issue with resolve_flats (https://github.com/mdbartos/pysheds/issues/90), I would like to do the following: find the resolve_flats raster values that are missing while they are present in dem. replace the missing values with very high values to ensure they drain outward. I think I can do this manipulation with numpy. Just wondering if there is a better way for doing this. thank you in advance.

grid.flowdir(data='inflated_dem', out_name='dir', dirmap=dirmap)

# find the direction that is set to -1 and the DEM is not the missing value (missing value here -9999)
# those cells are orrupted by inflated_dem and we set them to a very high value in order to fix them
temp=(np.array(grid.dir)==-1)*(np.array(grid.dem)!=-9999)
ID = np.where(temp==1) # location where the DEM value exists but corrupted by resolves_flats
inflated_dem_new = np.array(grid.inflated_dem)
inflated_dem_new[ID] = 10000 # an arbitary high value in m
# Add new inflated_dem to grid
grid.add_gridded_data(inflated_dem_new, data_name='inflated_dem_new', affine=grid.affine,
                      shape=grid.dem.shape, crs=grid.crs, nodata=np.nan)

# Again, specifying flow direction given the directional map
grid.flowdir(data='inflated_dem_new', out_name='dir', dirmap=dirmap)```