Add GDALGrid binding and Grid go function (#113)

errors.go

@@ -150,6 +150,7 @@ func ErrLogger(fn ErrorHandler) interface {
 	StatisticsOption
 	SetStatisticsOption
 	ClearStatisticsOption
+	GridOption
 } {
 	return errorCallback{fn}
 }
@@ -364,6 +365,9 @@ func (ec errorCallback) setClearStatisticsOpt(o *clearStatisticsOpt) {
 func (ec errorCallback) setGridCreateOpt(o *gridCreateOpts) {
 	o.errorHandler = ec.fn
 }
+func (ec errorCallback) setGridOpt(o *gridOpts) {
+	o.errorHandler = ec.fn
+}
 
 type multiError struct {
 	errs []error

godal.cpp

@@ -1725,3 +1725,25 @@ void godalGridCreate(cctx *ctx, char *pszAlgorithm, GDALGridAlgorithm eAlgorithm
 
 	godalUnwrap();
 }
+
+GDALDatasetH godalGrid(cctx *ctx, const char *pszDest, GDALDatasetH hSrcDS, char **switches) {
+	godalWrap(ctx);
+
+	GDALGridOptions *gridopts = GDALGridOptionsNew(switches,nullptr);
+	if(failed(ctx)) {
+		GDALGridOptionsFree(gridopts);
+		godalUnwrap();
+		return nullptr;
+	}
+
+
+	int usageErr=0;
+	GDALDatasetH ret = GDALGrid(pszDest, hSrcDS, gridopts, &usageErr);
+	GDALGridOptionsFree(gridopts);
+	if(ret==nullptr || usageErr!=0) {
+		forceError(ctx);
+	}
+
+	godalUnwrap();
+	return ret;
+}

godal.go

@@ -3814,6 +3814,41 @@ func GridCreate(pszAlgorithm string,
 	return nil
 }
 
+// Grid runs the library version of gdal_grid.
+// See the gdal_grid doc page to determine the valid flags/opts that can be set in switches.
+//
+// Example switches :
+//
+//	[]string{"-a", "maximum", "-txe", "0", "1"}
+//
+// Creation options and driver may be set in the switches slice with
+//
+//	switches:=[]string{"-co","TILED=YES","-of","GTiff"}
+//
+// NOTE: Some switches are NOT compatible with this binding, as a `nullptr` is passed to a later call to
+// `GDALGridOptionsNew()` (as the 2nd argument). Those switches are: "-oo", "-q", "-quiet"
+func (ds *Dataset) Grid(destPath string, switches []string, opts ...GridOption) (*Dataset, error) {
+	gridOpts := gridOpts{}
+	for _, opt := range opts {
+		opt.setGridOpt(&gridOpts)
+	}
+
+	cswitches := sliceToCStringArray(switches)
+	defer cswitches.free()
+
+	dest := unsafe.Pointer(C.CString(destPath))
+	cgc := createCGOContext(nil, gridOpts.errorHandler)
+	var dsRet C.GDALDatasetH
+	defer C.free(unsafe.Pointer(dest))
+
+	dsRet = C.godalGrid(cgc.cPointer(), (*C.char)(dest), ds.handle(), cswitches.cPointer())
+	if err := cgc.close(); err != nil {
+		return nil, err
+	}
+
+	return &Dataset{majorObject{C.GDALMajorObjectH(dsRet)}}, nil
+}
+
 type cgoContext struct {
 	cctx *C.cctx
 	opts cStringArray

godal.h

@@ -152,6 +152,7 @@ extern "C" {
     int godalGetRasterStatistics(cctx *ctx, GDALRasterBandH bnd, int bApproxOK, double *pdfMin, double *pdfMax, double *pdfMean, double *pdfStdDev);
     void godalSetRasterStatistics(cctx *ctx, GDALRasterBandH bnd, double dfMin, double dfMax, double dfMean, double dfStdDev);
 	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);
 #ifdef __cplusplus
 }
 #endif

godal_test.go

@@ -3979,6 +3979,95 @@ func TestStatistics(t *testing.T) {
 	assert.Error(t, err)
 }
 
+func TestGridLinear(t *testing.T) {
+	var (
+		err      error
+		outXSize = 256
+		outYSize = 256
+	)
+
+	vrtDs, err := CreateVector(Memory, "")
+	if err != nil {
+		t.Error(err)
+		return
+	}
+
+	geom, err := NewGeometryFromWKT("POLYGON((0 0 0, 0 1 1, 1 1 0, 1 0 1))", nil)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	_, 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 := 1
+	yMax := 0
+	argsString := fmt.Sprintf("-a linear -txe 0 1 -tye %d %d -outsize %d %d -ot Float64", yMin, yMax, outXSize, outYSize)
+	fname := "/vsimem/test.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
+		}
+	}
+	defer func() { _ = VSIUnlink(fname) }()
+	defer gridDs.Close()
+
+	var gridBindingPoints = make([]float64, outXSize*outYSize)
+	err = gridDs.Read(0, 0, gridBindingPoints, outXSize, outYSize)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+
+	var (
+		topLeftIndex     = 0
+		topRightIndex    = outXSize - 1
+		bottomLeftIndex  = outXSize * (outYSize - 1)
+		bottomRightIndex = (outXSize * outYSize) - 1
+		imageCentreIndex = outYSize*(outXSize/2) - 1
+	)
+
+	// For linear interpolation, we expect z-values of corners to match the input coordinates
+	// and the centre value to be the average of the 4 corner values
+	// 	TL (0, 0, 0), EXPECTED OUTPUT Z-VAL  = 0
+	// 	TR (1, 0, 1), EXPECTED OUTPUT Z-VAL  = 1
+	// 	BL (0, 1, 0), EXPECTED OUTPUT Z-VAL  = 0
+	// 	BR (1, 1, 1), EXPECTED OUTPUT Z-VAL  = 1
+	//  CR (0.5, 0.5), EXPECTED OUTPUT Z-VAL = (0 + 1 + 0 + 1) / 4 = 0.5
+	//
+	// NOTE: The input X and Y coords are offset slightly in GDAL before they're passed into a a gridding algorithm.
+	// This is why "TR" and "BL" below are expected to be 0.00390625 and NOT 0.
+	// See the `dfXPoint` and `dfYPoint` values in `GDALGridJobProcess()` for how these points are calculated
+	// TL
+	assert.Equal(t, 1.0, gridBindingPoints[topLeftIndex])
+	// TR
+	assert.Equal(t, 0.00390625, gridBindingPoints[topRightIndex])
+	// BL
+	assert.Equal(t, 0.00390625, gridBindingPoints[bottomLeftIndex])
+	// BR
+	assert.Equal(t, 1.0, gridBindingPoints[bottomRightIndex])
+	// Center
+	assert.Equal(t, 0.5, gridBindingPoints[imageCentreIndex])
+}
+
 func TestGridCreateLinear(t *testing.T) {
 	var (
 		err error
@@ -4035,6 +4124,75 @@ func TestGridCreateLinear(t *testing.T) {
 	assert.Equal(t, 0.5, gridCreateBindingPoints[imageCentreIndex])
 }
 
+func TestGridMaximum(t *testing.T) {
+	var (
+		err      error
+		outXSize = 256
+		outYSize = 256
+	)
+
+	vrtDs, err := CreateVector(Memory, "")
+	if err != nil {
+		t.Error(err)
+		return
+	}
+
+	geom, err := NewGeometryFromWKT("POLYGON((0 0 0, 0 1 1, 1 1 0, 1 0 1))", nil)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	_, 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
+	}
+
+	// NOTE: Flipping the arguments after `-tye` here, to account for `ProcessLayer` (in `GDALGrid`) flipping the coords "north up"
+	argsString := fmt.Sprintf("-a maximum -txe 0 1 -tye 0 1 -outsize %d %d -ot Float64", outXSize, outYSize)
+	fname := "/vsimem/test.tiff"
+
+	gridDs, err := vrtDs.Grid(fname, strings.Split(argsString, " "))
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	defer func() { _ = VSIUnlink(fname) }()
+	defer gridDs.Close()
+
+	var gridBindingPoints = make([]float64, outXSize*outYSize)
+	err = gridDs.Read(0, 0, gridBindingPoints, outXSize, outYSize)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+
+	var (
+		topLeftIndex     = 0
+		topRightIndex    = outXSize - 1
+		bottomLeftIndex  = outXSize * (outYSize - 1)
+		bottomRightIndex = (outXSize * outYSize) - 1
+		imageCentreIndex = outYSize*(outXSize/2) - 1
+	)
+
+	// All sampled values are expected to match the "maximum" value in the grid i.e. 1
+	// TL
+	assert.Equal(t, 1.0, gridBindingPoints[topLeftIndex])
+	// TR
+	assert.Equal(t, 1.0, gridBindingPoints[topRightIndex])
+	// BL
+	assert.Equal(t, 1.0, gridBindingPoints[bottomLeftIndex])
+	// BR
+	assert.Equal(t, 1.0, gridBindingPoints[bottomRightIndex])
+	// Center
+	assert.Equal(t, 1.0, gridBindingPoints[imageCentreIndex])
+}
+
 func TestGridCreateMaximum(t *testing.T) {
 	var (
 		err error
@@ -4072,3 +4230,32 @@ func TestGridCreateMaximum(t *testing.T) {
 	// Center
 	assert.Equal(t, 1.0, gridCreateBindingPoints[imageCentreIndex])
 }
+
+func TestGridInvalidSwitch(t *testing.T) {
+	vrtDs, err := CreateVector(Memory, "")
+	if err != nil {
+		t.Error(err)
+		return
+	}
+
+	geom, err := NewGeometryFromWKT("POLYGON((0 0 0, 0 1 1, 1 1 0, 1 0 1))", nil)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	_, 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
+	}
+	_, err = vrtDs.Grid("/vsimem/test.tiff", []string{"-invalidswitch"})
+	assert.Error(t, err)
+	ehc := eh()
+	_, err = vrtDs.Grid("/vsimem/test.tiff", []string{"-invalidswitch"}, ErrLogger(ehc.ErrorHandler))
+	assert.Error(t, err)
+}

options.go

@@ -1236,6 +1236,15 @@ type GridCreateOption interface {
 	setGridCreateOpt(mo *gridCreateOpts)
 }
 
+type gridOpts struct {
+	errorHandler ErrorHandler
+}
+
+// GridOption is an option that can be passed to Dataset.Grid()
+type GridOption interface {
+	setGridOpt(gOpt *gridOpts)
+}
+
 // RasterizeGeometryOption is an option that can be passed tp Dataset.RasterizeGeometry()
 type RasterizeGeometryOption interface {
 	setRasterizeGeometryOpt(o *rasterizeGeometryOpts)