Add GDALDEMProcessing binding (#117)

errors.go

@@ -152,6 +152,7 @@ func ErrLogger(fn ErrorHandler) interface {
 	ClearStatisticsOption
 	GridOption
 	NearblackOption
+	DemOption
 	SetGCPsOption
 	RegisterPluginOption
 } {
@@ -374,6 +375,9 @@ func (ec errorCallback) setGridOpt(o *gridOpts) {
 func (ec errorCallback) setNearblackOpt(o *nearBlackOpts) {
 	o.errorHandler = ec.fn
 }
+func (ec errorCallback) setDemOpt(o *demOpts) {
+	o.errorHandler = ec.fn
+}
 func (ec errorCallback) setSetGCPsOpt(o *setGCPsOpts) {
 	o.errorHandler = ec.fn
 }

godal.cpp

@@ -1788,6 +1788,27 @@ GDALDatasetH godalNearblack(cctx *ctx, const char *pszDest, GDALDatasetH hDstDS,
 	return ret;
 }
 
+GDALDatasetH godalDem(cctx *ctx, const char *pszDest, const char *pszProcessing, const char *pszColorFilename, GDALDatasetH hSrcDS, char **switches) {
+	godalWrap(ctx);
+
+	GDALDEMProcessingOptions *demopts = GDALDEMProcessingOptionsNew(switches,nullptr);
+	if(failed(ctx)) {
+		GDALDEMProcessingOptionsFree(demopts);
+		godalUnwrap();
+		return nullptr;
+	}
+
+	int usageErr=0;
+	GDALDatasetH ret = GDALDEMProcessing(pszDest, hSrcDS, pszProcessing, pszColorFilename, demopts, &usageErr);
+	GDALDEMProcessingOptionsFree(demopts);
+	if(ret==nullptr || usageErr!=0) {
+		forceError(ctx);
+	}
+
+	godalUnwrap();
+	return ret;
+}
+
 OGRSpatialReferenceH godalGetGCPSpatialRef(GDALDatasetH hSrcDS) {
 	return GDALGetGCPSpatialRef(hSrcDS);
 }

godal.go

@@ -3877,6 +3877,47 @@ func (ds *Dataset) Grid(destPath string, switches []string, opts ...GridOption)
 	return &Dataset{majorObject{C.GDALMajorObjectH(dsRet)}}, nil
 }
 
+// Dem runs the library version of gdaldem.
+// See the gdaldem doc page to determine the valid flags/opts that can be set in switches.
+//
+// Example switches (for "hillshade", switches differ per mode):
+//
+//	[]string{"-s", "111120", "-alt", "45"}
+//
+// Creation options and driver may be set in the switches slice with
+//
+//	switches:=[]string{"-co","TILED=YES","-of","GTiff"}
+//
+// NOTE: `colorFilename` is a "text-based color configuration file" that MUST ONLY be
+// provided when `processingMode` == "color-relief"
+func (ds *Dataset) Dem(destPath, processingMode string, colorFilename string, switches []string, opts ...DemOption) (*Dataset, error) {
+	demOpts := demOpts{}
+	for _, opt := range opts {
+		opt.setDemOpt(&demOpts)
+	}
+
+	cswitches := sliceToCStringArray(switches)
+	defer cswitches.free()
+
+	dest := unsafe.Pointer(C.CString(destPath))
+	defer C.free(unsafe.Pointer(dest))
+	alg := unsafe.Pointer(C.CString(processingMode))
+	defer C.free(unsafe.Pointer(alg))
+	var colorFn *C.char
+	if colorFilename != "" {
+		colorFn = C.CString(colorFilename)
+		defer C.free(unsafe.Pointer(colorFn))
+	}
+
+	cgc := createCGOContext(nil, demOpts.errorHandler)
+	dsRet := C.godalDem(cgc.cPointer(), (*C.char)(dest), (*C.char)(alg), colorFn, ds.handle(), cswitches.cPointer())
+	if err := cgc.close(); err != nil {
+		return nil, err
+	}
+
+	return &Dataset{majorObject{C.GDALMajorObjectH(dsRet)}}, nil
+}
+
 // Nearblack runs the library version of nearblack
 //
 // See the nearblack doc page to determine the valid flags/opts that can be set in switches.

godal.h

@@ -156,6 +156,7 @@ extern "C" {
 	void godalGridCreate(cctx *ctx, char *pszAlgorithm, GDALGridAlgorithm eAlgorithm, GUInt32 nPoints, const double *padfX, const double *padfY, const double *padfZ, double dfXMin, double dfXMax, double dfYMin, double dfYMax, GUInt32 nXSize, GUInt32 nYSize, GDALDataType eType, void *pData);
 	GDALDatasetH godalGrid(cctx *ctx, const char *pszDest, GDALDatasetH hSrcDS, char **switches);
 	GDALDatasetH godalNearblack(cctx *ctx, const char *pszDest, GDALDatasetH hDstDS, GDALDatasetH hSrcDS, char **switches);
+	GDALDatasetH godalDem(cctx *ctx, const char *pszDest, const char *pszProcessing, const char *pszColorFilename, GDALDatasetH hSrcDS, char **switches);
 
 	typedef struct {
 		const GDAL_GCP *gcpList;

godal_test.go

@@ -4798,3 +4798,285 @@ func TestSetGCPs2InvalidDataset(t *testing.T) {
 	err = vrtDs.SetGCPs([]GCP{}, GCPSpatialRef(&SpatialRef{}))
 	assert.Error(t, err)
 }
+
+func TestDemHillshade(t *testing.T) {
+	// 1. Create an image, linearly interpolated, from dark (on the left) to white (on the right), using `Grid()`
+	var (
+		outXSize = 2048
+		outYSize = 2048
+	)
+	vrtDs, err := CreateVector(Memory, "")
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	defer vrtDs.Close()
+	geom, err := NewGeometryFromWKT("POLYGON((500000 500000 10, 500000 600000 10, 600000 600000 2026, 600000 500000 2026))", nil)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	defer geom.Close()
+	_, err = vrtDs.CreateLayer("grid", nil, GTPolygon)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	_, err = vrtDs.Layers()[0].NewFeature(geom)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+
+	// As of GDAL v3.6, `GDALGrid` will swap `yMin` and `yMax` if `yMin` < `yMax`. In order to make the output of
+	// earlier GDAL versions (< 3.6) consistent with this, we're setting `yMin` > `yMax`.
+	yMin := 600000
+	yMax := 500000
+	argsString := fmt.Sprintf("-a linear -txe 500000 600000 -tye %d %d -outsize %d %d -ot Int16", yMin, yMax, outXSize, outYSize)
+	fname := "/vsimem/grid.tiff"
+	gridDs, err := vrtDs.Grid(fname, strings.Split(argsString, " "))
+	if err != nil {
+		// Handles QHull error differently here, as it's a compatibility issue not a gridding error
+		isQhullError := strings.HasSuffix(err.Error(), "without QHull support")
+		if isQhullError {
+			t.Log(`Skipping test, GDAL was built without "Delaunay triangulation" support which is required for the "Linear" gridding algorithm`)
+			return
+		} else {
+			t.Error(err)
+			return
+		}
+	}
+	gridDs.SetProjection("epsg:32632")
+	defer func() { _ = VSIUnlink(fname) }()
+	defer gridDs.Close()
+
+	fname2 := "/vsimem/dem.tif"
+	demDs, err := gridDs.Dem(fname2, "hillshade", "", []string{})
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	defer func() { _ = VSIUnlink(fname2) }()
+	defer demDs.Close()
+
+	demBuf := make([]uint8, outXSize*outYSize)
+	err = demDs.Read(0, 0, demBuf, outXSize, outYSize)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+
+	// Expected value of border is NODATA (0): https://gdal.org/programs/gdaldem.html
+	var nodata uint8 = 0
+	assert.Equal(t, nodata, demBuf[0])
+	assert.Equal(t, nodata, demBuf[outXSize-1])
+	assert.Equal(t, nodata, demBuf[(outXSize*(outYSize-1))])
+	assert.Equal(t, nodata, demBuf[(outXSize*outYSize)-1])
+
+	// Values in each "column" should be equal
+	// NOTE: Skips the outer pixel border since it's set to the NODATA value
+	for x := 1; x < 256; x++ {
+		var topValInColumn = demBuf[(1*outXSize)+x]
+		for y := 1; y < 2047; y++ {
+			thisColVal := demBuf[(y*outXSize)+x]
+			assert.Equal(t, topValInColumn, thisColVal)
+		}
+	}
+
+	// Iterate through lines on the middle row of the output dataset checking that
+	// each line has equal spacing between it and the next line
+	var (
+		expLineThickness    = 2
+		thisLineThickness   = 0
+		expInterLineSpaces  = 62
+		thisInterLineSpaces = 0
+		row                 = (outYSize / 2)
+
+		expSpaceVal uint8 = 183
+		expLineVal  uint8 = 182
+	)
+	for x := 1; x < outXSize-1; x++ {
+		thisCoordVal := demBuf[(row*outYSize)+x]
+		switch thisCoordVal {
+		case expSpaceVal:
+			if thisLineThickness > 0 {
+				assert.Equal(t, expLineThickness, thisLineThickness)
+				thisLineThickness = 0
+			}
+			thisInterLineSpaces++
+		case expLineVal:
+			if thisInterLineSpaces > 0 {
+				assert.Equal(t, expInterLineSpaces, thisInterLineSpaces)
+				thisInterLineSpaces = 0
+			}
+			thisLineThickness++
+		default:
+			t.Errorf("found coordinate with value not in: [%d, %d]", expSpaceVal, expLineVal)
+			return
+		}
+	}
+}
+
+func TestDemInvalidSwitch(t *testing.T) {
+	fname := "/vsimem/test.tiff"
+	vrtDs, err := Create(Memory, fname, 1, Byte, 256, 256)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	defer func() { _ = VSIUnlink(fname) }()
+	defer vrtDs.Close()
+
+	_, err = vrtDs.Dem("test", "hillshade", "", []string{"-invalidswitch"})
+	assert.Error(t, err)
+	ehc := eh()
+	_, err = vrtDs.Dem("test", "hillshade", "", []string{"-invalidswitch"}, ErrLogger(ehc.ErrorHandler))
+	assert.Error(t, err)
+}
+
+func TestDemSlope(t *testing.T) {
+	// 1. Create an image, linearly interpolated, from black (on the left) to white (on the right), using `Grid()`
+	var (
+		outXSize = 2048
+		outYSize = 2048
+	)
+	vrtDs, err := CreateVector(Memory, "")
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	defer vrtDs.Close()
+	geom, err := NewGeometryFromWKT("POLYGON((500000 500000 10, 500000 600000 10, 600000 600000 2026, 600000 500000 2026))", nil)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	defer geom.Close()
+	_, err = vrtDs.CreateLayer("grid", nil, GTPolygon)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	_, err = vrtDs.Layers()[0].NewFeature(geom)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+
+	// As of GDAL v3.6, `GDALGrid` will swap `yMin` and `yMax` if `yMin` < `yMax`. In order to make the output of
+	// earlier GDAL versions (< 3.6) consistent with this, we're setting `yMin` > `yMax`.
+	yMin := 600000
+	yMax := 500000
+	argsString := fmt.Sprintf("-a linear -txe 500000 600000 -tye %d %d -outsize %d %d -ot Int16", yMin, yMax, outXSize, outYSize)
+	fname := "/vsimem/grid.tiff"
+	gridDs, err := vrtDs.Grid(fname, strings.Split(argsString, " "))
+	if err != nil {
+		// Handles QHull error differently here, as it's a compatibility issue not a gridding error
+		isQhullError := strings.HasSuffix(err.Error(), "without QHull support")
+		if isQhullError {
+			t.Log(`Skipping test, GDAL was built without "Delaunay triangulation" support which is required for the "Linear" gridding algorithm`)
+			return
+		} else {
+			t.Error(err)
+			return
+		}
+	}
+	gridDs.SetProjection("epsg:32632")
+	defer func() { _ = VSIUnlink(fname) }()
+	defer gridDs.Close()
+
+	fname2 := "/vsimem/dem.tif"
+	demDs, err := gridDs.Dem(fname2, "slope", "", []string{"-p"})
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	defer func() { _ = VSIUnlink(fname2) }()
+	defer demDs.Close()
+
+	demBuf := make([]float32, outXSize*outYSize)
+	err = demDs.Read(0, 0, demBuf, outXSize, outYSize)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+
+	// Expected value of border is NODATA (0): https://gdal.org/programs/gdaldem.html
+	var nodata float32 = -9999
+	assert.Equal(t, nodata, demBuf[0])
+	assert.Equal(t, nodata, demBuf[outXSize-1])
+	assert.Equal(t, nodata, demBuf[(outXSize*(outYSize-1))])
+	assert.Equal(t, nodata, demBuf[(outXSize*outYSize)-1])
+
+	// Values in each "column" should be equal
+	// NOTE: Skips the outer pixel border since it's set to the NODATA value
+	for x := 1; x < 256; x++ {
+		var firstVal = demBuf[(1*outXSize)+x]
+		for y := 1; y < 2047; y++ {
+			thisCoord := demBuf[(y*outXSize)+x]
+			assert.Equal(t, firstVal, thisCoord)
+		}
+	}
+
+	// Iterate through lines on the middle row of the output dataset checking that
+	// each line has equal spacing between it and the next line
+	var (
+		expLineThickness    = 2
+		thisLineThickness   = 0
+		expInterLineSpaces  = 62
+		thisInterLineSpaces = 0
+		row                 = (outYSize / 2)
+
+		expSpaceVal float32 = 2.048
+		expLineVal  float32 = 1.024
+	)
+	for x := 1; x < outXSize-1; x++ {
+		thisCoordVal := demBuf[(row*outYSize)+x]
+		switch thisCoordVal {
+		case expSpaceVal:
+			if thisLineThickness > 0 {
+				assert.Equal(t, expLineThickness, thisLineThickness)
+				thisLineThickness = 0
+			}
+			thisInterLineSpaces++
+		case expLineVal:
+			if thisInterLineSpaces > 0 {
+				assert.Equal(t, expInterLineSpaces, thisInterLineSpaces)
+				thisInterLineSpaces = 0
+			}
+			thisLineThickness++
+		default:
+			t.Errorf("found coordinate with value not in: [%f, %f]", expSpaceVal, expLineVal)
+			return
+		}
+	}
+}
+
+func TestDemColorReliefNilFilename(t *testing.T) {
+	fname := "/vsimem/test.tiff"
+	vrtDs, err := Create(Memory, fname, 1, Byte, 256, 256)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	defer func() { _ = VSIUnlink(fname) }()
+	defer vrtDs.Close()
+
+	_, err = vrtDs.Dem("/vsimem/out.tiff", "color-relief", "", []string{})
+	assert.Error(t, err)
+}
+
+func TestDemColorReliefInvalidFilename(t *testing.T) {
+	fname := "/vsimem/test.tiff"
+	vrtDs, err := Create(Memory, fname, 1, Byte, 256, 256)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	defer func() { _ = VSIUnlink(fname) }()
+	defer vrtDs.Close()
+
+	invalidColorReliefFilename := "invalid_file"
+	_, err = vrtDs.Dem("/vsimem/out.tiff", "color-relief", invalidColorReliefFilename, []string{})
+	assert.Error(t, err)
+}