Merge upstream gdal

.gitignore

@@ -0,0 +1,5 @@
+examples/tiff/tiff
+examples/translate/translate
+examples/warp/warp
+examples/grid/grid
+.idea

README

@@ -16,7 +16,7 @@ Installation
 Compatibility
 -------------
 
-This software has been tested most recently on Ubuntu 18.10, GDAL version 2.3.2.
+This software has been tested most recently on Debian 11 Bullseye, GDAL version 3.2.2.
 
 -------------
 Examples

algorithms.go

@@ -3,14 +3,10 @@ package gdal
 /*
 #include "go_gdal.h"
 #include "gdal_version.h"
-
-#cgo linux  pkg-config: gdal
-#cgo darwin pkg-config: gdal
-#cgo windows LDFLAGS: -Lc:/gdal/release-1600-x64/lib -lgdal_i
-#cgo windows CFLAGS: -IC:/gdal/release-1600-x64/include
 */
 import "C"
 import (
+	"errors"
 	"fmt"
 	"unsafe"
 )
@@ -302,5 +298,242 @@ func (src RasterBand) SieveFilter(
 /* Gridding functions                            */
 /* --------------------------------------------- */
 
-//Unimplemented: CreateGrid
+// GridAlgorithm represents Grid Algorithm code
+type GridAlgorithm int
+
+//
+const (
+	GA_InverseDistancetoAPower                = GridAlgorithm(C.GGA_InverseDistanceToAPower)
+	GA_MovingAverage                          = GridAlgorithm(C.GGA_MovingAverage)
+	GA_NearestNeighbor                        = GridAlgorithm(C.GGA_NearestNeighbor)
+	GA_MetricMinimum                          = GridAlgorithm(C.GGA_MetricMinimum)
+	GA_MetricMaximum                          = GridAlgorithm(C.GGA_MetricMaximum)
+	GA_MetricRange                            = GridAlgorithm(C.GGA_MetricRange)
+	GA_MetricCount                            = GridAlgorithm(C.GGA_MetricCount)
+	GA_MetricAverageDistance                  = GridAlgorithm(C.GGA_MetricAverageDistance)
+	GA_MetricAverageDistancePts               = GridAlgorithm(C.GGA_MetricAverageDistancePts)
+	GA_Linear                                 = GridAlgorithm(C.GGA_Linear)
+	GA_InverseDistanceToAPowerNearestNeighbor = GridAlgorithm(C.GGA_InverseDistanceToAPowerNearestNeighbor)
+)
+
+// GridLinearOptions: Linear method control options.
+type GridLinearOptions struct {
+	// Radius: in case the point to be interpolated does not fit into a triangle of the Delaunay triangulation,
+	// use that maximum distance to search a nearest neighbour, or use nodata otherwise. If set to -1, the search
+	// distance is infinite. If set to 0, nodata value will be always used.
+	Radius float64
+	// NoDataValue: no data marker to fill empty points.
+	NoDataValue float64
+}
+
+// GridInverseDistanceToAPowerOptions: Inverse distance to a power method control options.
+type GridInverseDistanceToAPowerOptions struct {
+	// Power: Weighting power
+	Power float64
+	// Smoothing: Smoothing parameter
+	Smoothing float64
+	// AnisotropyRatio: Reserved for future use
+	AnisotropyRatio float64
+	// AnisotropyAngle: Reserved for future use
+	AnisotropyAngle float64
+	// Radius1: The first radius (X axis if rotation angle is 0) of search ellipse.
+	Radius1 float64
+	// Radius2: The second radius (Y axis if rotation angle is 0) of search ellipse.
+	Radius2 float64
+	// Angle: Angle of ellipse rotation in degrees. Ellipse rotated counter clockwise.
+	Angle float64
+	// MaxPoints: Maximum number of data points to use.
+	// Do not search for more points than this number. If less amount of points found the grid node
+	// considered empty and will be filled with NODATA marker.
+	MaxPoints uint32
+	// MinPoints: Minimum number of data points to use.
+	// If less amount of points found the grid node considered empty and will be filled with NODATA marker.
+	MinPoints uint32
+	// NoDataValue: No data marker to fill empty points.
+	NoDataValue float64
+}
+
+// GridInverseDistanceToAPowerNearestNeighborOptions: Inverse distance to a power, with nearest neighbour search,
+// control options
+type GridInverseDistanceToAPowerNearestNeighborOptions struct {
+	// Power: Weighting power
+	Power float64
+	// Radius: The radius of search circle
+	Radius float64
+	// Smoothing: Smoothing parameter
+	Smoothing float64
+	// MaxPoints: Maximum number of data points to use.
+	// Do not search for more points than this number. If less amount of points found the grid node
+	// considered empty and will be filled with NODATA marker.
+	MaxPoints uint32
+	// MinPoints: Minimum number of data points to use.
+	// If less amount of points found the grid node considered empty and will be filled with NODATA marker.
+	MinPoints uint32
+	// NoDataValue: No data marker to fill empty points.
+	NoDataValue float64
+}
+
+// GridMovingAverageOptions: Moving average method control options
+type GridMovingAverageOptions struct {
+	// Radius1: The first radius (X axis if rotation angle is 0) of search ellipse.
+	Radius1 float64
+	// Radius2: The second radius (Y axis if rotation angle is 0) of search ellipse.
+	Radius2 float64
+	// Angle: Angle of ellipse rotation in degrees. Ellipse rotated counter clockwise.
+	Angle float64
+	// MinPoints: Minimum number of data points to use.
+	// If less amount of points found the grid node considered empty and will be filled with NODATA marker.
+	MinPoints uint32
+	// NoDataValue: No data marker to fill empty points.
+	NoDataValue float64
+}
+
+// GridNearestNeighborOptions: Nearest neighbor method control options.
+type GridNearestNeighborOptions struct {
+	// Radius1: The first radius (X axis if rotation angle is 0) of search ellipse.
+	Radius1 float64
+	// Radius2: The second radius (Y axis if rotation angle is 0) of search ellipse.
+	Radius2 float64
+	// Angle: Angle of ellipse rotation in degrees. Ellipse rotated counter clockwise.
+	Angle float64
+	// NoDataValue: No data marker to fill empty points.
+	NoDataValue float64
+}
+
+// GridDataMetricsOptions: Data metrics method control options
+type GridDataMetricsOptions struct {
+	// Radius1: The first radius (X axis if rotation angle is 0) of search ellipse.
+	Radius1 float64
+	// Radius2: The second radius (Y axis if rotation angle is 0) of search ellipse.
+	Radius2 float64
+	// Angle: Angle of ellipse rotation in degrees. Ellipse rotated counter clockwise.
+	Angle float64
+	// MinPoints: Minimum number of data points to use.
+	// If less amount of points found the grid node considered empty and will be filled with NODATA marker.
+	MinPoints uint32
+	// NoDataValue: No data marker to fill empty points.
+	NoDataValue float64
+}
+
+var errInvalidOptionsTypeWasPassed = errors.New("invalid options type was passed")
+
+// GridCreate: Create regular grid from the scattered data.
+// This function takes the arrays of X and Y coordinates and corresponding Z values as input and computes
+// regular grid (or call it a raster) from these scattered data. You should supply geometry and extent of the
+// output grid.
+func GridCreate(
+	algorithm GridAlgorithm,
+	options interface{},
+	x, y, z []float64,
+	xMin, xMax, yMin, yMax float64,
+	nX, nY uint,
+	progress ProgressFunc,
+	data interface{},
+) ([]float64, error) {
+	if len(x) != len(y) || len(x) != len(z) {
+		return nil, errors.New("lengths of x, y, z should equal")
+	}
+
+	poptions := unsafe.Pointer(nil)
+	switch algorithm {
+	case GA_InverseDistancetoAPower:
+		soptions, ok := options.(GridInverseDistanceToAPowerOptions)
+		if !ok {
+			return nil, errInvalidOptionsTypeWasPassed
+		}
+		poptions = unsafe.Pointer(&C.GDALGridInverseDistanceToAPowerOptions{
+			dfPower:           C.double(soptions.Power),
+			dfSmoothing:       C.double(soptions.Smoothing),
+			dfAnisotropyRatio: C.double(soptions.AnisotropyRatio),
+			dfAnisotropyAngle: C.double(soptions.AnisotropyAngle),
+			dfRadius1:         C.double(soptions.Radius1),
+			dfRadius2:         C.double(soptions.Radius2),
+			dfAngle:           C.double(soptions.Angle),
+			nMaxPoints:        C.uint(soptions.MaxPoints),
+			nMinPoints:        C.uint(soptions.MinPoints),
+			dfNoDataValue:     C.double(soptions.NoDataValue),
+		})
+	case GA_InverseDistanceToAPowerNearestNeighbor:
+		soptions, ok := options.(GridInverseDistanceToAPowerNearestNeighborOptions)
+		if !ok {
+			return nil, errInvalidOptionsTypeWasPassed
+		}
+		poptions = unsafe.Pointer(&C.GDALGridInverseDistanceToAPowerNearestNeighborOptions{
+			dfPower:       C.double(soptions.Power),
+			dfRadius:      C.double(soptions.Radius),
+			dfSmoothing:   C.double(soptions.Smoothing),
+			nMaxPoints:    C.uint(soptions.MaxPoints),
+			nMinPoints:    C.uint(soptions.MinPoints),
+			dfNoDataValue: C.double(soptions.NoDataValue),
+		})
+	case GA_MovingAverage:
+		soptions, ok := options.(GridMovingAverageOptions)
+		if !ok {
+			return nil, errInvalidOptionsTypeWasPassed
+		}
+		poptions = unsafe.Pointer(&C.GDALGridMovingAverageOptions{
+			dfRadius1:     C.double(soptions.Radius1),
+			dfRadius2:     C.double(soptions.Radius2),
+			dfAngle:       C.double(soptions.Angle),
+			nMinPoints:    C.uint(soptions.MinPoints),
+			dfNoDataValue: C.double(soptions.NoDataValue),
+		})
+	case GA_NearestNeighbor:
+		soptions, ok := options.(GridNearestNeighborOptions)
+		if !ok {
+			return nil, errInvalidOptionsTypeWasPassed
+		}
+		poptions = unsafe.Pointer(&C.GDALGridNearestNeighborOptions{
+			dfRadius1:     C.double(soptions.Radius1),
+			dfRadius2:     C.double(soptions.Radius2),
+			dfAngle:       C.double(soptions.Angle),
+			dfNoDataValue: C.double(soptions.NoDataValue),
+		})
+	case GA_MetricMinimum, GA_MetricMaximum, GA_MetricCount, GA_MetricRange,
+		GA_MetricAverageDistance, GA_MetricAverageDistancePts:
+		soptions, ok := options.(GridDataMetricsOptions)
+		if !ok {
+			return nil, errInvalidOptionsTypeWasPassed
+		}
+		poptions = unsafe.Pointer(&C.GDALGridDataMetricsOptions{
+			dfRadius1:     C.double(soptions.Radius1),
+			dfRadius2:     C.double(soptions.Radius2),
+			dfAngle:       C.double(soptions.Angle),
+			nMinPoints:    C.uint(soptions.MinPoints),
+			dfNoDataValue: C.double(soptions.NoDataValue),
+		})
+	case GA_Linear:
+		soptions, ok := options.(GridLinearOptions)
+		if !ok {
+			return nil, errInvalidOptionsTypeWasPassed
+		}
+		poptions = unsafe.Pointer(&C.GDALGridLinearOptions{
+			dfRadius:      C.double(soptions.Radius),
+			dfNoDataValue: C.double(soptions.NoDataValue),
+		})
+	}
+
+	buffer := make([]float64, nX*nY)
+	arg := &goGDALProgressFuncProxyArgs{progress, data}
+	err := C.GDALGridCreate(
+		C.GDALGridAlgorithm(algorithm),
+		poptions,
+		C.uint(uint(len(x))),
+		(*C.double)(unsafe.Pointer(&x[0])),
+		(*C.double)(unsafe.Pointer(&y[0])),
+		(*C.double)(unsafe.Pointer(&z[0])),
+		C.double(xMin),
+		C.double(xMax),
+		C.double(yMin),
+		C.double(yMax),
+		C.uint(nX),
+		C.uint(nY),
+		C.GDALDataType(Float64),
+		unsafe.Pointer(&buffer[0]),
+		C.goGDALProgressFuncProxyB(),
+		unsafe.Pointer(arg),
+	).Err()
+	return buffer, err
+}
+
 //Unimplemented: ComputeMatchingPoints

algorithms_test.go

@@ -0,0 +1,80 @@
+package gdal
+
+import (
+	"errors"
+	"fmt"
+	"io/ioutil"
+	"math"
+	"strconv"
+	"strings"
+	"testing"
+)
+
+func readGridFile(filename string) (x, y, z []float64, err error) {
+	var b []byte
+	if b, err = ioutil.ReadFile(filename); err != nil {
+		return
+	}
+	arr := strings.Split(string(b), "\n")
+	x, y, z = make([]float64, len(arr)), make([]float64, len(arr)), make([]float64, len(arr))
+	for i, el := range arr {
+		xyz := strings.Split(el, ",")
+		if len(xyz) != 3 {
+			err = errors.New("wrong input file format, should be CSV with 3 columns: y,x,z")
+		}
+
+		if y[i], err = strconv.ParseFloat(xyz[0], 64); err != nil {
+			return
+		}
+		if x[i], err = strconv.ParseFloat(xyz[1], 64); err != nil {
+			return
+		}
+		if z[i], err = strconv.ParseFloat(xyz[2], 64); err != nil {
+			return
+		}
+	}
+	return
+}
+
+func TestGridCreate(t *testing.T) {
+	x, y, z, err := readGridFile("testdata/grid.csv")
+	if err != nil {
+		t.Fatalf("failed to readGridFile: %v", err)
+		return
+	}
+	var nX, nY uint = 420, 470
+
+	// finding max and min values
+	var xMin, xMax, yMin, yMax = math.MaxFloat64, -math.MaxFloat64, math.MaxFloat64, -math.MaxFloat64
+	for i := range x {
+		if x[i] < xMin {
+			xMin = x[i]
+		}
+		if x[i] > xMax {
+			xMax = x[i]
+		}
+		if y[i] < yMin {
+			yMin = y[i]
+		}
+		if y[i] > yMax {
+			yMax = y[i]
+		}
+	}
+
+	fmt.Println("Calling gdal.GridCreate")
+	data, err := GridCreate(
+		GA_Linear,
+		GridLinearOptions{Radius: -1, NoDataValue: 0},
+		x, y, z,
+		xMin, xMax, yMin, yMax,
+		nX, nY,
+		DummyProgress,
+		nil,
+	)
+	if err != nil {
+		t.Errorf("GridCreate: %v", err)
+	}
+	if expectedDataLen := int(nX * nY); len(data) != expectedDataLen {
+		t.Errorf("expected length of data equal to %d", expectedDataLen)
+	}
+}