I have a DEM that is slightly too short for one of my shapefiles. I would like to extend the DEM to the North so it covers entirely the shapefile, and then clip the DEM to the extent of the shapefile. I would like to have the "filled" artificial values of the DEM to be equal to 2000. How could I do that ?
I found a way to do it on Python: (PS: I USED A DIFFERENT DEM FOR THE EXAMPLE, TO DEMONSTRATE HOW TO DEAL WITH BORDERS)
# Open the DEM and extract all the info about it src = rio.open('IfSAR_merged.tif') data = src.read() height = data.shape width = data.shape cols, rows = np.meshgrid(np.arange(width), np.arange(height)) xs, ys = rio.transform.xy(src.transform, rows, cols) lons= np.array(xs) lats = np.array(ys) out_meta = src.meta
# There is a "border" in our DEM, we want to fill the northernmost part of it with 2000s (we can do that because we know there are mountains and no ocean there) # Recover the lowest row of our northermost border filled with 0s for i in range(0, 400): for j in range(0,data.shape): if data[0,i,j] == 0: data[0,i,j] = 2000
# Calculate the step in each direction step_lon = lons[0,1] - lons[0,0] # We don't modify the amount or columns in our matrix, just the rows (lats) min_lon = lons[0,0] # Latitude step is negative because we put the northernmost latitude in our metadata step_lat = lats[1,0] - lats[0,1] # Our new DEM is "higher" than the previous one, so we have to modify the highest latitude in the metadata (we added 1000 rows) max_lat = lats[0,0] + 1000*(-step_lat) # We add 1000 rows to our DEM extension = np.ones((1,1000,data.shape))*2000 new_dem = np.hstack((extension, data)) # Modify the metadata with the new dimensions, modify the Affine function out_meta['width'] = new_dem.shape out_meta['height'] = new_dem.shape out_meta['transform'] = Affine(step_lon, 0.0, min_lon, 0.0, step_lat, max_lat) # Write our new DEM with rio.open('Extended_DEM.tif', "w", **out_meta) as dest: dest.write(new_dem)