Creating DEM from mesh?

I have a mesh ( with points and triangles that define it), and I will need to create a DEM model out of it ( in geotiff format if possible for space consideration, if not, then in ASCII grid format). I need the DEM for watershed and catchment delineation purpose.

Instead of using existing software like ArcGIS, I want to write my code to do it ( in C++ or C#), so that I can later integrate them into my own GIS software. is there anything in GDAL or other open source source code libraries that do this well? Or do I have to write from scratch?

PS:

1. there are a lot of posts like Seeking tool to generate Mesh from DTM? on how to do it the other way round, ie: from DEM grid points to mesh, but they are not what I want.


2. I have a similar question as Getting grid points for DEM purpose that is bounded by surveyor points? on how to get DEM grid points from point clouds, but that is fundamentally different from this one. I don't need the mesh interpolation part for this question, as I already have the mesh.


3. I can always create the bounding box for the mesh, loop through each grid points, get the z, and then construct the grid, but this is crude and can be slow. I'm unsure whether there is a faster way of doing it?

Answer

Something similar to Python Scipy will do the job. I can share one of my scripts which uses Scipy to interpolate a raster image from vertices ans saves it as a geotiff file.

Extracting vertices from the mesh should be easy, I just load a csv file to x, y, z variables where each one is an array:

import numpy
import scipy.interpolate as inp

from osgeo import gdal , osr

x, y, z = numpy.loadtxt ( ’ points.csvv ’ , skiprows =1 ,
delimiter = " , " , unpack = True )

Next some useful statistics:

xmin = min ( x )
xmax = max ( x )
ymin = min ( y )
ymax = max ( y )


# number of pixels with 1m resolution
nx = (int(xmax - xmin + 1))
ny = (int(ymax - ymin + 1))

Now a grid, which is also a bounding box, it is represented by nx pixels between xmin and xmax and the same for Y direction:

xi = numpy.linspace(xmin, xmax, nx)
yi = numpy.linspace(ymin, ymax, ny)
xi, yi = numpy.meshgrid(xi, yi) 

And the most important function, in Scipy it is just simple like that but if you can't find any library for C++ or C# you can always write it from scratch, some of these interpolations are quite simple (ex. inverse distance weighting).

# used method: nearest neighbour; there are also cubic and something else
zi = inp.griddata((x, y), z, (xi, yi), method='nearest')

zi variable represents the corresponding Z values so all we have to do is write it to a file (geotiff in this case).

rows,cols = numpy.shape(zi)
sizex = (xmax-xmin)/float(cols)
sizey = (ymax-ymin)/float(rows)

driver = gdal.GetDriverByName('GTiff')

output_raster = driver.Create('raster.tif', rows, cols, 1, gdal.GDT_Float32)

output_raster.GetRasterBand(1).WriteArray(zi)

# georeference
georeference = (xmin, sizex, 0, ymin, 0, sizey) 
output_raster.SetGeoTransform(georeference)
srs = osr.SpatialReference().ImportFromEPSG(2180)
output_raster.SetProjection(srs.ExportToWkt())

And thats all, the example raster.tif file created by this script (zmin red, zmax green):

Here are the docs, you can check out these methods, everything is described:

https://docs.scipy.org/doc/