diff --git a/build/package/Dockerfile b/build/package/Dockerfile index e7daab1..7de75fe 100644 --- a/build/package/Dockerfile +++ b/build/package/Dockerfile @@ -1,5 +1,8 @@ # Build context is the repo root: docker build -f build/package/Dockerfile . -FROM golang:1.26.1-alpine AS builder +FROM golang:1.26-alpine AS builder + +# godal (GDAL bindings) needs cgo, a C toolchain and the GDAL headers. +RUN apk add --no-cache gdal-dev gcc musl-dev pkgconfig WORKDIR /src @@ -7,13 +10,13 @@ COPY go.mod go.sum ./ RUN go mod download COPY . . -RUN CGO_ENABLED=0 GOOS=linux go build -trimpath -o /out/gis ./cmd/gis +RUN CGO_ENABLED=1 GOOS=linux go build -trimpath -o /out/gis ./cmd/gis -FROM alpine:3.20 +FROM alpine:3.22 -# gdal-tools provides gdal_translate / gdalinfo for raster COG conversion and -# ogr2ogr for extracting vector feature geometry (worker). -RUN apk add --no-cache ca-certificates tzdata gdal-tools \ +# gdal: shared library for godal (raster math in the worker). +# gdal-tools: gdal_translate / gdalinfo (COG conversion) and ogr2ogr. +RUN apk add --no-cache ca-certificates tzdata gdal gdal-tools \ && adduser -D -u 10001 app WORKDIR /app diff --git a/go.mod b/go.mod index 27b678d..e06f590 100644 --- a/go.mod +++ b/go.mod @@ -5,6 +5,7 @@ go 1.26.1 toolchain go1.26.2 require ( + github.com/airbusgeo/godal v0.0.18 github.com/caarlos0/env/v11 v11.4.1 github.com/go-chi/chi/v5 v5.3.0 github.com/go-playground/validator/v10 v10.30.3 diff --git a/go.sum b/go.sum index c5340a8..f51c393 100644 --- a/go.sum +++ b/go.sum @@ -1,3 +1,5 @@ +github.com/airbusgeo/godal v0.0.18 h1:NQJ9P7OgSs5zVtlWH/tOta7QHYbbtha7Nttf63YM01Q= +github.com/airbusgeo/godal v0.0.18/go.mod h1:Mpn/6exBzEbIx4DaG2rIjCmAiF9Ww4SJN2ZtwVsSwCo= github.com/caarlos0/env/v11 v11.4.1 h1:fYwH0sWEsBSMPG7t4e/PEfTFzrWrpjyygXyUnWiSwEw= github.com/caarlos0/env/v11 v11.4.1/go.mod h1:qupehSf/Y0TUTsxKywqRt/vJjN5nz6vauiYEUUr8P4U= github.com/cespare/xxhash/v2 v2.3.0 h1:UL815xU9SqsFlibzuggzjXhog7bL6oX9BbNZnL2UFvs= diff --git a/internal/raster/godal.go b/internal/raster/godal.go new file mode 100644 index 0000000..904f44d --- /dev/null +++ b/internal/raster/godal.go @@ -0,0 +1,133 @@ +package raster + +import ( + "context" + "fmt" + "math" + "sync" + + "github.com/airbusgeo/godal" +) + +var registerOnce sync.Once + +// GodalIO implements IO with the GDAL bindings (github.com/airbusgeo/godal). +// GDAL datasets are not goroutine-safe: callers must not share a Reader or +// Writer between goroutines without external synchronisation. +type GodalIO struct{} + +// NewGodalIO registers GDAL drivers (once) and returns a GodalIO. +func NewGodalIO() *GodalIO { + registerOnce.Do(godal.RegisterAll) + return &GodalIO{} +} + +// DefaultWriterOptions returns the Landscape Program output profile. +func DefaultWriterOptions() WriterOptions { + return WriterOptions{NoData: -9999, TileSize: 512, Compression: "DEFLATE"} +} + +type godalReader struct { + ds *godal.Dataset + grid Grid +} + +func (g *GodalIO) Open(_ context.Context, path string) (Reader, error) { + ds, err := godal.Open(path) + if err != nil { + return nil, fmt.Errorf("open raster %s: %w", path, err) + } + st := ds.Structure() + if st.NBands < 1 { + ds.Close() + return nil, fmt.Errorf("open raster %s: no bands", path) + } + gt, err := ds.GeoTransform() + if err != nil { + ds.Close() + return nil, fmt.Errorf("open raster %s: geotransform: %w", path, err) + } + return &godalReader{ds: ds, grid: Grid{Width: st.SizeX, Height: st.SizeY, GeoTransform: gt, ProjectionWKT: ds.Projection()}}, nil +} + +func (r *godalReader) Grid() Grid { return r.grid } + +func (r *godalReader) NoData() (float64, bool) { return r.ds.Bands()[0].NoData() } + +func (r *godalReader) ReadFloat64(w Window, dst []float64) error { + if len(dst) < w.Width*w.Height { + return fmt.Errorf("read window: dst has %d elements, need %d", len(dst), w.Width*w.Height) + } + return r.ds.Bands()[0].Read(w.X, w.Y, dst[:w.Width*w.Height], w.Width, w.Height) +} + +func (r *godalReader) Close() error { return r.ds.Close() } + +type godalWriter struct{ ds *godal.Dataset } + +func (g *GodalIO) Create(_ context.Context, path string, grid Grid, o WriterOptions) (Writer, error) { + tile := fmt.Sprint(o.TileSize) + ds, err := godal.Create(godal.GTiff, path, 1, godal.Int16, grid.Width, grid.Height, + godal.CreationOption("TILED=YES", "BLOCKXSIZE="+tile, "BLOCKYSIZE="+tile, + "COMPRESS="+o.Compression, "PREDICTOR=2", "BIGTIFF=IF_SAFER")) + if err != nil { + return nil, fmt.Errorf("create raster %s: %w", path, err) + } + if err := ds.SetGeoTransform(grid.GeoTransform); err != nil { + ds.Close() + return nil, fmt.Errorf("create raster %s: geotransform: %w", path, err) + } + if grid.ProjectionWKT != "" { + if err := ds.SetProjection(grid.ProjectionWKT); err != nil { + ds.Close() + return nil, fmt.Errorf("create raster %s: projection: %w", path, err) + } + } + if err := ds.Bands()[0].SetNoData(float64(o.NoData)); err != nil { + ds.Close() + return nil, fmt.Errorf("create raster %s: nodata: %w", path, err) + } + return &godalWriter{ds: ds}, nil +} + +func (w *godalWriter) WriteInt16(win Window, src []int16) error { + if len(src) < win.Width*win.Height { + return fmt.Errorf("write window: src has %d elements, need %d", len(src), win.Width*win.Height) + } + return w.ds.Bands()[0].Write(win.X, win.Y, src[:win.Width*win.Height], win.Width, win.Height) +} + +func (w *godalWriter) Close() error { return w.ds.Close() } + +// SameGrid reports whether two grids describe the same pixels: equal size, +// geotransform within 1e-9 (absolute), and the same CRS (GDAL IsSame; an +// empty projection only matches an empty projection). +// +// ProjectionWKT is parsed with godal.NewSpatialRef rather than +// NewSpatialRefFromWKT: despite the field name, callers may pass any user +// projection string GDAL's SetProjection accepts (e.g. "EPSG:32642"), not +// only WKT, and NewSpatialRefFromWKT rejects non-WKT input. +func (g *GodalIO) SameGrid(a, b Grid) bool { + if a.Width != b.Width || a.Height != b.Height { + return false + } + for i := range a.GeoTransform { + if math.Abs(a.GeoTransform[i]-b.GeoTransform[i]) > 1e-9 { + return false + } + } + if a.ProjectionWKT == "" || b.ProjectionWKT == "" { + return a.ProjectionWKT == b.ProjectionWKT + } + sa, err := godal.NewSpatialRef(a.ProjectionWKT) + if err != nil { + return false + } + defer sa.Close() + sb, err := godal.NewSpatialRef(b.ProjectionWKT) + if err != nil { + return false + } + defer sb.Close() + return sa.IsSame(sb) +} diff --git a/internal/raster/godal_test.go b/internal/raster/godal_test.go new file mode 100644 index 0000000..783d826 --- /dev/null +++ b/internal/raster/godal_test.go @@ -0,0 +1,62 @@ +package raster + +import ( + "context" + "path/filepath" + "testing" +) + +func TestGodalIO_RoundTrip(t *testing.T) { + rio := NewGodalIO() + ctx := context.Background() + path := filepath.Join(t.TempDir(), "out.tif") + grid := Grid{Width: 10, Height: 7, GeoTransform: [6]float64{100, 30, 0, 200, 0, -30}, ProjectionWKT: "EPSG:32642"} + + w, err := rio.Create(ctx, path, grid, DefaultWriterOptions()) + if err != nil { + t.Fatal(err) + } + full := make([]int16, 10*7) + for i := range full { + full[i] = int16(i%6) - 1 // -1..4 (-1 stands in for garbage, not nodata) + } + full[0] = -9999 + if err := w.WriteInt16(Window{0, 0, 10, 7}, full); err != nil { + t.Fatal(err) + } + if err := w.Close(); err != nil { + t.Fatal(err) + } + + r, err := rio.Open(ctx, path) + if err != nil { + t.Fatal(err) + } + defer r.Close() + g := r.Grid() + if g.Width != 10 || g.Height != 7 || g.GeoTransform != grid.GeoTransform || g.ProjectionWKT == "" { + t.Fatalf("grid mismatch: %+v", g) + } + if nd, ok := r.NoData(); !ok || nd != -9999 { + t.Fatalf("nodata = %v %v", nd, ok) + } + dst := make([]float64, 3*2) + if err := r.ReadFloat64(Window{X: 8, Y: 5, Width: 2, Height: 2}, dst[:4]); err != nil { + t.Fatal(err) + } + // pixel (8,5) has index 58 -> 58%6-1 = 3 ; (9,5)=59 -> 4 ; (8,6)=68 -> 1 ; (9,6)=69 -> 2 + want := []float64{3, 4, 1, 2} + for i, v := range want { + if dst[i] != v { + t.Fatalf("pixel %d = %v want %v", i, dst[i], v) + } + } + if !rio.SameGrid(g, grid) { + t.Fatal("SameGrid should be true for identical grids") + } + other := grid + other.GeoTransform[0] += 1 + if rio.SameGrid(g, other) { + t.Fatal("SameGrid should be false for shifted origin") + } +} diff --git a/internal/raster/io.go b/internal/raster/io.go new file mode 100644 index 0000000..05783af --- /dev/null +++ b/internal/raster/io.go @@ -0,0 +1,77 @@ +package raster + +import "context" + +// Grid describes a raster's pixel grid: its size, geotransform, and +// projection. Two grids describe the same pixels when IO.SameGrid reports +// true for them. +type Grid struct { + Width, Height int + GeoTransform [6]float64 + ProjectionWKT string +} + +// Window is a rectangular block of pixels within a raster, in pixel/line +// coordinates with (0,0) at the top-left corner. +type Window struct { + X, Y, Width, Height int +} + +// Reader reads pixel data and metadata from an existing raster. +type Reader interface { + // Grid returns the raster's pixel grid. + Grid() Grid + // NoData returns the raster's nodata value and whether one is set. + NoData() (float64, bool) + // ReadFloat64 reads the pixels in w into dst, row-major, converting the + // underlying pixel type to float64. len(dst) must be at least + // w.Width*w.Height. + ReadFloat64(w Window, dst []float64) error + // Close releases the underlying dataset handle. + Close() error +} + +// WriterOptions configures a raster created by IO.Create. +type WriterOptions struct { + // NoData is the nodata value written to the output band. + NoData int16 + // TileSize is the internal tile size (both block dimensions) in pixels. + TileSize int + // Compression is the GDAL COMPRESS creation option value. + Compression string +} + +// Writer writes pixel data to a raster created by IO.Create. +type Writer interface { + // WriteInt16 writes src, row-major, into the pixels described by w. + // len(src) must be at least w.Width*w.Height. + WriteInt16(w Window, src []int16) error + // Close flushes and releases the underlying dataset handle. + Close() error +} + +// IO is an engine-agnostic block I/O abstraction over raster storage. It lets +// callers open and create rasters without depending on a particular GDAL +// binding. +type IO interface { + // Open opens the raster at path for reading. + Open(ctx context.Context, path string) (Reader, error) + // Create creates a new single-band int16 raster at path with the given + // grid and writer options, ready for WriteInt16 calls. + Create(ctx context.Context, path string, g Grid, o WriterOptions) (Writer, error) + // SameGrid reports whether a and b describe the same pixel grid: equal + // size, geotransform, and coordinate reference system. + SameGrid(a, b Grid) bool +} + +// Windows tiles a w×h raster into row-major windows of at most size×size. +func Windows(w, h, size int) []Window { + var out []Window + for y := 0; y < h; y += size { + wh := min(size, h-y) + for x := 0; x < w; x += size { + out = append(out, Window{X: x, Y: y, Width: min(size, w-x), Height: wh}) + } + } + return out +} diff --git a/internal/raster/io_test.go b/internal/raster/io_test.go new file mode 100644 index 0000000..b0ddf1a --- /dev/null +++ b/internal/raster/io_test.go @@ -0,0 +1,61 @@ +package raster + +import ( + "reflect" + "testing" +) + +func TestWindows(t *testing.T) { + tests := []struct { + name string + w, h, size int + want []Window + }{ + { + name: "10x10 size4 tiles into 9 windows", + w: 10, h: 10, size: 4, + want: []Window{ + {X: 0, Y: 0, Width: 4, Height: 4}, + {X: 4, Y: 0, Width: 4, Height: 4}, + {X: 8, Y: 0, Width: 2, Height: 4}, + {X: 0, Y: 4, Width: 4, Height: 4}, + {X: 4, Y: 4, Width: 4, Height: 4}, + {X: 8, Y: 4, Width: 2, Height: 4}, + {X: 0, Y: 8, Width: 4, Height: 2}, + {X: 4, Y: 8, Width: 4, Height: 2}, + {X: 8, Y: 8, Width: 2, Height: 2}, + }, + }, + { + name: "5x3 size8 fits in one window", + w: 5, h: 3, size: 8, + want: []Window{ + {X: 0, Y: 0, Width: 5, Height: 3}, + }, + }, + } + + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + got := Windows(tt.w, tt.h, tt.size) + if !reflect.DeepEqual(got, tt.want) { + t.Fatalf("Windows(%d, %d, %d) = %+v, want %+v", tt.w, tt.h, tt.size, got, tt.want) + } + if len(got) > 0 { + last := got[len(got)-1] + wantLast := tt.want[len(tt.want)-1] + if last != wantLast { + t.Fatalf("last window = %+v, want %+v", last, wantLast) + } + } + }) + } +} + +func TestDefaultWriterOptions(t *testing.T) { + got := DefaultWriterOptions() + want := WriterOptions{NoData: -9999, TileSize: 512, Compression: "DEFLATE"} + if got != want { + t.Fatalf("DefaultWriterOptions() = %+v, want %+v", got, want) + } +}