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?

Graviton
  • 928
  • 9
  • 28
  • You can do this in QGIS https://gis.stackexchange.com/questions/118243/create-a-dem-from-point-data, if you're keen you can look at the source which will be in C or C++. – Michael Stimson Dec 03 '18 at 02:15

2 Answers2

2

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/

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
dmh126
  • 6,732
  • 2
  • 21
  • 36
1

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
Graviton
  • 928
  • 9
  • 28