Example of reading multi-dimensional NetCDF files using ``fhs/go-netcdf`` library

main.go

package main

import (
	"fmt"
	"github.com/fhs/go-netcdf/netcdf"
	"iter"
	"math"
	"os"
	"strings"
)

const ncLength = 7200
const ncWidth = 2000
const ncMinLat = float32(-180.0)
const ncMaxLat = float32(180)
const ncMinLon = float32(-50)
const ncMaxLon = float32(50)
const ncLonStep = (ncMaxLon - ncMinLon) / (ncWidth - 1)
const ncLatStep = (ncMaxLat - ncMinLat) / (ncLength - 1)

func ncCoordsIndex(lon, lat float32) (int, int, error) {
	if lon < ncMinLon || lon > ncMaxLon || lat < ncMinLat || lat > ncMaxLat {
		return 0, 0, fmt.Errorf("coordinates %f, %f are out of bounds", lon, lat)
	}

	lonIndex := int(math.Round(float64((lon - ncMinLon) / ncLonStep)))
	latIndex := int(math.Round(float64((lat - ncMinLat) / ncLatStep)))

	return lonIndex, latIndex, nil
}

type PrecipData struct {
	Date   int
	Precip float32
	Lon    float32
	Lat    float32
}

func ncReadTimeRange(nc *netcdf.Dataset) (int, error) {
	dTime, err := nc.Dim("time")
	if err != nil {
		return 0, err
	}

	timeRange, err := dTime.Len()
	if err != nil {
		return 0, err
	}

	return int(timeRange), nil
}

func ncReadPrecip(nc *netcdf.Dataset, bottomLeftLon, bottomLeftLat, topRightLon, topRightLat float32) (iter.Seq[PrecipData], error) {
	timeRange, err := ncReadTimeRange(nc)
	if err != nil {
		return nil, err
	}

	endLonIndex, startLatIndex, err := ncCoordsIndex(bottomLeftLon, bottomLeftLat)
	if err != nil {
		return nil, err
	}

	startLonIndex, endLatIndex, err := ncCoordsIndex(topRightLon, topRightLat)
	if err != nil {
		return nil, err
	}

	latRange := endLatIndex - startLatIndex + 1
	lonRange := endLonIndex - startLonIndex + 1
	dataSize := timeRange * latRange * lonRange

	precipVar, err := nc.Var("precip")
	if err != nil {
		return nil, err
	}

	data := make([]float32, dataSize)
	err = precipVar.ReadFloat32Slice(
		data,
		[]uint64{0, uint64(startLonIndex), uint64(startLatIndex)},
		[]uint64{uint64(timeRange), uint64(lonRange), uint64(latRange)},
	)

	if err != nil {
		return nil, err
	}

	return func(yield func(PrecipData) bool) {
		for i, d := range data {
			lat := i % latRange
			lon := (i / latRange) % lonRange
			date := i/latRange/lonRange + 1

			if d == -9999 || date != 1 {
				continue
			}

			data := PrecipData{
				Date:   date,
				Precip: d,
				Lon:    topRightLon + float32(lon)*ncLonStep,
				Lat:    bottomLeftLat + float32(lat)*ncLatStep,
			}

			if !yield(data) {
				return
			}
		}
	}, nil
}

func main() {
	// https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_daily/tifs/p05/1985/chirps-v2.0.1985.01.01.tif.gz
	file, err := netcdf.OpenFile("/home/nathan/Projects/gssat/apps/databundler/chirps-db/1985.nc", netcdf.NOWRITE)
	if err != nil {
		panic(err)
	}
	defer file.Close()

	iter, err := ncReadPrecip(&file, 0, 0, -20, 20)
	if err != nil {
		panic(err)
	}

	var sb strings.Builder
	sb.WriteString("{ \"type\": \"FeatureCollection\", \"features\": [")
	for d := range iter {
		sb.WriteString(fmt.Sprintf("{\"type\": \"Feature\", \"properties\": { \"p\": %f }, \"geometry\": {\"type\": \"Point\", \"coordinates\": [%f, %f]}}", d.Precip, d.Lat, d.Lon))
		sb.WriteString(",")
	}

	str := sb.String()[0 : sb.Len()-1]
	str += "]}"

	os.WriteFile("data.json", []byte(str), 0644)
}

nohardcode.go

package main

import (
	"fmt"
	"github.com/fhs/go-netcdf/netcdf"
	"iter"
	"math"
	"os"
	"strings"
)

func readNcSize(nc *netcdf.Dataset) (length int, width int, days int, err error) {
	time, err := nc.Dim("time")
	if err != nil {
		return 0, 0, 0, err
	}
	t, err := time.Len()
	if err != nil {
		return 0, 0, 0, err
	}

	// Lon and lat are swapped on purpose
	lon, err := nc.Dim("latitude")
	if err != nil {
		return 0, 0, 0, err
	}
	i, err := lon.Len()
	if err != nil {
		return 0, 0, 0, err
	}

	// Lon and lat are swapped on purpose
	lat, err := nc.Dim("longitude")
	if err != nil {
		return 0, 0, 0, err
	}
	j, err := lat.Len()
	if err != nil {
		return 0, 0, 0, err
	}

	return int(i), int(j), int(t), nil
}

func readNcEnvelope(nc *netcdf.Dataset) (minLon float32, maxLon float32, minLat float32, maxLat float32, err error) {
	precipVar, err := nc.Var("precip")
	if err != nil {
		return 0, 0, 0, 0, err
	}

	// Typo is intentional, for some reason the CHIRPS netcdf files are full of typos
	attrs := []netcdf.Attr{
		precipVar.Attr("geostatial_lat_min"),
		precipVar.Attr("geostatial_lat_max"),
		precipVar.Attr("geostatial_lon_min"),
		precipVar.Attr("geostatial_lon_max"),
	}

	buf := make([]float32, 4)
	for i, attr := range attrs {
		if err := attr.ReadFloat32s(buf[i : i+1]); err != nil {
			return 0, 0, 0, 0, err
		}
	}

	// Lon and lat are swapped on purpose
	return buf[0], buf[1], buf[2], buf[3], nil
}

func ncCoordsIndex(nc *netcdf.Dataset, lon, lat float32) (int, int, error) {
	minLon, maxLon, minLat, maxLat, err := readNcEnvelope(nc)
	if err != nil {
		return 0, 0, err
	}

	width, length, _, err := readNcSize(nc)
	if err != nil {
		return 0, 0, err
	}

	if lon < minLon || lon > maxLon || lat < minLat || lat > maxLat {
		return 0, 0, fmt.Errorf("coordinates %f, %f are out of bounds", lon, lat)
	}

	ncLonStep := (maxLon - minLon) / float32(width-1)
	ncLatStep := (maxLat - minLat) / float32(length-1)

	lonIndex := int(math.Round(float64((lon - minLon) / ncLonStep)))
	latIndex := int(math.Round(float64((lat - minLat) / ncLatStep)))

	return lonIndex, latIndex, nil
}

type PrecipData struct {
	Date   int
	Precip float32
	Lon    float32
	Lat    float32
}

func ncReadPrecip(nc *netcdf.Dataset, bottomLeftLon, bottomLeftLat, topRightLon, topRightLat float32) (iter.Seq[PrecipData], error) {
	minLon, maxLon, minLat, maxLat, err := readNcEnvelope(nc)
	if err != nil {
		return nil, err
	}

	width, length, timeRange, err := readNcSize(nc)
	if err != nil {
		return nil, err
	}

	endLonIndex, startLatIndex, err := ncCoordsIndex(nc, bottomLeftLon, bottomLeftLat)
	if err != nil {
		return nil, err
	}

	startLonIndex, endLatIndex, err := ncCoordsIndex(nc, topRightLon, topRightLat)
	if err != nil {
		return nil, err
	}

	latRange := endLatIndex - startLatIndex + 1
	lonRange := endLonIndex - startLonIndex + 1
	dataSize := timeRange * latRange * lonRange

	precipVar, err := nc.Var("precip")
	if err != nil {
		return nil, err
	}

	data := make([]float32, dataSize)
	err = precipVar.ReadFloat32Slice(
		data,
		[]uint64{0, uint64(startLonIndex), uint64(startLatIndex)},
		[]uint64{uint64(timeRange), uint64(lonRange), uint64(latRange)},
	)

	if err != nil {
		return nil, err
	}

	ncLonStep := (maxLon - minLon) / float32(width-1)
	ncLatStep := (maxLat - minLat) / float32(length-1)

	return func(yield func(PrecipData) bool) {
		for i, d := range data {
			lat := i % latRange
			lon := (i / latRange) % lonRange
			date := i/latRange/lonRange + 1

			if d == -9999 || date != 1 {
				continue
			}

			data := PrecipData{
				Date:   date,
				Precip: d,
				Lon:    topRightLon + float32(lon)*ncLonStep,
				Lat:    bottomLeftLat + float32(lat)*ncLatStep,
			}

			if !yield(data) {
				return
			}
		}
	}, nil
}

func main() {
	// https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_daily/tifs/p05/1985/chirps-v2.0.1985.01.01.tif.gz
	file, err := netcdf.OpenFile("/home/nathan/Projects/gssat/apps/databundler/chirps-db/1985.nc", netcdf.NOWRITE)
	if err != nil {
		panic(err)
	}
	defer file.Close()

	iter, err := ncReadPrecip(&file, 0, 0, -20, 20)
	if err != nil {
		panic(err)
	}

	var sb strings.Builder
	sb.WriteString("{ \"type\": \"FeatureCollection\", \"features\": [")
	for d := range iter {
		sb.WriteString(fmt.Sprintf("{\"type\": \"Feature\", \"properties\": { \"p\": %f }, \"geometry\": {\"type\": \"Point\", \"coordinates\": [%f, %f]}}", d.Precip, d.Lat, d.Lon))
		sb.WriteString(",")
	}

	str := sb.String()[0 : sb.Len()-1]
	str += "]}"

	os.WriteFile("data.json", []byte(str), 0644)
}