1

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 ?

DEM too short

Metadata shapefile: enter image description here

Metadata Raster: enter image description here

2
  • You can't extend a DEM because there is no data to extend with. Have you tried merging the DEM to the North and then cropping?
    – Pointdump
    May 16 at 8:27
  • I tried but unfortunately it give me the same error. I kind of managed to "extend" the DEM to the North by opening it as an array on Python, extending this array by a 1000 rows to the North, and filling them with values superior to 2000. I'm cleaning the code to put it as an answer in case somebody needs that
    – vdc
    May 16 at 17:13

1 Answer 1

1

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[1]
width = data.shape[2]
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

DEM with borders

# 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[2]):
        if data[0,i,j] == 0:
            data[0,i,j] = 2000

Dem without border in the North

# 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[2]))*2000
new_dem = np.hstack((extension, data))

# Modify the metadata with the new dimensions, modify the Affine function
out_meta['width'] = new_dem.shape[2]
out_meta['height'] = new_dem.shape[1]
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)

What we obtain is a DEM that can be clipped to the extents of the raster. New DEM with fake data in the North

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.