????????概要

概要

????????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%，不過完成這個庫的工作量也是巨大，這篇博客也是拋磚引玉，也希望能吸引到志同道合開發者的加入，有興趣的可以后臺聯系，