使用godal进行地理空间数据处理

godal是Go语言对GDAL(Geospatial Data Abstraction Library)的绑定，用于处理各种栅格和矢量地理空间数据格式。下面我将介绍godal的基本用法和示例代码。

安装

首先需要安装GDAL库，然后安装godal：

# Ubuntu/Debian
sudo apt-get install gdal-bin libgdal-dev

# macOS
brew install gdal

# 然后安装godal
go get -u github.com/airbusgeo/godal

基本用法

1. 打开和读取栅格数据

package main

import (
	"fmt"
	"github.com/airbusgeo/godal"
)

func main() {
	// 初始化GDAL
	godal.RegisterAll()

	// 打开栅格文件
	ds, err := godal.Open("example.tif")
	if err != nil {
		panic(err)
	}
	defer ds.Close()

	// 获取栅格信息
	fmt.Printf("Raster size: %dx%d\n", ds.RasterXSize(), ds.RasterYSize())
	fmt.Printf("Number of bands: %d\n", ds.RasterCount())

	// 读取第一个波段
	band := ds.Bands()[0]
	blockX, blockY := band.BlockSize()
	fmt.Printf("Band block size: %dx%d\n", blockX, blockY)
	fmt.Printf("Data type: %s\n", band.DataType().Name())
}

2. 读取矢量数据

package main

import (
	"fmt"
	"github.com/airbusgeo/godal"
)

func main() {
	godal.RegisterAll()

	// 打开矢量文件
	ds, err := godal.Open("example.shp")
	if err != nil {
		panic(err)
	}
	defer ds.Close()

	// 获取图层
	layer := ds.Layers()[0]
	fmt.Printf("Layer name: %s\n", layer.Name())
	fmt.Printf("Feature count: %d\n", layer.FeatureCount())

	// 遍历要素
	layer.ResetReading()
	for {
		feature := layer.NextFeature()
		if feature == nil {
			break
		}
		geom := feature.Geometry()
		if geom != nil {
			fmt.Printf("Geometry type: %s\n", geom.Type().Name())
		}
	}
}

3. 创建和写入栅格数据

package main

import (
	"github.com/airbusgeo/godal"
)

func main() {
	godal.RegisterAll()

	// 创建驱动
	driver, err := godal.GetDriverByName("GTiff")
	if err != nil {
		panic(err)
	}

	// 创建新数据集
	ds, err := driver.Create("output.tif", 100, 100, 1, godal.Byte, nil)
	if err != nil {
		panic(err)
	}
	defer ds.Close()

	// 设置地理变换和投影
	transform := [6]float64{440720.0, 60.0, 0.0, 3751320.0, 0.0, -60.0}
	ds.SetGeoTransform(transform[:])
	ds.SetProjection(`PROJCS["NAD27 / UTM zone 11N",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Clarke 1866",6378206.4,294.9786982138982,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG","6267"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","26711"]]`)

	// 写入数据
	band := ds.Bands()[0]
	data := make([]byte, 100*100)
	for i := range data {
		data[i] = byte(i % 256)
	}
	err = band.IO(godal.Write, 0, 0, 100, 100, data, 100, 100)
	if err != nil {
		panic(err)
	}
}

4. 栅格数据处理示例

package main

import (
	"fmt"
	"github.com/airbusgeo/godal"
)

func main() {
	godal.RegisterAll()

	// 打开输入文件
	srcDS, err := godal.Open("input.tif")
	if err != nil {
		panic(err)
	}
	defer srcDS.Close()

	// 创建输出文件
	driver, _ := godal.GetDriverByName("GTiff")
	dstDS, err := driver.Create("output.tif", srcDS.RasterXSize(), srcDS.RasterYSize(), srcDS.RasterCount(), srcDS.Bands()[0].DataType(), nil)
	if err != nil {
		panic(err)
	}
	defer dstDS.Close()

	// 复制地理信息
	gt := srcDS.GeoTransform()
	dstDS.SetGeoTransform(gt[:])
	dstDS.SetProjection(srcDS.Projection())

	// 处理每个波段
	for i := 1; i <= srcDS.RasterCount(); i++ {
		srcBand := srcDS.Bands()[i-1]
		dstBand := dstDS.Bands()[i-1]

		// 读取数据
		data := make([]float32, srcDS.RasterXSize()*srcDS.RasterYSize())
		err = srcBand.IO(godal.Read, 0, 0, srcDS.RasterXSize(), srcDS.RasterYSize(), data, srcDS.RasterXSize(), srcDS.RasterYSize())
		if err != nil {
			panic(err)
		}

		// 简单处理数据 - 这里示例乘以2
		for j := range data {
			data[j] *= 2
		}

		// 写入数据
		err = dstBand.IO(godal.Write, 0, 0, srcDS.RasterXSize(), srcDS.RasterYSize(), data, srcDS.RasterXSize(), srcDS.RasterYSize())
		if err != nil {
			panic(err)
		}
	}

	fmt.Println("Processing completed!")
}

高级功能

1. 重投影

package main

import (
	"github.com/airbusgeo/godal"
)

func main() {
	godal.RegisterAll()

	// 打开源数据集
	srcDS, err := godal.Open("input.tif")
	if err != nil {
		panic(err)
	}
	defer srcDS.Close()

	// 定义目标投影 (WGS84)
	dstSRS := godal.NewSpatialReference()
	dstSRS.FromEPSG(4326)

	// 创建重投影选项
	options := []string{
		"RESAMPLING=NEAREST",
		"COMPRESS=DEFLATE",
	}

	// 执行重投影
	err = godal.Warp("output.tif", []godal.Dataset{srcDS}, options, dstSRS)
	if err != nil {
		panic(err)
	}
}

2. 栅格计算

package main

import (
	"github.com/airbusgeo/godal"
)

func main() {
	godal.RegisterAll()

	// 打开两个输入栅格
	ds1, err := godal.Open("input1.tif")
	if err != nil {
		panic(err)
	}
	defer ds1.Close()

	ds2, err := godal.Open("input2.tif")
	if err != nil {
		panic(err)
	}
	defer ds2.Close()

	// 创建计算选项
	options := []string{
		"FORMAT=GTiff",
		"OUTPUT_TYPE=Float32",
		"NO_DATA=-9999",
	}

	// 执行栅格计算 (这里示例相加)
	err = godal.Calc("output.tif", []godal.Dataset{ds1, ds2}, "A+B", options)
	if err != nil {
		panic(err)
	}
}

注意事项

1. 确保在使用前调用godal.RegisterAll()注册所有驱动
2. 记得关闭数据集(defer ds.Close())
3. 处理大数据时注意内存使用，可以使用分块处理
4. 错误处理很重要，GDAL操作可能会因为各种原因失败

godal提供了对GDAL功能的全面访问，可以处理几乎所有常见的地理空间数据格式。更多高级功能可以参考GDAL官方文档和godal的API文档。