当前位置：首页 > news >正文

Go+Gdal 完成高性能GIS数据空间分析

news 2025/8/28 13:04:37

概要

环境准备

技术流程

一、在golang中如何调用gdal

二、读取数据

三、执行空间分析

四、性能提升

小结

概要

Gdal库可以说是所有gis软件的基础，基本上现在所有的工业gis软件都是基于gdal开发的，其主要包括了栅格处理、矢量处理、坐标系处理所涉及的各类基础api。本研究主要使用go的高并发能力和能使用原生C语言的混合开发的能力来实现对硬件性能的极致开发，来实现对大体量gis数据的快速空间分析。

环境准备

Golang的版本我选择的1.24.2，Gdal的版本的3.11.3(兼容3.8.3-3.11.3之间的版本)，然后就是GCC我使用的是mingw64，具体版本是winlibs personal build version gcc-11.3.0-llvm-14.0.3-mingw-w64msvcrt-10.0.0-r3，目前我在windows，linux，android这三个环境都成功实现了编译。

Golang调用Gdal并进行项目打包有目前有两种方案，一种是使用动态链接库，该方法优点就是简单高效，将gdal库作为外部链接库调用，打包时间短，方便好用不折腾。缺点么就是需要把动态库环境一起打包，系统架构略显臃肿。另一种就是静态链接库，简单来说就是把gdal中的相关代码进行源码编译，并使用go语言的gc将相关代码一并打包到exe中，优点就是耦合性强，打包出来的二进制可执行程序直接就自带了gdal的全部功能，缺点也非常明显，打包时间过于长了，把gdal的源码编译进go中一般要几个小时，每次打包都是一种煎熬，还有就是打包环境非常之严格，稍微错一点你都得折腾半天。所以还是推荐使用动态链接库的方案，本文也主要接受动态链接库的方案。

windows系统中使用最简单的方案就是通过osgeo4w-setup.exe来安装gdal的动态库，下面是具体操作流程：

官网地址：Making sure you're not a bot!

下载安装包：

执行安装步骤无限下一步即可

然后如图所示选中，无限下一步等待gdal的动态库下载完成即可

Linux环境中则直接sudo apt install gdal-bin libgdal-dev

如果是arm架构的linux系统则直接 pkg install gdal 即可

安卓架构比较复杂暂时不做介绍

安装完成后，命令行输入 gdainfo --version 能正确弹出版本号信息说明gdal环境就已经搞定了。

技术流程

一、在golang中如何调用gdal

下面是项目架构，目前只完成了一小部分，后面完成后会正式上线github

在介绍之前做个吐槽，其实gdal官方对几种热门语言是单独做了api的

同时有一些其他的开发者也开发了其他语言的api接口，官方也进行了收录

但是搞笑的来了，第一个就是go语言的，我进了github仓库看了下，这东西几年前就没维护了，而且做的功能也非常少，基本就是在go中使用原生c语言接口调用go来封装的一系列函数，但是这些函数都是很基础的函数，实用性太低，并且这东西兼容的gdal版本是2.3.2，这都是8年前的版本了，代码拉下来完全调用不起，各种报错。基于这种情况，我直接完全放弃了这个库，选择了全部手搓这条不归路，哈哈哈。

在cgo_header 定义gdal的环境，并让其能兼容多种系统

下面是具体代码

package OSGEO/*
#cgo windows CFLAGS: -IC:/OSGeo4W/include -IC:/OSGeo4W/include/gdal
#cgo windows LDFLAGS: -LC:/OSGeo4W/lib -lgdal_i -lstdc++ -static-libgcc -static-libstdc++
#cgo linux CFLAGS: -I/usr/include/gdal
#cgo linux LDFLAGS: -L/usr/lib -lgdal -lstdc++
#cgo android CFLAGS: -I/data/data/com.termux/files/usr/include
#cgo android LDFLAGS: -L/data/data/com.termux/files/usr/lib -lgdal
#include <gdal.h>
#include <gdal_alg.h>
#include <ogr_api.h>
#include <ogr_srs_api.h>
#include <cpl_error.h>
#include <cpl_conv.h>
#include <cpl_string.h>
#include <stdlib.h>// 初始化GDAL并修复PROJ配置问题
void initializeGDALWithProjFix(const char* projDataPath, const char* shapeEncoding) {// 动态设置PROJ数据路径，支持自定义路径if (projDataPath && strlen(projDataPath) > 0) {CPLSetConfigOption("PROJ_DATA", projDataPath);      // 设置PROJ数据目录CPLSetConfigOption("PROJ_LIB", projDataPath);       // 设置PROJ库路径（兼容性）}// 强制使用传统的GIS轴顺序（经度，纬度），避免坐标轴混乱CPLSetConfigOption("OSR_DEFAULT_AXIS_MAPPING_STRATEGY", "TRADITIONAL_GIS_ORDER");// 动态设置Shapefile编码，支持不同字符集if (shapeEncoding && strlen(shapeEncoding) > 0) {CPLSetConfigOption("SHAPE_ENCODING", shapeEncoding); // 设置Shapefile文件编码}// 注册所有GDAL驱动程序，启用栅格数据支持GDALAllRegister();// 注册所有OGR驱动程序，启用矢量数据支持OGRRegisterAll();
}// 创建EPSG:4490坐标系并设置正确的轴顺序
OGRSpatialReferenceH createEPSG4490WithCorrectAxis() {// 创建新的空间参考系统对象OGRSpatialReferenceH hSRS = OSRNewSpatialReference(NULL);if (!hSRS) {return NULL;  // 内存分配失败时返回NULL}// 导入EPSG:4490坐标系定义（中国大地坐标系2000）if (OSRImportFromEPSG(hSRS, 4490) != OGRERR_NONE) {OSRDestroySpatialReference(hSRS);  // 失败时清理资源return NULL;                        // 返回NULL表示创建失败}// 设置传统的GIS轴顺序（经度在前，纬度在后）OSRSetAxisMappingStrategy(hSRS, OAMS_TRADITIONAL_GIS_ORDER);return hSRS;  // 返回成功创建的空间参考系统
}// 从十六进制WKB字符串创建几何对象
OGRGeometryH createGeometryFromWKBHex(const char* wkbHex) {// 验证输入参数有效性if (!wkbHex || strlen(wkbHex) == 0) {return NULL;  // 空输入时返回NULL}int hexLen = strlen(wkbHex);  // 获取十六进制字符串长度// 验证十六进制字符串长度必须为偶数if (hexLen % 2 != 0) {return NULL;  // 奇数长度无效}int wkbLen = hexLen / 2;  // 计算WKB二进制数据长度// 分配内存存储WKB二进制数据unsigned char* wkbData = (unsigned char*)malloc(wkbLen);if (!wkbData) {return NULL;  // 内存分配失败}// 将十六进制字符串转换为字节数组for (int i = 0; i < wkbLen; i++) {char hex[3] = {wkbHex[i*2], wkbHex[i*2+1], '\0'};  // 提取两个十六进制字符wkbData[i] = (unsigned char)strtol(hex, NULL, 16);  // 转换为字节值}OGRGeometryH hGeometry = NULL;  // 初始化几何对象指针OGRErr err;                     // 错误代码变量// 检查是否是EWKB格式（至少需要9字节：1字节序+4字节类型+4字节SRID）if (wkbLen >= 9) {uint32_t geomType;  // 几何类型变量// 根据字节序读取几何类型if (wkbData[0] == 1) { // 小端序（LSB）geomType = wkbData[1] | (wkbData[2] << 8) | (wkbData[3] << 16) | (wkbData[4] << 24);} else { // 大端序（MSB）geomType = (wkbData[1] << 24) | (wkbData[2] << 16) | (wkbData[3] << 8) | wkbData[4];}// 检查是否包含SRID信息（EWKB格式标志位）if (geomType & 0x20000000) {  // 包含SRID的EWKB格式// 创建不包含SRID的标准WKB数据unsigned char* standardWkb = (unsigned char*)malloc(wkbLen - 4);if (standardWkb) {// 复制字节序标识standardWkb[0] = wkbData[0];// 复制几何类型（移除SRID标志位）uint32_t cleanGeomType = geomType & ~0x20000000;if (wkbData[0] == 1) { // 小端序写入standardWkb[1] = cleanGeomType & 0xFF;standardWkb[2] = (cleanGeomType >> 8) & 0xFF;standardWkb[3] = (cleanGeomType >> 16) & 0xFF;standardWkb[4] = (cleanGeomType >> 24) & 0xFF;} else { // 大端序写入standardWkb[1] = (cleanGeomType >> 24) & 0xFF;standardWkb[2] = (cleanGeomType >> 16) & 0xFF;standardWkb[3] = (cleanGeomType >> 8) & 0xFF;standardWkb[4] = cleanGeomType & 0xFF;}// 跳过SRID（4字节），复制剩余几何数据memcpy(standardWkb + 5, wkbData + 9, wkbLen - 9);// 从标准WKB创建几何对象err = OGR_G_CreateFromWkb(standardWkb, NULL, &hGeometry, wkbLen - 4);free(standardWkb);  // 释放临时缓冲区}} else {// 标准WKB格式，直接解析err = OGR_G_CreateFromWkb(wkbData, NULL, &hGeometry, wkbLen);}} else {// 数据长度不足，尝试直接解析（可能是简单几何）err = OGR_G_CreateFromWkb(wkbData, NULL, &hGeometry, wkbLen);}free(wkbData);  // 释放WKB数据缓冲区// 检查几何对象创建是否成功if (err != OGRERR_NONE) {if (hGeometry) {OGR_G_DestroyGeometry(hGeometry);  // 失败时清理已分配的几何对象}return NULL;  // 返回NULL表示创建失败}return hGeometry;  // 返回成功创建的几何对象
}// 从WKB数据中提取几何类型信息
uint32_t getGeometryTypeFromWKBData(const unsigned char* wkbData, int wkbLen) {// 验证输入参数（至少需要5字节：1字节序+4字节类型）if (!wkbData || wkbLen < 5) {return 0;  // 无效输入返回0}// 读取字节序标识（第一个字节）unsigned char byteOrder = wkbData[0];// 根据字节序读取几何类型（接下来4个字节）uint32_t geomType;if (byteOrder == 1) { // 小端序（LSB first）geomType = wkbData[1] | (wkbData[2] << 8) | (wkbData[3] << 16) | (wkbData[4] << 24);} else { // 大端序（MSB first）geomType = (wkbData[1] << 24) | (wkbData[2] << 16) | (wkbData[3] << 8) | wkbData[4];}// 移除SRID标志位，返回纯几何类型return geomType & ~0x20000000;
}// 清理GDAL资源的辅助函数
void cleanupGDAL() {GDALDestroyDriverManager();  // 清理GDAL驱动管理器OGRCleanupAll();            // 清理所有OGR资源
}void goErrorHandler(CPLErr eErrClass, int err_no, const char *msg) {// 可以在这里处理GDAL错误
}
*/
import "C"import ("errors" // 用于创建错误对象"log""os"      // 用于环境变量操作"runtime" // 用于运行时信息获取"unsafe"  // 用于C指针操作
)// SpatialReference 表示空间参考系统的Go包装器
type SpatialReference struct {cPtr C.OGRSpatialReferenceH // C语言空间参考系统指针
}// Geometry 表示几何对象的Go包装器
type Geometry struct {cPtr C.OGRGeometryH // C语言几何对象指针
}// InitializeGDAL 初始化GDAL库，支持自定义配置// shapeEncoding: Shapefile编码，为空时使用默认编码
func InitializeGDAL() error {// 如果未指定PROJ路径，根据操作系统设置默认路径projDataPath := ""if runtime.GOOS == "windows" {// Windows下优先使用环境变量，否则使用默认路径if envPath := os.Getenv("PROJ_DATA"); envPath != "" {projDataPath = envPath} else {projDataPath = "C:/OSGeo4W/share/proj" // Windows默认路径}} else if runtime.GOOS == "linux" {// Linux下优先使用环境变量，否则使用系统路径if envPath := os.Getenv("PROJ_DATA"); envPath != "" {projDataPath = envPath} else {projDataPath = "/usr/share/proj" // Linux默认路径}} else if runtime.GOOS == "android" {if envPath := os.Getenv("PROJ_DATA"); envPath != "" {projDataPath = envPath} else {projDataPath = "/data/data/com.termux/files/usr/share/proj" // Linux默认路径}}// 如果未指定编码，使用默认编码shapeEncoding := "GBK"// 转换Go字符串为C字符串cProjPath := C.CString(projDataPath)cEncoding := C.CString(shapeEncoding)// 确保C字符串资源被释放defer C.free(unsafe.Pointer(cProjPath))defer C.free(unsafe.Pointer(cEncoding))// 调用C函数初始化GDALC.initializeGDALWithProjFix(cProjPath, cEncoding)return nil // 初始化成功
}// 返回空间参考系统对象和可能的错误
func CreateEPSG4490WithCorrectAxis() (*SpatialReference, error) {// 调用C函数创建空间参考系统cPtr := C.createEPSG4490WithCorrectAxis()if cPtr == nil {return nil, errors.New("failed to create EPSG:4490 spatial reference system") // 创建失败}// 创建Go包装器对象srs := &SpatialReference{cPtr: cPtr}// 设置终结器，确保C资源被正确释放runtime.SetFinalizer(srs, (*SpatialReference).destroy)return srs, nil // 返回成功创建的对象
}// CreateGeometryFromWKBHex 从十六进制WKB字符串创建几何对象
// wkbHex: 十六进制格式的WKB字符串
// 返回几何对象和可能的错误
func CreateGeometryFromWKBHex(wkbHex string) (*Geometry, error) {// 验证输入参数if wkbHex == "" {return nil, errors.New("WKB hex string cannot be empty") // 空字符串错误}// 转换Go字符串为C字符串cWkbHex := C.CString(wkbHex)defer C.free(unsafe.Pointer(cWkbHex)) // 确保C字符串被释放// 调用C函数创建几何对象cPtr := C.createGeometryFromWKBHex(cWkbHex)if cPtr == nil {return nil, errors.New("failed to create geometry from WKB hex string") // 创建失败}// 创建Go包装器对象geom := &Geometry{cPtr: cPtr}return geom, nil // 返回成功创建的对象
}// GetGeometryTypeFromWKBData 从WKB数据中提取几何类型
// wkbData: WKB二进制数据
// 返回几何类型代码
func GetGeometryTypeFromWKBData(wkbData []byte) uint32 {// 验证输入数据长度if len(wkbData) < 5 {return 0 // 数据长度不足}// 调用C函数提取几何类型return uint32(C.getGeometryTypeFromWKBData((*C.uchar)(&wkbData[0]), C.int(len(wkbData))))
}// CleanupGDAL 清理GDAL资源，程序退出前调用
func CleanupGDAL() {C.cleanupGDAL() // 调用C函数清理资源
}// destroy 销毁空间参考系统的C资源（终结器函数）
func (srs *SpatialReference) destroy() {if srs.cPtr != nil {C.OSRDestroySpatialReference(srs.cPtr) // 释放C资源srs.cPtr = nil                         // 避免重复释放}
}// destroy 销毁几何对象的C资源（终结器函数）
func (geom *Geometry) destroy() {if geom.cPtr != nil {// 添加错误恢复机制defer func() {if r := recover(); r != nil {log.Printf("警告: 销毁几何对象时发生错误: %v", r)}geom.cPtr = nil}()C.OGR_G_DestroyGeometry(geom.cPtr)geom.cPtr = nil}
}// Destroy 手动销毁空间参考系统资源
func (srs *SpatialReference) Destroy() {runtime.SetFinalizer(srs, nil) // 移除终结器srs.destroy()                  // 立即释放资源
}// Destroy 手动销毁几何对象资源
func (geom *Geometry) Destroy() {runtime.SetFinalizer(geom, nil) // 移除终结器geom.destroy()                  // 立即释放资源
}

go语言中是能够直接使用原生的C语言的，写法是在通过/*/*符号来定义C语言的代码块，通过#include来引入gdal的环境变量，并写一些基础函数，在go语言的代码块中使用C.xxx 来调用相关的C语言的方法。需要注意的数据的共享，C语言的GO语言的数据结构并不互通，字符串需要使用C.CString 来进行转换，还有一个很大的坑就是，一定要在go语言的函数中手动释放C语言的相关对象，因为go的GC内存回收机制只会回收go相关的变量，如果不回收C相关的，就会造成内存泄漏，一不小心内存就直接干满了。

二、读取数据

这一步非常关键，需要将GIS数据无损转换为Gdal的OGRLayerH数据，比如shp、geodatabase、postgis等数据，Gdal是有接口能直接读取的，这里我将其封装为了函数来实现数据读取功能。以下面读取pg数据表为例：

package OSGEO/*
#include <gdal.h>
#include <gdal_alg.h>
#include <ogr_api.h>
#include <ogr_srs_api.h>
#include <cpl_error.h>
#include <cpl_conv.h>
#include <cpl_string.h>
#include <stdlib.h>
*/
import "C"import (config2 "awesomeProject/config""fmt""runtime""unsafe"
)// PostGISConfig PostGIS连接配置
type PostGISConfig struct {Host     stringPort     stringDatabase stringUser     stringPassword stringSchema   stringTable    string
}// GDALLayer 包装GDAL图层
type GDALLayer struct {layer   C.OGRLayerHdataset C.OGRDataSourceHdriver  C.OGRSFDriverH
}// PostGISReader PostGIS读取器
type PostGISReader struct {config *PostGISConfig
}// NewPostGISReader 创建新的PostGIS读取器
func NewPostGISReader(config *PostGISConfig) *PostGISReader {return &PostGISReader{config: config,}
}// ReadGeometryTable 读取PostGIS几何表数据
func (r *PostGISReader) ReadGeometryTable() (*GDALLayer, error) {// 初始化GDALInitializeGDAL()// 构建连接字符串connStr := fmt.Sprintf("PG:host=%s port=%s dbname=%s user=%s password=%s",r.config.Host, r.config.Port, r.config.Database,r.config.User, r.config.Password)cConnStr := C.CString(connStr)defer C.free(unsafe.Pointer(cConnStr))// 获取PostgreSQL驱动driver := C.OGRGetDriverByName(C.CString("PostgreSQL"))if driver == nil {return nil, fmt.Errorf("无法获取PostgreSQL驱动")}// 打开数据源dataset := C.OGROpen(cConnStr, C.int(0), nil) // 0表示只读if dataset == nil {return nil, fmt.Errorf("无法连接到PostGIS数据库: %s", connStr)}// 构建图层名称（包含schema）var layerName stringif r.config.Schema != "" {layerName = fmt.Sprintf("%s.%s", r.config.Schema, r.config.Table)} else {layerName = r.config.Table}cLayerName := C.CString(layerName)defer C.free(unsafe.Pointer(cLayerName))// 获取图层layer := C.OGR_DS_GetLayerByName(dataset, cLayerName)if layer == nil {C.OGR_DS_Destroy(dataset)return nil, fmt.Errorf("无法找到图层: %s", layerName)}gdalLayer := &GDALLayer{layer:   layer,dataset: dataset,driver:  driver,}// 设置finalizer以确保资源清理runtime.SetFinalizer(gdalLayer, (*GDALLayer).cleanup)return gdalLayer, nil
}// GetFeatureCount 获取要素数量
func (gl *GDALLayer) GetFeatureCount() int {return int(C.OGR_L_GetFeatureCount(gl.layer, C.int(1))) // 1表示强制计算
}// GetLayerDefn 获取图层定义
func (gl *GDALLayer) GetLayerDefn() C.OGRFeatureDefnH {return C.OGR_L_GetLayerDefn(gl.layer)
}// GetFieldCount 获取字段数量
func (gl *GDALLayer) GetFieldCount() int {defn := gl.GetLayerDefn()return int(C.OGR_FD_GetFieldCount(defn))
}// GetFieldName 获取字段名称
func (gl *GDALLayer) GetFieldName(index int) string {defn := gl.GetLayerDefn()fieldDefn := C.OGR_FD_GetFieldDefn(defn, C.int(index))if fieldDefn == nil {return ""}return C.GoString(C.OGR_Fld_GetNameRef(fieldDefn))
}// GetGeometryType 获取几何类型
func (gl *GDALLayer) GetGeometryType() string {defn := gl.GetLayerDefn()geomType := C.OGR_FD_GetGeomType(defn)return C.GoString(C.OGRGeometryTypeToName(geomType))
}// GetSpatialRef 获取空间参考系统
func (gl *GDALLayer) GetSpatialRef() C.OGRSpatialReferenceH {return C.OGR_L_GetSpatialRef(gl.layer)
}// ResetReading 重置读取位置
func (gl *GDALLayer) ResetReading() {C.OGR_L_ResetReading(gl.layer)
}// GetNextFeature 获取下一个要素
func (gl *GDALLayer) GetNextFeature() C.OGRFeatureH {return C.OGR_L_GetNextFeature(gl.layer)
}// PrintLayerInfo 打印图层信息
func (gl *GDALLayer) PrintLayerInfo() {fmt.Printf("图层信息:\n")fmt.Printf("  要素数量: %d\n", gl.GetFeatureCount())fmt.Printf("  几何类型: %s\n", gl.GetGeometryType())fmt.Printf("  字段数量: %d\n", gl.GetFieldCount())fmt.Printf("  字段列表:\n")for i := 0; i < gl.GetFieldCount(); i++ {fmt.Printf("    %d: %s\n", i, gl.GetFieldName(i))}// 打印空间参考系统信息srs := gl.GetSpatialRef()if srs != nil {var projStr *C.charC.OSRExportToProj4(srs, &projStr)if projStr != nil {fmt.Printf("  投影: %s\n", C.GoString(projStr))}}
}// IterateFeatures 遍历所有要素
func (gl *GDALLayer) IterateFeatures(callback func(feature C.OGRFeatureH)) {gl.ResetReading()for {feature := gl.GetNextFeature()if feature == nil {break}callback(feature)// 释放要素C.OGR_F_Destroy(feature)}
}// cleanup 清理资源
func (gl *GDALLayer) cleanup() {if gl.dataset != nil {C.OGR_DS_Destroy(gl.dataset)gl.dataset = nil}
}// Close 手动关闭资源
func (gl *GDALLayer) Close() {gl.cleanup()runtime.SetFinalizer(gl, nil)
}func MakePGReader(table string) *PostGISReader {con := config2.MainConfigconfig := &PostGISConfig{Host:     con.Host,Port:     con.Port,Database: con.Dbname,User:     con.Username,Password: con.Password,Schema:   "public", // 可选，默认为publicTable:    table,}// 创建读取器reader := NewPostGISReader(config)return reader
}

我还写了一些其他的数据格式转换函数，比如gdb转geojson，shp和geojson互转、dxf、kml和geojson互转、gdb和shp转Layer，gdb和shp直接转pg数据表等，这里就不贴代码，主要是太多了几千行，而且分的很散，等后面git仓库建好后大家自行下载即可。

三、执行空间分析

目前只完成了相交分析，就是两个gdal图层执行相交分析，根据参数来实现相交的字段类型，下面是具体流程：

package OSGEO/*
#include <gdal.h>
#include <gdal_alg.h>
#include <ogr_api.h>
#include <ogr_srs_api.h>
#include <cpl_error.h>
#include <cpl_conv.h>
#include <cpl_string.h>
#include <stdlib.h>
// 声明外部函数，避免重复定义
extern int handleProgressUpdate(double, char*, void*);// 创建内存图层用于存储相交结果
static OGRLayerH createMemoryLayer(const char* layerName, OGRwkbGeometryType geomType, OGRSpatialReferenceH srs) {// 创建内存驱动OGRSFDriverH memDriver = OGRGetDriverByName("MEM");if (!memDriver) {return NULL;}// 创建内存数据源OGRDataSourceH memDS = OGR_Dr_CreateDataSource(memDriver, "", NULL);if (!memDS) {return NULL;}// 创建图层OGRLayerH layer = OGR_DS_CreateLayer(memDS, layerName, srs, geomType, NULL);return layer;
}// 添加字段到图层
static int addFieldToLayer(OGRLayerH layer, const char* fieldName, OGRFieldType fieldType) {OGRFieldDefnH fieldDefn = OGR_Fld_Create(fieldName, fieldType);if (!fieldDefn) {return 0;}OGRErr err = OGR_L_CreateField(layer, fieldDefn, 1); // 1表示强制创建OGR_Fld_Destroy(fieldDefn);return (err == OGRERR_NONE) ? 1 : 0;
}// 复制字段值
static void copyFieldValue(OGRFeatureH srcFeature, OGRFeatureH dstFeature, int srcFieldIndex, int dstFieldIndex) {if (OGR_F_IsFieldSet(srcFeature, srcFieldIndex)) {OGRFieldDefnH fieldDefn = OGR_F_GetFieldDefnRef(srcFeature, srcFieldIndex);OGRFieldType fieldType = OGR_Fld_GetType(fieldDefn);switch (fieldType) {case OFTInteger:OGR_F_SetFieldInteger(dstFeature, dstFieldIndex, OGR_F_GetFieldAsInteger(srcFeature, srcFieldIndex));break;case OFTReal:OGR_F_SetFieldDouble(dstFeature, dstFieldIndex, OGR_F_GetFieldAsDouble(srcFeature, srcFieldIndex));break;case OFTString:OGR_F_SetFieldString(dstFeature, dstFieldIndex, OGR_F_GetFieldAsString(srcFeature, srcFieldIndex));break;default:// 其他类型转为字符串OGR_F_SetFieldString(dstFeature, dstFieldIndex, OGR_F_GetFieldAsString(srcFeature, srcFieldIndex));break;}}
}// 进度回调函数 - 这个函数会被GDAL调用
static int progressCallback(double dfComplete, const char *pszMessage, void *pProgressArg) {// pProgressArg 包含Go回调函数的信息if (pProgressArg != NULL) {// 调用Go函数处理进度更新return handleProgressUpdate(dfComplete, (char*)pszMessage, pProgressArg);}return 1; // 继续执行
}// 执行带进度监测的相交分析
static OGRErr performIntersectionWithProgress(OGRLayerH inputLayer,OGRLayerH methodLayer,OGRLayerH resultLayer,char **options,void *progressData) {return OGR_L_Intersection(inputLayer, methodLayer, resultLayer, options,progressCallback, progressData);
}*/
import "C"import ("fmt""runtime""sync""unsafe"
)// ProgressCallback 进度回调函数类型
// 返回值：true继续执行，false取消执行
type ProgressCallback func(complete float64, message string) bool// ProgressData 进度数据结构，用于在C和Go之间传递信息
type ProgressData struct {callback  ProgressCallbackcancelled boolmutex     sync.RWMutex
}// 全局进度数据映射，用于在C回调中找到对应的Go回调
var (progressDataMap   = make(map[uintptr]*ProgressData)progressDataMutex sync.RWMutex
)// handleProgressUpdate 处理来自C的进度更新
// 这个函数被C代码调用，必须使用export注释
//
//export handleProgressUpdate
func handleProgressUpdate(complete C.double, message *C.char, progressArg unsafe.Pointer) C.int {// 线程安全地获取进度数据progressDataMutex.RLock()data, exists := progressDataMap[uintptr(progressArg)]progressDataMutex.RUnlock()if !exists || data == nil {return 1 // 继续执行}// 检查是否已被取消data.mutex.RLock()if data.cancelled {data.mutex.RUnlock()return 0 // 取消执行}callback := data.callbackdata.mutex.RUnlock()if callback != nil {// 转换消息字符串msg := ""if message != nil {msg = C.GoString(message)}// 调用Go回调函数shouldContinue := callback(float64(complete), msg)if !shouldContinue {// 用户取消操作data.mutex.Lock()data.cancelled = truedata.mutex.Unlock()return 0 // 取消执行}}return 1 // 继续执行
}// IntersectionResult 相交分析结果
type IntersectionResult struct {OutputLayer *GDALLayerResultCount int
}// FieldsInfo 字段信息结构
type FieldsInfo struct {Name      stringType      C.OGRFieldTypeFromTable string // 标记字段来源表
}// FieldMergeStrategy 字段合并策略枚举
type FieldMergeStrategy intconst (// UseTable1Fields 只使用第一个表的字段UseTable1Fields FieldMergeStrategy = iota// UseTable2Fields 只使用第二个表的字段UseTable2Fields// MergePreferTable1 合并字段，冲突时优先使用table1MergePreferTable1// MergePreferTable2 合并字段，冲突时优先使用table2MergePreferTable2// MergeWithPrefix 合并字段，使用前缀区分来源MergeWithPrefix
)func (s FieldMergeStrategy) String() string {switch s {case UseTable1Fields:return "只使用表1字段"case UseTable2Fields:return "只使用表2字段"case MergePreferTable1:return "合并字段(优先表1)"case MergePreferTable2:return "合并字段(优先表2)"case MergeWithPrefix:return "合并字段(使用前缀区分)"default:return "未知策略"}
}// SpatialIntersectionAnalysis 执行带进度监控的空间相交分析
func SpatialIntersectionAnalysis(table1, table2 string, strategy FieldMergeStrategy, progressCallback ProgressCallback) (*IntersectionResult, error) {// 读取两个几何表reader1 := MakePGReader(table1)layer1, err := reader1.ReadGeometryTable()if err != nil {return nil, fmt.Errorf("读取表 %s 失败: %v", table1, err)}defer layer1.Close()reader2 := MakePGReader(table2)layer2, err := reader2.ReadGeometryTable()if err != nil {return nil, fmt.Errorf("读取表 %s 失败: %v", table2, err)}defer layer2.Close()// 执行相交分析，传递进度回调resultLayer, err := performIntersectionAnalysisWithStrategy(layer1, layer2, table1, table2, strategy, progressCallback)if err != nil {return nil, fmt.Errorf("执行相交分析失败: %v", err)}// 计算结果数量resultCount := resultLayer.GetFeatureCount()fmt.Printf("使用字段策略: %s\n", strategy.String())return &IntersectionResult{OutputLayer: resultLayer,ResultCount: resultCount,}, nil
}// performIntersectionAnalysisWithStrategy 执行相交分析的核心逻辑
func performIntersectionAnalysisWithStrategy(layer1, layer2 *GDALLayer, table1Name, table2Name string, strategy FieldMergeStrategy, progressCallback ProgressCallback) (*GDALLayer, error) {// 获取空间参考系统（使用第一个图层的SRS）srs := layer1.GetSpatialRef()// 动态确定结果几何类型（基于输入图层）defn1 := layer1.GetLayerDefn()geomType1 := C.OGR_FD_GetGeomType(defn1)defn2 := layer2.GetLayerDefn()geomType2 := C.OGR_FD_GetGeomType(defn2)// 选择更复杂的几何类型作为结果类型var resultGeomType C.OGRwkbGeometryTypeif geomType1 == C.wkbPolygon || geomType2 == C.wkbPolygon ||geomType1 == C.wkbMultiPolygon || geomType2 == C.wkbMultiPolygon {resultGeomType = C.wkbMultiPolygon} else if geomType1 == C.wkbLineString || geomType2 == C.wkbLineString ||geomType1 == C.wkbMultiLineString || geomType2 == C.wkbMultiLineString {resultGeomType = C.wkbMultiLineString} else {resultGeomType = C.wkbMultiPoint}// 创建内存图层存储结果layerName := C.CString(fmt.Sprintf("intersection_%s_%s", table1Name, table2Name))defer C.free(unsafe.Pointer(layerName))resultLayerPtr := C.createMemoryLayer(layerName, resultGeomType, srs)if resultLayerPtr == nil {return nil, fmt.Errorf("创建结果图层失败")}// 创建结果图层包装器resultLayer := &GDALLayer{layer:   resultLayerPtr,dataset: nil,driver:  nil,}runtime.SetFinalizer(resultLayer, (*GDALLayer).cleanup)// 根据策略设置结果图层的字段结构并执行相交分析err := executeIntersectionWithStrategy(layer1, layer2, resultLayer, table1Name, table2Name, strategy, progressCallback)if err != nil {return nil, fmt.Errorf("执行相交分析失败: %v", err)}// 获取结果数量resultCount := resultLayer.GetFeatureCount()fmt.Printf("相交分析完成，共生成 %d 个相交要素\n", resultCount)return resultLayer, nil
}// executeIntersectionWithStrategy 根据策略执行相交分析
func executeIntersectionWithStrategy(layer1, layer2, resultLayer *GDALLayer, table1Name, table2Name string, strategy FieldMergeStrategy, progressCallback ProgressCallback) error {var options **C.chardefer func() {if options != nil {C.CSLDestroy(options)}}()switch strategy {case UseTable1Fields:// 只使用第一个图层的字段 - 先设置结果图层字段结构err := copyLayerSchema(resultLayer, layer1)if err != nil {return fmt.Errorf("复制图层1字段结构失败: %v", err)}case UseTable2Fields:// 只使用第二个图层的字段 - 先设置结果图层字段结构err := copyLayerSchema(resultLayer, layer2)if err != nil {return fmt.Errorf("复制图层2字段结构失败: %v", err)}case MergePreferTable1:// 合并字段，冲突时优先使用table1 - 让GDAL自动处理，table1优先return executeGDALIntersectionWithProgress(layer2, layer1, resultLayer, nil, progressCallback)case MergePreferTable2:// 合并字段，冲突时优先使用table2 - 让GDAL自动处理return executeGDALIntersectionWithProgress(layer1, layer2, resultLayer, nil, progressCallback)case MergeWithPrefix:// 使用前缀区分字段来源table1Prefix := C.CString(fmt.Sprintf("INPUT_PREFIX=%s_", table1Name))table2Prefix := C.CString(fmt.Sprintf("METHOD_PREFIX=%s_", table2Name))defer C.free(unsafe.Pointer(table1Prefix))defer C.free(unsafe.Pointer(table2Prefix))options = C.CSLAddString(options, table1Prefix)options = C.CSLAddString(options, table2Prefix)return executeGDALIntersectionWithProgress(layer1, layer2, resultLayer, options, progressCallback)default:return fmt.Errorf("不支持的字段合并策略: %v", strategy)}// 对于UseTable1Fields和UseTable2Fields，需要特殊处理return executeIntersectionWithFieldFilterAndProgress(layer1, layer2, resultLayer, strategy, progressCallback)
}// executeGDALIntersectionWithProgress 带进度监测的GDAL相交分析
func executeGDALIntersectionWithProgress(inputLayer, methodLayer, resultLayer *GDALLayer, options **C.char, progressCallback ProgressCallback) error {var progressData *ProgressDatavar progressArg unsafe.Pointer// 启用多线程处理C.CPLSetConfigOption(C.CString("GDAL_NUM_THREADS"), C.CString("ALL_CPUS"))defer C.CPLSetConfigOption(C.CString("GDAL_NUM_THREADS"), nil)// 如果有进度回调，设置进度数据if progressCallback != nil {progressData = &ProgressData{callback:  progressCallback,cancelled: false,}// 将进度数据存储到全局映射中progressArg = unsafe.Pointer(progressData)progressDataMutex.Lock()progressDataMap[uintptr(progressArg)] = progressDataprogressDataMutex.Unlock()// 确保在函数结束时清理进度数据defer func() {progressDataMutex.Lock()delete(progressDataMap, uintptr(progressArg))progressDataMutex.Unlock()}()}// 调用C函数执行相交分析err := C.performIntersectionWithProgress(inputLayer.layer,methodLayer.layer,resultLayer.layer,options,progressArg,)if err != C.OGRERR_NONE {// 检查是否是用户取消导致的错误if progressData != nil {progressData.mutex.RLock()cancelled := progressData.cancelledprogressData.mutex.RUnlock()if cancelled {return fmt.Errorf("操作被用户取消")}}return fmt.Errorf("GDAL相交分析失败，错误代码: %d", int(err))}return nil
}// executeIntersectionWithFieldFilterAndProgress 执行带字段过滤和进度监测的相交分析
func executeIntersectionWithFieldFilterAndProgress(layer1, layer2, resultLayer *GDALLayer, strategy FieldMergeStrategy, progressCallback ProgressCallback) error {// 创建临时图层用于GDAL相交分析tempLayerName := C.CString("temp_intersection")defer C.free(unsafe.Pointer(tempLayerName))tempLayerPtr := C.createMemoryLayer(tempLayerName,C.OGR_FD_GetGeomType(C.OGR_L_GetLayerDefn(resultLayer.layer)),layer1.GetSpatialRef(),)if tempLayerPtr == nil {return fmt.Errorf("创建临时图层失败")}defer C.OGR_L_Dereference(tempLayerPtr)tempLayer := &GDALLayer{layer: tempLayerPtr}// 执行GDAL相交分析到临时图层（让GDAL自动处理所有字段）err := executeGDALIntersectionWithProgress(layer1, layer2, tempLayer, nil, progressCallback)if err != nil {return fmt.Errorf("临时相交分析失败: %v", err)}// 从临时图层复制要素到结果图层，只复制需要的字段return copyFeaturesWithFieldFilter(tempLayer, resultLayer, strategy)
}// copyLayerSchema 复制图层结构
func copyLayerSchema(targetLayer, sourceLayer *GDALLayer) error {sourceDefn := sourceLayer.GetLayerDefn()fieldCount := int(C.OGR_FD_GetFieldCount(sourceDefn))for i := 0; i < fieldCount; i++ {fieldDefn := C.OGR_FD_GetFieldDefn(sourceDefn, C.int(i))// 创建字段定义的副本newFieldDefn := C.OGR_Fld_Create(C.OGR_Fld_GetNameRef(fieldDefn),C.OGR_Fld_GetType(fieldDefn),)if newFieldDefn == nil {return fmt.Errorf("创建字段定义失败")}defer C.OGR_Fld_Destroy(newFieldDefn)// 复制字段属性C.OGR_Fld_SetWidth(newFieldDefn, C.OGR_Fld_GetWidth(fieldDefn))C.OGR_Fld_SetPrecision(newFieldDefn, C.OGR_Fld_GetPrecision(fieldDefn))// 添加字段到目标图层if C.OGR_L_CreateField(targetLayer.layer, newFieldDefn, 1) != C.OGRERR_NONE {fieldName := C.GoString(C.OGR_Fld_GetNameRef(fieldDefn))return fmt.Errorf("创建字段 %s 失败", fieldName)}}return nil
}// copyFeaturesWithFieldFilter 复制要素并过滤字段
func copyFeaturesWithFieldFilter(sourceLayer, targetLayer *GDALLayer, strategy FieldMergeStrategy) error {sourceLayer.ResetReading()targetDefn := targetLayer.GetLayerDefn()sourceLayer.IterateFeatures(func(sourceFeature C.OGRFeatureH) {// 创建新要素targetFeature := C.OGR_F_Create(targetDefn)if targetFeature == nil {return}defer C.OGR_F_Destroy(targetFeature)// 复制几何体geom := C.OGR_F_GetGeometryRef(sourceFeature)if geom != nil {geomClone := C.OGR_G_Clone(geom)C.OGR_F_SetGeometry(targetFeature, geomClone)C.OGR_G_DestroyGeometry(geomClone)}// 复制字段值targetFieldCount := int(C.OGR_FD_GetFieldCount(targetDefn))for i := 0; i < targetFieldCount; i++ {targetFieldDefn := C.OGR_FD_GetFieldDefn(targetDefn, C.int(i))targetFieldName := C.GoString(C.OGR_Fld_GetNameRef(targetFieldDefn))// 在源要素中查找对应的字段fieldNameCStr := C.CString(targetFieldName)sourceFieldIndex := C.OGR_F_GetFieldIndex(sourceFeature, fieldNameCStr)C.free(unsafe.Pointer(fieldNameCStr))if sourceFieldIndex >= 0 && C.OGR_F_IsFieldSet(sourceFeature, sourceFieldIndex) == 1 {// 复制字段值copyFieldValue(sourceFeature, targetFeature, int(sourceFieldIndex), i)}}// 添加要素到目标图层C.OGR_L_CreateFeature(targetLayer.layer, targetFeature)})return nil
}// copyFieldValue 复制字段值
func copyFieldValue(srcFeature, dstFeature C.OGRFeatureH, srcIndex, dstIndex int) {srcFieldDefn := C.OGR_F_GetFieldDefnRef(srcFeature, C.int(srcIndex))fieldType := C.OGR_Fld_GetType(srcFieldDefn)switch fieldType {case C.OFTInteger:value := C.OGR_F_GetFieldAsInteger(srcFeature, C.int(srcIndex))C.OGR_F_SetFieldInteger(dstFeature, C.int(dstIndex), value)case C.OFTReal:value := C.OGR_F_GetFieldAsDouble(srcFeature, C.int(srcIndex))C.OGR_F_SetFieldDouble(dstFeature, C.int(dstIndex), value)case C.OFTString:value := C.OGR_F_GetFieldAsString(srcFeature, C.int(srcIndex))C.OGR_F_SetFieldString(dstFeature, C.int(dstIndex), value)default:// 其他类型转为字符串value := C.OGR_F_GetFieldAsString(srcFeature, C.int(srcIndex))C.OGR_F_SetFieldString(dstFeature, C.int(dstIndex), value)}
}

设计了5种字段合并策略，并实现了线程安全、进度管理、以及取消机制。

我采用这版本进行测试将一个12万要素的图层和一个2.3万要素的图层进行分析（一个区县的变更调查数据和永久基本农田数据），执行时间竟然长到了3分钟，如果是用fme进行分析，26秒就能出结果

下面是fme和go的性能对比

fme中29秒完成了计算

go语言中测试竟然花费了240秒，慢了数倍，我是直接调用的原生GDAL函数OGR_L_Intersection直接对两个图层进行计算的，按理来说不应该如此慢，运行的时候cpu调用只有百分之2，可以说几乎数据都积压在了内存中。后面我看了gdal的源码终于找到了原因，因为这玩意儿在gdal内部是单线程执行的，要想提速必须进行数据分块。

四、性能提升

提升性能这地方的坑真的是太多了，最开始我想到的是通过建立空间索引来提升效率，但是折腾了一番效率提升0，因为我的数据都是在内存中，查询本来就足够快速，空间索引起到的帮助几乎可以忽略不记。那么就只剩下分块并行处理这一个办法了，但是中途又遇到一个大坑，就是你直接用go语言的goroutine来对gdal函数进行并发是会报错的，因为gdal本身内部就以及写好了多线程操作，所以对应他本身是不存在线程安全这个说法的，必须使用互斥锁才行。还有一个天坑就是数据分块后，在边界上的矢量会进行重复计算，输出的数据会出现重叠，所以预先需要对边界进行去重、同时因为并发执行，你最终的数据汇总类型、以及你的进度回调函数怎么保证进度不乱飞这都是坑，不过好歹还是折腾出来了。下面是2.0版本的执行效率：

分块+执行+数据写入到pg中一共只花费了16.5秒，其中执行只花了5秒，这个速度可以说是非常满意了。同时我将数据导入到fme中进行验证

数据重叠性验证

有一个要素出现了重叠，检查后发现是dltb本身的重叠导致的，所以去重效果是非常明显的

然后是数据包含性验证

可以看到相交分析的数据不管是和dltb 还是yjjbnt都是包含关系，说明相交是成功的。

通过inspector查看数据，属性也完美执行！

小结

做底层开发是真难啊，一不小心就内存泄漏，一不小心就是各种编译问题找不到原因，这个研究算是搞定了Go语言做强gis开发的瓶颈问题，理论上这个库完成了，可以作为当前最强性能的gis分析引擎，我将所有空间分析都做了个线程设置，只要CPU足够给力，瞬间能够把占用拉到100%，不过完成这个库的工作量也是巨大，这篇博客也是抛砖引玉，也希望能吸引到志同道合开发者的加入，有兴趣的可以后台联系，