1

I have a list of surveyor points, 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 tools that do this well? Or do I have to write from scratch?

1

I'm pasting my answer from Creating DEM from mesh?

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/

0

Actually there is a more direct way of doing it, by using gdal_grid ( in C++ command line), or Gdal.wrapper_GDALGrid ( in GDAL.Net). Here's the code sample:

string vrtFile =@"vrtFile.vrt";
string tiffFile= @"rasterized.tif";
var ds = Gdal.OpenEx(vrtFile, 0, null, null, null);
Gdal.wrapper_GDALGrid(tiffFile, ds, null, null, string.Empty));

A sample of the content in vrt file:

<OGRVRTDataSource>
<OGRVRTLayer name="vrtFile">
    <SrcDataSource>vrtFile.csv</SrcDataSource>
    <GeometryType>wkbPoint</GeometryType>
    <GeometryField encoding="PointFromColumns" x="long" y="lat" z="depth"/>
</OGRVRTLayer>
</OGRVRTDataSource>

A sample of content in vrtFile.csv, this is your scattered point list

long,lat,depth
565650,5121960,1048
565680,5121960,1043  

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.