Merge pull request lukeroth#51 from mtfelian/master

Implemented GDALGridCreate CGo wrapper for all available algorithms

.gitignore

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

algorithms.go

@@ -11,6 +11,7 @@ package gdal
 */
 import "C"
 import (
+	"errors"
 	"fmt"
 	"unsafe"
 )
@@ -302,5 +303,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)
+	}
+}

doc.go

@@ -24,7 +24,7 @@ A simple program to create a georeferenced blank 256x256 GeoTIFF:
 		flag.Parse()
 		filename := flag.Arg(0)
 		if filename == "" {
-			fmt.Printf("Usage: test_tiff [filename]\n")
+			fmt.Printf("Usage: tiff [filename]\n")
 			return
 		}
 		buffer := make([]uint8, 256 * 256)

examples/grid/go.mod

@@ -0,0 +1,7 @@
+module github.com/lukeroth/gdal/examples/grid
+
+go 1.13
+
+replace github.com/lukeroth/gdal => ../..
+
+require github.com/lukeroth/gdal v0.0.0-00010101000000-000000000000