Skip to content

Commit

Permalink
Add Grid bindings and GridCreate go function (#111)
Browse files Browse the repository at this point in the history
  • Loading branch information
pericles-tpt authored Sep 8, 2023
1 parent 71f3ed8 commit 5f05eb4
Show file tree
Hide file tree
Showing 8 changed files with 1,555 additions and 41 deletions.
3 changes: 3 additions & 0 deletions errors.go
Original file line number Diff line number Diff line change
Expand Up @@ -361,6 +361,9 @@ func (ec errorCallback) setSetStatisticsOpt(o *setStatisticsOpt) {
func (ec errorCallback) setClearStatisticsOpt(o *clearStatisticsOpt) {
o.errorHandler = ec.fn
}
func (ec errorCallback) setGridCreateOpt(o *gridCreateOpts) {
o.errorHandler = ec.fn
}

type multiError struct {
errs []error
Expand Down
18 changes: 14 additions & 4 deletions go.mod
Original file line number Diff line number Diff line change
Expand Up @@ -3,10 +3,20 @@ module github.com/airbusgeo/godal
go 1.15

require (
cloud.google.com/go/storage v1.17.0
cloud.google.com/go/iam v1.1.2 // indirect
cloud.google.com/go/storage v1.32.0
github.com/airbusgeo/cogger v0.0.7
github.com/airbusgeo/osio v0.1.2
github.com/airbusgeo/errs v0.0.3 // indirect
github.com/airbusgeo/osio v0.1.3
github.com/google/s2a-go v0.1.7 // indirect
github.com/google/uuid v1.3.1 // indirect
github.com/hashicorp/golang-lru v1.0.2 // indirect
github.com/spf13/cobra v1.2.1
github.com/stretchr/testify v1.7.0
google.golang.org/api v0.58.0
github.com/stretchr/testify v1.8.3
golang.org/x/oauth2 v0.12.0 // indirect
google.golang.org/api v0.138.0
google.golang.org/appengine v1.6.8 // indirect
google.golang.org/genproto v0.0.0-20230822172742-b8732ec3820d // indirect
google.golang.org/genproto/googleapis/api v0.0.0-20230822172742-b8732ec3820d // indirect
google.golang.org/grpc v1.58.0 // indirect
)
1,353 changes: 1,316 additions & 37 deletions go.sum

Large diffs are not rendered by default.

33 changes: 33 additions & 0 deletions godal.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@

#include <gdal_utils.h>
#include <gdal_alg.h>
#include <gdalgrid.h>

extern "C" {
extern long long int _gogdalSizeCallback(char* key, char** errorString);
Expand Down Expand Up @@ -1692,3 +1693,35 @@ void test_godal_error_handling(cctx *ctx) {
godalUnwrap();
}

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) {
godalWrap(ctx);
CPLErr ret;

if (!GDALHasTriangulation() && eAlgorithm == GGA_Linear) {
CPLError(CE_Failure, CPLE_AppDefined, "unable to run GGA_Linear algorithm, since GDAL built without QHull support");
godalUnwrap();
return;
}

void *ppOptions;
#if GDAL_VERSION_NUM >= GDAL_COMPUTE_VERSION(3, 7, 0)
ret = GDALGridParseAlgorithmAndOptions(pszAlgorithm, &eAlgorithm, &ppOptions);
#else
ret = ParseAlgorithmAndOptions(pszAlgorithm, &eAlgorithm, &ppOptions);
#endif
if(ret!=0) {
forceCPLError(ctx, ret);
godalUnwrap();
return;
}

ret = GDALGridCreate(eAlgorithm, ppOptions, nPoints, padfX, padfY,
padfZ, dfXMin, dfXMax, dfYMin, dfYMax, nXSize, nYSize,
eType, pData, nullptr, nullptr);
if(ret!=0) {
forceCPLError(ctx, ret);
}

godalUnwrap();
}
81 changes: 81 additions & 0 deletions godal.go
Original file line number Diff line number Diff line change
Expand Up @@ -1710,6 +1710,35 @@ func (ra ResamplingAlg) rioAlg() (C.GDALRIOResampleAlg, error) {
}
}

func gridAlgFromString(str string) (C.GDALGridAlgorithm, error) {
switch str {
case "invdist":
return C.GGA_InverseDistanceToAPower, nil
case "average":
return C.GGA_MovingAverage, nil
case "nearest":
return C.GGA_NearestNeighbor, nil
case "minimum":
return C.GGA_MetricMinimum, nil
case "maximum":
return C.GGA_MetricMaximum, nil
case "range":
return C.GGA_MetricRange, nil
case "count":
return C.GGA_MetricCount, nil
case "average_distance":
return C.GGA_MetricAverageDistance, nil
case "average_distance_pts":
return C.GGA_MetricAverageDistancePts, nil
case "linear":
return C.GGA_Linear, nil
case "invdistnn":
return C.GGA_InverseDistanceToAPowerNearestNeighbor, nil
default:
return C.GGA_InverseDistanceToAPower, fmt.Errorf("unknown gridding algorithm %s", str)
}
}

func bufferType(buffer interface{}) DataType {
switch buffer.(type) {
case []byte:
Expand Down Expand Up @@ -3733,6 +3762,58 @@ func BuildVRT(dstVRTName string, sourceDatasets []string, switches []string, opt
return &Dataset{majorObject{C.GDALMajorObjectH(hndl)}}, nil
}

// GridCreate, creates a grid from scattered data, given provided gridding parameters as a string (pszAlgorithm)
// and the arguments required for `godalGridCreate()` (binding for GDALGridCreate)
//
// NOTE: For valid gridding algorithm strings see: https://gdal.org/programs/gdal_grid.html#interpolation-algorithms
func GridCreate(pszAlgorithm string,
xCoords []float64,
yCoords []float64,
zCoords []float64,
dfXMin float64,
dfXMax float64,
dfYMin float64,
dfYMax float64,
nXSize int,
nYSize int,
buffer interface{},
opts ...GridCreateOption,
) error {
if len(xCoords) != len(yCoords) || len(yCoords) != len(zCoords) {
return errors.New("`xCoords`, `yCoords` and `zCoords` are not all equal length")
}

gco := gridCreateOpts{}
for _, o := range opts {
o.setGridCreateOpt(&gco)
}

griddingAlgStr := strings.Split(pszAlgorithm, ":")[0]
algCEnum, err := gridAlgFromString(griddingAlgStr)
if err != nil {
return err
}

var (
params = unsafe.Pointer(C.CString(pszAlgorithm))
cgc = createCGOContext(nil, gco.errorHandler)
)
defer C.free(params)

var (
dtype = bufferType(buffer)
dsize = dtype.Size()
numGridBytes = C.int(nXSize * nYSize * dsize)
cBuf = cBuffer(buffer, int(numGridBytes)/dsize)
)
cgc = createCGOContext(nil, gco.errorHandler)
C.godalGridCreate(cgc.cPointer(), (*C.char)(params), algCEnum, C.uint(len(xCoords)), cDoubleArray(xCoords), cDoubleArray(yCoords), cDoubleArray(zCoords), C.double(dfXMin), C.double(dfXMax), C.double(dfYMin), C.double(dfYMax), C.uint(nXSize), C.uint(nYSize), C.GDALDataType(dtype), cBuf)
if err := cgc.close(); err != nil {
return err
}
return nil
}

type cgoContext struct {
cctx *C.cctx
opts cStringArray
Expand Down
2 changes: 2 additions & 0 deletions godal.h
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@

#define _GNU_SOURCE 1
#include <gdal.h>
#include <gdal_alg.h>
#include <ogr_srs_api.h>
#include <cpl_conv.h>
#include "cpl_port.h"
Expand Down Expand Up @@ -150,6 +151,7 @@ extern "C" {
void godalComputeRasterStatistics(cctx *ctx, GDALRasterBandH bnd, int bApproxOK, double *pdfMin, double *pdfMax, double *pdfMean, double *pdfStdDev);
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);
#ifdef __cplusplus
}
#endif
Expand Down
95 changes: 95 additions & 0 deletions godal_test.go
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ import (
"os"
"path"
"path/filepath"
"strings"
"syscall"
"testing"
"time"
Expand Down Expand Up @@ -3977,3 +3978,97 @@ func TestStatistics(t *testing.T) {
_, _, err = bnd.GetStatistics()
assert.Error(t, err)
}

func TestGridCreateLinear(t *testing.T) {
var (
err error

xCoords = []float64{0, 1, 0, 1}
yCoords = []float64{0, 0, 1, 1}
zCoords = []float64{1, 0, 0, 1}
outXSize = 256
outYSize = 256
)

var gridCreateBindingPoints = make([]float64, outXSize*outYSize)
err = GridCreate("linear", xCoords, yCoords, zCoords, 0, 1, 0, 1, outXSize, outYSize, gridCreateBindingPoints)
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
}
}

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, gridCreateBindingPoints[topLeftIndex])
// TR
assert.Equal(t, 0.00390625, gridCreateBindingPoints[topRightIndex])
// BL
assert.Equal(t, 0.00390625, gridCreateBindingPoints[bottomLeftIndex])
// BR
assert.Equal(t, 1.0, gridCreateBindingPoints[bottomRightIndex])
// Center
assert.Equal(t, 0.5, gridCreateBindingPoints[imageCentreIndex])
}

func TestGridCreateMaximum(t *testing.T) {
var (
err error

xCoords = []float64{0, 1, 0, 1}
yCoords = []float64{0, 0, 1, 1}
zCoords = []float64{1, 0, 0, 1}
outXSize = 256
outYSize = 256
)

var gridCreateBindingPoints = make([]float64, outXSize*outYSize)
err = GridCreate("maximum", xCoords, yCoords, zCoords, 0, 1, 0, 1, outXSize, outYSize, gridCreateBindingPoints)
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, gridCreateBindingPoints[topLeftIndex])
// TR
assert.Equal(t, 1.0, gridCreateBindingPoints[topRightIndex])
// BL
assert.Equal(t, 1.0, gridCreateBindingPoints[bottomLeftIndex])
// BR
assert.Equal(t, 1.0, gridCreateBindingPoints[bottomRightIndex])
// Center
assert.Equal(t, 1.0, gridCreateBindingPoints[imageCentreIndex])
}
11 changes: 11 additions & 0 deletions options.go
Original file line number Diff line number Diff line change
Expand Up @@ -1225,6 +1225,17 @@ type rasterizeGeometryOpts struct {
errorHandler ErrorHandler
}

type gridCreateOpts struct {
errorHandler ErrorHandler
}

// GridCreateOption is an option that can be passed to GridCreate
// Available options are:
// - none yet
type GridCreateOption interface {
setGridCreateOpt(mo *gridCreateOpts)
}

// RasterizeGeometryOption is an option that can be passed tp Dataset.RasterizeGeometry()
type RasterizeGeometryOption interface {
setRasterizeGeometryOpt(o *rasterizeGeometryOpts)
Expand Down

0 comments on commit 5f05eb4

Please sign in to comment.