4

I am quite new to GDAL and Python but want to use it to create groundwater models in FloPy.

I have a raster with rivers (roughly drawn) and a DEM raster. The resolution of the river-raster is higher than that of the DEM. Thereby the river-raster only covers a part of the total area of the DEM-map:

raster with roughly the rivers and the DEM in the back

What I want is the following: 1) Apply the resolution of the DEM to the river-raster 2) Extend the river-raster so it has the same extent as the DEM by filling it with no-data values within the boundaries of the DEM.

In QGIS it works easily by saving the river raster with a different extent and resolution but since I want to do this for whole batches it would be nice if I learn how to do this with GDAL in Python.

2

There's another way to do this, also. You can load in your DEM as a numpy array using

dem_obj = gdal.Open(path_to_DEM)
I = dem_obj.ReadAsArray()

Now set all the values equal to nan:

I = I * np.nan

Now you need to load your rivers raster:

rr_obj = gdal.Open(path_to_riv_raster)
Ir = rr_obj.ReadAsArray()

I would resample the river raster next, using skimage.transform's resize function. You need to know the size of your resampled river raster, though, which you can obtain through the geotransforms of both:

gt_dem = dem_obj.GetGeoTransform()
new_rr_shape = (Ir.shape[0]/gt_dem[4], Ir.shape[1]/gt_dem[1])

Now resize the river raster so that the pixels are the same size as the DEM:

import skimage.transform as st 
Ir = st.resize(Ir, new_rr_shape, mode='constant')

Now you have your river raster at the same resolution as your DEM. The next step is to place it in the correct position within the np.nan array we created from the DEM. To do that, you need to use the geotransforms of both images to align the top-left-most pixel.

gt_rr = rr_obj.GetGeoTransform()
ncols_offset = (gt_dem[0] - gt_rr[0]) / gt_dem[1]
nrows_offset = (gt_dem[3] - gt_rr[3]) / gt_dem[4]

You'll notice that your ncols and nrows _offsets are not integers, but you'll need them to be. You can use either round() or int() (or other choices) to make them integers; the choice you make will shift the river raster no more than a pixel in either direction.

Now that you have your river raster resampled to the same resolution, and an empty (nan-filled, really) DEM template raster to put them into, and the row, col offsets, you can just put the river raster in with numpy:

import numpy as np
I[nrows_offset:nrows_offset+Ir.shape[1], ncols_offset:ncols_offset+Ir.shape[0]] = Ir

Finally, you'll want to save your resampled, rescaled river raster as a geotiff. I will post code here you can use as a template:

def write_geotiff(raster, gt, wkt, outputpath, dtype=gdal.GDT_UInt16, options=['COMPRESS=LZW'], color_table=0, nbands=1, nodata=False):

    width = np.shape(raster)[1]
    height = np.shape(raster)[0]

    # Prepare destination file
    driver = gdal.GetDriverByName("GTiff")
    if options != 0:
        dest = driver.Create(outputpath, width, height, nbands, dtype, options)
    else:
        dest = driver.Create(outputpath, width, height, nbands, dtype)

    # Write output raster
    if color_table != 0:
        dest.GetRasterBand(1).SetColorTable(color_table)

    dest.GetRasterBand(1).WriteArray(raster)

    if nodata is not False:
        dest.GetRasterBand(1).SetNoDataValue(nodata)

    # Set transform and projection
    dest.SetGeoTransform(gt)
    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    dest.SetProjection(srs.ExportToWkt())

    # Close output raster dataset 
    dest = None

Your call should look something like:

write_geotiff(I, gt_dem, dem_obj.GetProjection(), path_to_output_file, nodata=np.nan)
  • Disclaimer: I didn't test any of this, so there's a chance I switched rows/cols or missed adding a 1 when putting the resampled river raster into the nan raster. But this is a workflow you could use to accomplish what you want.

  • An alternative method would have you first pad your river raster so that its extents matched the DEM's extents, then resample it to the same grid size as the DEM. This would take more memory, but it might align the rasters slightly better. All the tools you need to do it that way are given above, you just need to reorder them and change the inputs accordingly.

0

One possible solution is to create an empty raster using the DEM as the template and gdal_calc (see here). This will ensure the same resolution, geotransform, pixel alignment, etc. Then polygonize your river raster using gdal_polygonize.py, and finally burn the polygons back into the empty raster with gdal_rasterize.

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.