1

This question already has an answer here:

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?

marked as duplicate by PolyGeo Dec 13 '18 at 1:26

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • I managed to interpolate rasters from csv file (but technically there is no difference between these and a set of vertices extracted from the mesh) using Python Scipy, which has nice interpolation modules. If you are interested I can share the code later. – dmh126 Dec 13 '18 at 0:34
  • @dmh126 , The CSV file is the point cloud, or already in mesh format? Anyway in both cases I will be interested to see the code, regardless of language – Graviton Dec 13 '18 at 0:38
  • Point cloud. Just X,Y,Z. – dmh126 Dec 13 '18 at 0:40
  • @dmh126 , in that case, can you post the answer to this question? – Graviton Dec 13 '18 at 0:51
  • Okay, I posted here and also pasted it to your other question. Just close one of them because it's a duplicate. If you got any question I can answer them tomorrow. – dmh126 Dec 13 '18 at 1:03
1

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): enter image description here

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

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

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