test_geom.go

package main

// #include "gdal.h"
// #include "gdal_alg.h"
// #include "ogr_api.h"
// #include "ogr_srs_api.h"
// #include "cpl_string.h"
// #cgo LDFLAGS: -lgdal
import "C"

import (
	"fmt"
	"math"
	"unsafe"
)

var cWGS84WKT *C.char = C.CString(`GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AUTHORITY["EPSG","4326"]]","proj4":"+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs `)

func envelopePolygon(hDS C.GDALDatasetH) C.OGRGeometryH {
        geoTrans := make([]float64, 6)
        C.GDALGetGeoTransform(hDS, (*C.double)(&geoTrans[0]))

        var ulX, ulY C.double
        C.GDALApplyGeoTransform((*C.double)(&geoTrans[0]), C.double(0), C.double(0), &ulX, &ulY)
        var lrX, lrY C.double
        C.GDALApplyGeoTransform((*C.double)(&geoTrans[0]), C.double(C.GDALGetRasterXSize(hDS)), C.double(C.GDALGetRasterYSize(hDS)), &lrX, &lrY)

        polyWKT := fmt.Sprintf("POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))", ulX, ulY,
                ulX, lrY,
                lrX, lrY,
                lrX, ulY,
                ulX, ulY)

        ppszData := C.CString(polyWKT)
        defer C.free(unsafe.Pointer(ppszData))

        var hGeom C.OGRGeometryH
        hSRS := C.OSRNewSpatialReference(C.GDALGetProjectionRef(hDS))
        _ = C.OGR_G_CreateFromWkt(&ppszData, hSRS, &hGeom)
        return hGeom
}

func getDrillFileDescriptor(ds C.GDALDatasetH, g C.OGRGeometryH) {
        gCopy := C.OGR_G_Clone(g)

        if C.GoString(C.GDALGetProjectionRef(ds)) != "" {
                desSRS := C.OSRNewSpatialReference(C.GDALGetProjectionRef(ds))
                defer C.OSRDestroySpatialReference(desSRS)
                srcSRS := C.OSRNewSpatialReference(cWGS84WKT)
                defer C.OSRDestroySpatialReference(srcSRS)
                trans := C.OCTNewCoordinateTransformation(srcSRS, desSRS)
                C.OGR_G_Transform(gCopy, trans)
        }

        fileEnv := envelopePolygon(ds)
        var fileWkt *C.char
        C.OGR_G_ExportToWkt(fileEnv, &fileWkt)
        inters := C.OGR_G_Intersection(gCopy, fileEnv)
        var intersWkt *C.char
        C.OGR_G_ExportToWkt(inters, &intersWkt)

        var env C.OGREnvelope
        C.OGR_G_GetEnvelope(inters, &env)

        geot := make([]float64, 6)
        C.GDALGetGeoTransform(ds, (*C.double)(&geot[0]))

        invGeot := make([]float64, 6)
        C.GDALInvGeoTransform((*C.double)(&geot[0]), (*C.double)(&invGeot[0]))

        var offMinX, offMinY, offMaxX, offMaxY C.double
        C.GDALApplyGeoTransform((*C.double)(&invGeot[0]), env.MinX, env.MinY, &offMinX, &offMinY)
        C.GDALApplyGeoTransform((*C.double)(&invGeot[0]), env.MaxX, env.MaxY, &offMaxX, &offMaxY)

        offsetX := int(math.Min(float64(offMinX), float64(offMaxX)))
        offsetY := int(math.Min(float64(offMinY), float64(offMaxY)))
        countX := int(math.Max(float64(offMinX), float64(offMaxX))) - offsetX
        countY := int(math.Max(float64(offMinY), float64(offMaxY))) - offsetY
        if countX == 0 {
                countX++
        }
        if countY == 0 {
                countY++
        }

        return
}

func main() {
	C.GDALAllRegister()
	path := "/g/data3/fr5/prl900/FC.v302.MCD43A4/FC.v302.MCD43A4.h12v03.2012.005.nc"
	cPath := C.CString(path)
	defer C.free(unsafe.Pointer(cPath))
	ds := C.GDALOpen(cPath, C.GDAL_OF_READONLY)
	if ds == nil {
		fmt.Printf("GDAL could not open dataset: %s", path)
		return
	}
	defer C.GDALClose(ds)

	geom := envelopePolygon(ds)

	getDrillFileDescriptor(ds, geom)
}