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)