Convert a raster layer to a vector layer, by creating geojson polygon features for each individual pixel's extent in the raster layer. Equivalent to QGIS' Raster pixels to polygons Algorithm

raster-pixels-to-polygons.go

package main

import (
	"errors"
	"flag"
	"math"
	"os"

	"github.com/lukeroth/gdal"
	geojson "github.com/paulmach/go.geojson"
	log "github.com/sirupsen/logrus"
)

func main() {

	var filePath string
	flag.StringVar(&filePath, "f", "", "Raster file")

	var band int
	flag.IntVar(&band, "b", 1, "Raster Band to extract")

	var outFile string
	flag.StringVar(&outFile, "o", "", "Output file")

	flag.Parse()

	if _, err := os.Stat(filePath); errors.Is(err, os.ErrNotExist) {
		if err != nil {
			log.Error(err)
			return
		}
	}

	if outFile == "" {
		log.Errorf("provide output file")
		return
	}

	log.Infof("Converting '%s' using band '%d' ", filePath, band)

	dataset, err := gdal.Open(filePath, gdal.ReadOnly)

	if err != nil {
		log.Error(err)
		return
	}

	defer dataset.Close()

	log.Infof("Opened dataset using '%s' driver", dataset.Driver().ShortName())

	bandsCount := dataset.RasterCount()

	if band > bandsCount {
		log.Errorf("selected band : %d is out of available bands range : 1 - %d ", band, bandsCount)
		return
	}

	xSize := dataset.RasterXSize()
	ySize := dataset.RasterYSize()

	geot := dataset.GeoTransform()

	rasterUnitsPerPixelX := math.Abs(geot[1])
	rasterUnitsPerPixelY := math.Abs(geot[5])

	xMin := math.Min(geot[0], geot[0]+float64(xSize)*geot[1])
	xMax := math.Max(geot[0], geot[0]+float64(xSize)*geot[1])
	yMin := math.Min(geot[3], geot[3]+float64(ySize)*geot[5])
	yMax := math.Max(geot[3], geot[3]+float64(ySize)*geot[5])
	_ = []float64{xMin, yMin, xMax, yMax}

	currentY := yMax - 0.5*rasterUnitsPerPixelY

	raster := dataset.RasterBand(1)
	noDataValue, _ := raster.NoDataValue()

	buffer := make([]float64, xSize*ySize)

	raster.IO(gdal.Read, 0, 0, xSize, ySize, buffer, xSize, ySize, 0, 0)

	featureCollection := geojson.NewFeatureCollection()

	for row := 0; row < ySize; row++ {
		currentX := xMin + 0.5*rasterUnitsPerPixelX

		for col := 0; col < xSize; col++ {
			pixelVal := buffer[col+row*xSize]

			if pixelVal != noDataValue {
				geometry := createGeometryForPixel(currentX, currentY, rasterUnitsPerPixelX, rasterUnitsPerPixelY)
				feature := geojson.NewFeature(geometry)
				feature.Properties = map[string]interface{}{
					"value": pixelVal,
				}
				featureCollection.AddFeature(feature)
			}

			currentX += rasterUnitsPerPixelX
		}

		currentY -= rasterUnitsPerPixelY
	}

	rawJSON, err := featureCollection.MarshalJSON()

	if err != nil {
		log.Error(err)
		return
	}

	log.Infof("Saving geojson output to '%s' ", outFile)

	f, err := os.Create(outFile)

	if err != nil {
		log.Error(err)
		return
	}

	_, err = f.Write(rawJSON)

	if err != nil {
		log.Error(err)
		return
	}

	log.Infof("Done! We have put some good stuff at '%s' ", outFile)
}

func createGeometryForPixel(centerX float64, centerY float64, pixelWidthX float64, pixelWidthY float64) *geojson.Geometry {

	hCellSizeX := pixelWidthX / 2.0
	hCellSizeY := pixelWidthY / 2.0

	xMin := centerX - hCellSizeX //west
	yMin := centerY - hCellSizeY // south
	xMax := centerX + hCellSizeX // east
	yMax := centerY + hCellSizeY // north

	lowLeft := []float64{xMin, yMin}
	topLeft := []float64{xMin, yMax}
	topRight := []float64{xMax, yMax}
	lowRight := []float64{xMax, yMin}

	polygon := [][][]float64{{lowLeft, lowRight, topRight, topLeft, lowLeft}}

	return geojson.NewPolygonGeometry(polygon)
}