概要
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函数初始化GDAL
C.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 string
Port string
Database string
User string
Password string
Schema string
Table string
}
// GDALLayer 包装GDAL图层
type GDALLayer struct {
layer C.OGRLayerH
dataset C.OGRDataSourceH
driver 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) {
// 初始化GDAL
InitializeGDAL()
// 构建连接字符串
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 string
if 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.char
C.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.MainConfig
config := &PostGISConfig{
Host: con.Host,
Port: con.Port,
Database: con.Dbname,
User: con.Username,
Password: con.Password,
Schema: "public", // 可选，默认为public
Table: 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 ProgressCallback
cancelled bool
mutex 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.callback
data.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 = true
data.mutex.Unlock()
return 0 // 取消执行
}
}
return 1 // 继续执行
}
// IntersectionResult 相交分析结果
type IntersectionResult struct {
OutputLayer *GDALLayer
ResultCount int
}
// FieldsInfo 字段信息结构
type FieldsInfo struct {
Name string
Type C.OGRFieldType
FromTable string // 标记字段来源表
}
// FieldMergeStrategy 字段合并策略枚举
type FieldMergeStrategy int
const (
// UseTable1Fields 只使用第一个表的字段
UseTable1Fields FieldMergeStrategy = iota
// UseTable2Fields 只使用第二个表的字段
UseTable2Fields
// MergePreferTable1 合并字段，冲突时优先使用table1
MergePreferTable1
// MergePreferTable2 合并字段，冲突时优先使用table2
MergePreferTable2
// 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.OGRwkbGeometryType
if 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.char
defer 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 *ProgressData
var 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)] = progressData
progressDataMutex.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.cancelled
progressData.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%，不过完成这个库的工作量也是巨大，这篇博客也是抛砖引玉，也希望能吸引到志同道合开发者的加入，有兴趣的可以后台联系，
扫一扫
3418
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
