Add GDALInvGeoTransform

[dzfranklin]

I'm very new to gdal. I want to lookup a wgs84 coordinate in a dataset already in the wgs84 projection. Based on reading the source of gdallocationinfo I think I need to invert the geotransform matrix.

Did I miss a better way to do this?

I couldn't find that function in godal. It's just some straightforward math so I think it makes sense to port it to go. Here's my first draft (not fully tested).

This isn't the most idiomatic api but since Dataset.GeoTransform returns a [6]float64 I don't think it can be a method.

/*
InvGeoTransform inverts a Geotransform.

This function will invert a standard 3x2 set of GeoTransform coefficients.
This converts the equation from being pixel to geo to being geo to pixel.

- gt_in: Input geotransform (six doubles - unaltered)

- gt_out: Output geotransform (six doubles - updated)

Panics if the equation is uninvertable.
*/
func InvGeoTransform(input [6]float64) [6]float64 {
    // Ported from GDALInvGeoTransform

    out := [6]float64{}

    // Special case - no rotation - to avoid computing determinate
    // and potential precision issues.
    if input[2] == 0.0 && input[4] == 0.0 && input[1] != 0.0 &&
        input[5] != 0.0 {
        /*X = input[0] + x * input[1]
          Y = input[3] + y * input[5]
          -->
          x = -input[0] / input[1] + (1 / input[1]) * X
          y = -input[3] / input[5] + (1 / input[5]) * Y
        */
        out[0] = -input[0] / input[1]
        out[1] = 1.0 / input[1]
        out[2] = 0.0
        out[3] = -input[3] / input[5]
        out[4] = 0.0
        out[5] = 1.0 / input[5]
        return out
    }

    // Assume a 3rd row that is [1 0 0].

    // Compute determinate.

    det := input[1]*input[5] - input[2]*input[4]
    magnitude := max(max(math.Abs(input[1]), math.Abs(input[2])),
        max(math.Abs(input[4]), math.Abs(input[5])))

    if math.Abs(det) <= 1e-10*magnitude*magnitude {
        panic("equation is uninvertable")
    }

    invDet := 1.0 / det

    // Compute adjoint, and divide by determinate.

    out[1] = input[5] * invDet
    out[4] = -input[4] * invDet

    out[2] = -input[2] * invDet
    out[5] = input[1] * invDet

    out[0] = (input[2]*input[3] - input[0]*input[5]) * invDet
    out[3] = (-input[1]*input[3] + input[0]*input[4]) * invDet

    return out
}

[tbonfort]

GDAL probably already an API invert a geotransform, and exposing that function in godal would be the way forward (as godal is a gdal wrapper, not a gdal reimplementation in go)