# Getting grid points for DEM purpose that is bounded by surveyor points?

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?

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): Here are the docs, you can check out these methods, everything is described:

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

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
``````