read_elevation_geotiff.go

package main

import (
	"fmt"

	"github.com/lukeroth/gdal"
	"github.com/paulmach/orb"
	"github.com/paulmach/orb/project"
)

func main() {
	// INPUTS
	fname := fmt.Sprintf("elevation_tiles/%v/%v/%v.tif", 10, 163, 395)
	point := orb.Point{-122.4472942, 37.7543107} // The top of Twin Peaks in SF

	// Open the GeoTiff
	dataset, err := gdal.Open(fname, gdal.ReadOnly)
	if err != nil {
		panic(err)
	}
	defer dataset.Close()

	// Read the first raster band
	band := dataset.RasterBand(1)

	// Read the transform values for this tile
	gt := dataset.GeoTransform()

	minx := gt[0]
	maxy := gt[3]
	pixelWidth := gt[1]

	// Calculate the pixel offsets from the top left of the dataset
	mercPoint := project.WGS84.ToMercator(point)
	xOff := int((mercPoint.X() - minx) / pixelWidth)
	yOff := int((maxy - mercPoint.Y()) / pixelWidth)

	// Read the pixel value at the calculated offset
	buf := make([]int16, 1)
	err = band.IO(gdal.Read, xOff, yOff, 1, 1, buf, 1, 1, 0, 0)
	if err != nil {
		panic(err)
	}

	// Print the answer
	fmt.Printf("elevation: meters=%v @ lon=%v lat=%v\n", buf[0], point.X(), point.Y())
}