Add GDALViewshedGenerate binding (#142)

driver.go

@@ -124,6 +124,10 @@ func (dn DriverName) setRasterizeOpt(to *rasterizeOpts) {
 	to.driver = dn
 }
 
+func (dn DriverName) setViewshedOpt(to *viewshedOpts) {
+	to.driver = dn
+}
+
 type driversOpt struct {
 	drivers []string
 }

errors.go

@@ -154,6 +154,7 @@ func ErrLogger(fn ErrorHandler) interface {
 	GridOption
 	NearblackOption
 	DemOption
+	ViewshedOption
 	SetGCPsOption
 	GCPsToGeoTransformOption
 	RegisterPluginOption
@@ -388,6 +389,9 @@ func (ec errorCallback) setNearblackOpt(o *nearBlackOpts) {
 func (ec errorCallback) setDemOpt(o *demOpts) {
 	o.errorHandler = ec.fn
 }
+func (ec errorCallback) setViewshedOpt(o *viewshedOpts) {
+	o.errorHandler = ec.fn
+}
 func (ec errorCallback) setSetGCPsOpt(o *setGCPsOpts) {
 	o.errorHandler = ec.fn
 }

godal.cpp

@@ -1905,6 +1905,24 @@ GDALDatasetH godalDem(cctx *ctx, const char *pszDest, const char *pszProcessing,
 	return ret;
 }
 
+GDALDatasetH godalViewshedGenerate(cctx *ctx, GDALRasterBandH bnd, const char *pszDriverName, const char *pszTargetRasterName, const char **papszCreationOptions, double dfObserverX,
+    					   double dfObserverY, double dfObserverHeight, double dfTargetHeight, double dfVisibleVal, double dfInvisibleVal, double dfOutOfRangeVal, 
+						   double dfNoDataVal, double dfCurvCoeff, GUInt32 eMode, double dfMaxDistance, GUInt32 heightMode) {
+	godalWrap(ctx);
+	GDALDatasetH ret = nullptr;
+#if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 1, 0)
+	ret = GDALViewshedGenerate(bnd, pszDriverName, pszTargetRasterName, papszCreationOptions, dfObserverX, dfObserverY, dfObserverHeight, dfTargetHeight, dfVisibleVal, dfInvisibleVal, dfOutOfRangeVal, 
+	dfNoDataVal, dfCurvCoeff, GDALViewshedMode(eMode), dfMaxDistance, nullptr, nullptr, GDALViewshedOutputType(heightMode), nullptr);
+#else
+	CPLError(CE_Failure, CPLE_AppDefined, "Viewshed not implemented in gdal < 3.1");
+#endif
+	if(ret == nullptr) {
+		forceError(ctx);
+	}
+	godalUnwrap();
+	return ret;
+}
+
 OGRSpatialReferenceH godalGetGCPSpatialRef(GDALDatasetH hSrcDS) {
 	return GDALGetGCPSpatialRef(hSrcDS);
 }

godal.go

@@ -4060,6 +4060,96 @@ func (ds *Dataset) Dem(destPath, processingMode string, colorFilename string, sw
 	return &Dataset{majorObject{C.GDALMajorObjectH(dsRet)}}, nil
 }
 
+// ViewshedMode is the "cell height calculation mode" for the viewshed process
+//
+// Source: https://github.com/OSGeo/gdal/blob/master/alg/viewshed/viewshed_types.h
+type ViewshedMode uint32
+
+const (
+	// MDiagonal is the "diagonal mode"
+	MDiagonal ViewshedMode = iota + 1
+	// MEdge is the "edge mode"
+	MEdge
+	// MMax is the "maximum value produced by Diagonal and Edge mode"
+	MMax
+	// MMin is the "minimum value produced by Diagonal and Edge mode"
+	MMin
+)
+
+// ViewshedOutputType sets the return type and information represented by the returned data
+//
+// Source: https://gdal.org/en/stable/programs/gdal_viewshed.html
+//
+// NOTE: "Cumulative (ACCUM)" mode not currently supported, as it's not available in the `GDALViewshedGenerate` function
+// (it's only used in the command line invocation of `viewshed`)
+type ViewshedOutputType uint32
+
+const (
+	// Normal returns a raster of type Byte containing visible locations
+	Normal ViewshedOutputType = iota + 1
+	// MinTargetHeightFromDem return a raster of type Float64 containing the minimum target height for target to be visible from the DEM surface
+	MinTargetHeightFromDem
+	// MinTargetHeightFromGround return a raster of type Float64 containing the minimum target height for target to be visible from ground level
+	MinTargetHeightFromGround
+)
+
+// Viewshed (binding for GDALViewshedGenerate), creates a viewshed from a raster DEM, these parameters (mostly) map to parameters for GDALViewshedGenerate
+// for more information see: https://gdal.org/en/stable/api/gdal_alg.html#_CPPv420GDALViewshedGenerate15GDALRasterBandHPKcPKc12CSLConstListddddddddd16GDALViewshedModed16GDALProgressFuncPv22GDALViewshedOutputType12CSLConstList
+//
+// Several Viewshed parameters have default values defined in GDAL, see `Options` in: https://github.com/OSGeo/gdal/blob/master/alg/viewshed/viewshed_types.h
+//
+// to define different values for these parameters, provide the corresponding options:
+//
+//   - DriverName(dn DriverName) sets the GDAL driver
+//   - TargetHeight(h float64) sets the target height above the DEM surface
+//   - VisibilityVals(vv float64, iv float64) sets the raster output value for visible and invisible pixels
+//   - OutOfRangeVal(ov float64) sets the raster output value for pixels outside of max distance
+//   - MaxDistance(d float64) sets the maximum number of pixels to search from the observer
+//   - CurveCoeff(cc float64) sets the coefficient for atmospheric refraction
+//   - NoDataVal(nd float64) sets the raster output value for pixels with no data
+//   - CellMode(cm ViewshedMode) sets mode of cell height calculation
+//   - HeightMode(hm ViewshedOutputType) sets the type of the output information
+//   - CreationOption(opts ...string) options to pass to a driver when creating a dataset
+//
+// WARNING: One of the Godal tests for this function is ported from a GDAL test finalised in this commit (tagged for 3.10.2RC1): https://github.com/OSGeo/gdal/commit/33dd00c63155250afce04092c77cb225570efa64
+//
+// Automated testing across different GDAL versions shows that GDALViewshedGenerate FAILS some of these tests on versions of GDAL < 3.10.0 as follows:
+//
+//   - version < 3.1.0: viewshed is not supported and will always throw a "not implemented" error
+//   - 3.1.0 <= version < 3.4.2: all tests fail
+//   - 3.4.2 <= version < 3.10.0: tests where heightMode == MinTargetHeightFromDem fail
+//   - version >= 3.10.0: all tests pass
+func (srcBand Band) Viewshed(targetRasterName string, observerX float64, observerY float64, observerHeight float64, opts ...ViewshedOption) (*Dataset, error) {
+	vso := viewshedOpts{
+		driver:     GTiff,
+		visibleVal: 255,
+		noDataVal:  -1,
+		curveCoeff: .85714,
+		cellMode:   MEdge,
+		heightMode: Normal,
+	}
+	for _, opt := range opts {
+		opt.setViewshedOpt(&vso)
+	}
+
+	copts := sliceToCStringArray(vso.creation)
+	defer copts.free()
+	driver := unsafe.Pointer(C.CString(string(vso.driver)))
+	defer C.free(unsafe.Pointer(driver))
+	targetRaster := unsafe.Pointer(C.CString(targetRasterName))
+	defer C.free(unsafe.Pointer(targetRaster))
+
+	cgc := createCGOContext(nil, vso.errorHandler)
+	dsRet := C.godalViewshedGenerate(cgc.cPointer(), srcBand.handle(), (*C.char)(driver), (*C.char)(targetRaster), copts.cPointer(), C.double(observerX),
+		C.double(observerY), C.double(observerHeight), C.double(vso.targetHeight), C.double(vso.visibleVal), C.double(vso.invisibleVal), C.double(vso.outOfRangeVal),
+		C.double(vso.noDataVal), C.double(vso.curveCoeff), C.uint(vso.cellMode), C.double(vso.maxDistance), C.uint(vso.heightMode))
+	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

@@ -171,6 +171,10 @@ extern "C" {
 	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);
+	GDALDatasetH godalViewshedGenerate(cctx *ctx, GDALRasterBandH bnd, const char *pszDriverName, const char *pszTargetRasterName, const char **papszCreationOptions,
+							   double dfObserverX, double dfObserverY, double dfObserverHeight, double dfTargetHeight, double dfVisibleVal, 
+							   double dfInvisibleVal, double dfOutOfRangeVal,  double dfNoDataVal, double dfCurvCoeff, GUInt32 eMode, 
+							   double dfMaxDistance, GUInt32 heightMode);
 
 	typedef struct {
 		const GDAL_GCP *gcpList;

godal_test.go

@@ -4440,6 +4440,170 @@ func TestGridInvalidSwitch(t *testing.T) {
 	assert.Error(t, err)
 }
 
+// Test Ported from: https://github.com/OSGeo/gdal/blob/6cdae8b8f7d09ecf67e24959e984d2e7bbe3ee62/autotest/cpp/test_viewshed.cpp#L98
+func TestViewshedSimpleHeight(t *testing.T) {
+	// setup common to all scopes
+	var (
+		ehc = eh()
+		drv = DriverName("MEM")
+
+		identity = [6]float64{0, 1, 0, 0, 0, 1}
+
+		xLen = 5
+		yLen = 5
+		in   = []int8{
+			-1, 0, 1, 0, -1,
+			-1, 2, 0, 4, -1,
+			-1, 1, 0, -1, -1,
+			0, 3, 0, 2, 0,
+			-1, 0, 0, 3, -1,
+		}
+		observable = []float64{
+			4, 2, 0, 4, 8,
+			3, 2, 0, 4, 3,
+			2, 1, 0, -1, -2,
+			4, 3, 0, 2, 1,
+			6, 3, 0, 2, 4,
+		}
+	)
+	vrtDs, err := Create(Memory, "", 1, Int8, xLen, yLen)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	defer vrtDs.Close()
+	err = vrtDs.SetGeoTransform(identity)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	err = vrtDs.Bands()[0].IO(IOWrite, 0, 0, in, xLen, yLen)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+
+	// from cpp scope 1: normal
+	// NOTE: This test fails in releases older than 3.4.2
+	if CheckMinVersion(3, 4, 2) {
+		rds, err := vrtDs.Bands()[0].Viewshed("none", 2, 2, 0, DriverName(drv), CurveCoeff(0), ErrLogger(ehc.ErrorHandler))
+		if err != nil {
+			t.Error(err)
+			return
+		}
+		defer rds.Close()
+
+		out := make([]int8, xLen*yLen)
+		err = rds.Bands()[0].IO(IORead, 0, 0, out, xLen, yLen)
+		if err != nil {
+			t.Error(err)
+			return
+		}
+
+		expected := make([]int8, xLen*yLen)
+		for i := 0; i < len(in); i++ {
+			var v int8 = 0
+			if in[i] >= int8(observable[i]) {
+				v = 127
+			}
+			expected[i] = v
+		}
+		assert.Equal(t, expected, out)
+	}
+
+	// from cpp scope 2: dem
+	// NOTE: This test fails in releases older than 3.10
+	if CheckMinVersion(3, 10, 0) {
+		rds, err := vrtDs.Bands()[0].Viewshed("none", 2, 2, 0, DriverName(drv), CurveCoeff(0), HeightMode(MinTargetHeightFromDem))
+		if err != nil {
+			t.Error(err)
+			return
+		}
+		defer rds.Close()
+
+		dem := make([]float64, xLen*yLen)
+		err = rds.Bands()[0].IO(IORead, 0, 0, dem, xLen, yLen)
+		if err != nil {
+			t.Error(err)
+			return
+		}
+
+		expected := make([]float64, xLen*yLen)
+		copy(expected, observable)
+		for i := 0; i < len(expected); i++ {
+			expected[i] = math.Max(0.0, expected[i])
+		}
+		assert.Equal(t, expected, dem)
+	}
+
+	// from cpp scope 3: ground
+	// NOTE: This test fails in releases older than 3.4.2
+	if CheckMinVersion(3, 4, 2) {
+		rds, err := vrtDs.Bands()[0].Viewshed("none", 2, 2, 0, DriverName(drv), CurveCoeff(0), HeightMode(MinTargetHeightFromGround))
+		if err != nil {
+			t.Error(err)
+			return
+		}
+		defer rds.Close()
+
+		ground := make([]float64, xLen*yLen)
+		err = rds.Bands()[0].IO(IORead, 0, 0, ground, xLen, yLen)
+		if err != nil {
+			t.Error(err)
+			return
+		}
+
+		expected := make([]float64, xLen*yLen)
+		for i := 0; i < len(expected); i++ {
+			expected[i] = math.Max(0.0, observable[i]-float64(in[i]))
+		}
+		assert.Equal(t, expected, ground)
+	}
+}
+
+func TestViewshedCreationOptions(t *testing.T) {
+
+	var (
+		driver  = GTiff
+		tmpname = tempfile()
+	)
+	defer os.Remove(tmpname)
+	vrtDs, err := Create(driver, tmpname, 1, Int8, 20, 20, CreationOption("TILED=YES", "BLOCKXSIZE=128", "BLOCKYSIZE=128"))
+	if err != nil {
+		t.Error(err)
+		return
+	}
+	defer vrtDs.Close()
+	identity := [6]float64{0, 1, 0, 0, 0, 1}
+	err = vrtDs.SetGeoTransform(identity)
+	if err != nil {
+		t.Error(err)
+		return
+	}
+
+	// Invalid - with error logger
+	ehc := eh()
+	ds, err := vrtDs.Bands()[0].Viewshed("none", 2, 2, 0, CreationOption("INVALID_OPT=BAR"), ErrLogger(ehc.ErrorHandler))
+	if err == nil {
+		ds.Close()
+	}
+	assert.Error(t, err)
+
+	// Invalid - no error logger
+	ds, err = vrtDs.Bands()[0].Viewshed("none", 2, 2, 0, CreationOption("INVALID_OPT=BAR"))
+	if err == nil {
+		ds.Close()
+	}
+	assert.Error(t, err)
+
+	// Valid
+	ds, err = vrtDs.Bands()[0].Viewshed("none", 2, 2, 0, CreationOption("TILED=YES", "BLOCKXSIZE=128", "BLOCKYSIZE=128"))
+	if err == nil {
+		ds.Close()
+	}
+	assert.NoError(t, err)
+}
+
 func TestNearblackBlack(t *testing.T) {
 	// 1. Create an image, linearly interpolated, from black (on the left) to white (on the right), using `Grid()`
 	var (