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

基于Go语言开发的GIS空间分析库Gogeo使用教程

news 2025/10/25 12:34:01

    • 概要

      配置库

      使用说明

      1. 数据读取

      2. 数据分析

      3. 数据导出

      使用示例

      小结

概要

        经过一段时间的完善和优化，目前基于Go语言和Gdal开发的GIS空间分析库已经正式开源至Github。该库目前能实现shp、gdb、postgis矢量数据的写入和写出，并能实现将数据读取为GoLayer图层，并能实现相交、裁剪、融合、标识、交集取反、擦除等图层空间叠加分析。结合了Golang的语言级别的高并发能力，将图层分块并行计算，能达到对CPU的高度压榨，经过测试性能远超大部分GIS桌面产品。

        Github地址：https://github.com/GrainArc/Gogeo

配置库

       该库需要提前配置Gdal环境。

        针对windows，我们打包了快捷安装包，能一键配置Gdal的环境，无限下一步安装即可。百度云盘连接：https://pan.baidu.com/s/1tSNRXKdkZdswruyN1Fl8kg?pwd=emui

        针对Linux环境就简单了，如果是arm版本就直接pkg install gdal 就行了，如果是amd环境的直接运行apt install gdal-bin和 apt install libgdal-dev命令即可。

使用说明

       以windows系统为例，如果是直接默认安装路径， 可以无视。如果是自定义了安装路径，需要在cgo_header.go代码中配置一下cgo的路径，将里面的路径修改为你的实际安装路径

#cgo windows CFLAGS: -IC:/OSGeo4W/include -IC:/OSGeo4W/include/gdal -Wno-unused-result
#cgo windows LDFLAGS: -LC:/OSGeo4W/lib -lgdal_i -lstdc++ -static-libgcc -static-libstdc++

        安装库指令

go get -u github.com/GrainArc/Gogeo

1. 数据读取

shp读取，先初始化一个读取器，然后调用ReadShapeFile函数实现数据读取

shpFile1 := "C:\\Users\\Administrator\\Desktop\\bb\\xxx.shp"	
reader1, _ := Gogeo.NewFileGeoReader(shpFile1)
layer1, _ := reader1.ReadShapeFile()

gdb数据读取，方法同上，换成ReadGDBFile函数即可

gdbFile1 := "C:\\Users\\Administrator\\Desktop\\bb\\xxx.gdb"	
reader1, _ := Gogeo.NewFileGeoReader(gdbFile1 )
layer1, _ := reader1.ReadGDBFile()

postgis数据库读取

config := &PostGISConfig{Host:     "localhost",  // 数据库主机地址Port:     "5432",       // PostgreSQL默认端口Database: "testdb",     // 数据库名称User:     "postgres",   // 用户名Password: "password",   // 密码Schema:   "public",     // 模式名称，默认为publicTable:    "test_table", // 要读取的表名}// 创建PostGIS读取器reader := NewPostGISReader(config)if reader == nil {t.Fatal("创建PostGIS读取器失败")}// 读取几何表数据layer, err := reader.ReadGeometryTable()if err != nil {t.Fatalf("读取PostGIS表失败: %v", err)}

2. 数据分析

相交分析参数：

inputLayer：输入图层
methodLayer：相交图层
config ：分析配置参数
strategy：分析字段保留策略
	shpFile1 := "C:\\Users\\Administrator\\Desktop\\bb\\aa.shp"shpFile2 := "C:\\Users\\Administrator\\Desktop\\bb\\bb.shp"fmt.Printf("开始时间: %s\n", startTime.Format("2006-01-02 15:04:05"))fmt.Println("----------------------------------------")fmt.Println("正在读取第一个shapefile...")reader1, _ := Gogeo.NewFileGeoReader(shpFile1)layer1, _ := reader1.ReadShapeFile()reader2, _ := Gogeo.NewFileGeoReader(shpFile2)layer2, _ := reader2.ReadShapeFile()// 3. 配置并行相交分析参数config := &Gogeo.ParallelGeosConfig{TileCount:        10,               // 4x4分块MaxWorkers:       16,               // 使用所有CPU核心BufferDistance:   0.0001,           // 分块缓冲距离IsMergeTile:      true,             // 合并分块结果ProgressCallback: progressCallback, // 进度回调函数PrecisionConfig: &Gogeo.GeometryPrecisionConfig{Enabled:       true,GridSize:      0.000000001, // 容差PreserveTopo:  true,        // 保持拓扑KeepCollapsed: false,       // 不保留退化几何},}// 修改你的调用代码，捕获具体的错误信息OUT, err := Gogeo.SpatialIntersectionAnalysis(layer1, layer2, config, 2)if err != nil {fmt.Printf("空间擦除分析失败: %v\n", err)return}

其他空间分析方案使用方式都是类似，具体参数配置查看源码

3. 数据导出

GDB写出，参数说明：

sourceLayer *GDALLayer ：输出图层

filePath string ： 写出文件路径

layerName string：写出gdb中的图层名

overwrite bool ：是否覆盖，如果不覆盖会追加要素至已有gdb中

err = Gogeo.WriteGDBLayer(OutputLayer, "C:\\Users\\Administrator\\Desktop\\bb\\layerName.gdb", "layerName", true)

shp写出同上，函数替换一下即可

err = Gogeo.WriteShapeFileLayer(OUT.OutputLayer, "C:\\Users\\Administrator\\Desktop\\bb\\layerName.shp", "layerName", true)

postgis数据库写入

import ("fmt""gorm.io/driver/postgres""gorm.io/gorm""testing""time"
)// TestSaveGDALLayerToPG 测试将GDALLayer保存到PostGIS
func TestSaveGDALLayerToPG(t *testing.T) {// 1. 配置数据库连接dsn := "host=localhost user=postgres password=password dbname=testdb port=5432 sslmode=disable TimeZone=Asia/Shanghai"db, err := gorm.Open(postgres.Open(dsn), &gorm.Config{})if err != nil {t.Fatalf("连接数据库失败: %v", err)}// 2. 读取源数据（从Shapefile）shpFilePath := "./testdata/test.shp"reader, err := NewFileGeoReader(shpFilePath)if err != nil {t.Fatalf("创建文件读取器失败: %v", err)}gdalLayer, err := reader.ReadShapeFile()if err != nil {t.Fatalf("读取Shapefile失败: %v", err)}defer gdalLayer.Close()// 3. 打印源数据信息fmt.Println("\n========== 源数据信息 ==========")gdalLayer.PrintLayerSummary()// 4. 保存到PostGIStableName := "test_layer"schema := "public"srid := 4326 // WGS84坐标系fmt.Println("\n========== 开始保存到PostGIS ==========")err = SaveGDALLayerToPG(db, gdalLayer, tableName, schema, srid)if err != nil {t.Fatalf("保存到PostGIS失败: %v", err)}// 5. 验证数据是否成功写入fmt.Println("\n========== 验证写入结果 ==========")verifyResult := verifyPostGISTable(t, db, tableName, schema)if !verifyResult {t.Error("数据验证失败")}fmt.Println("\n========== 测试完成 ==========")
}// TestSaveGDALLayerToPGBatch 测试批量保存GDALLayer到PostGIS
func TestSaveGDALLayerToPGBatch(t *testing.T) {// 1. 配置数据库连接dsn := "host=localhost user=postgres password=password dbname=testdb port=5432 sslmode=disable TimeZone=Asia/Shanghai"db, err := gorm.Open(postgres.Open(dsn), &gorm.Config{})if err != nil {t.Fatalf("连接数据库失败: %v", err)}// 2. 读取源数据shpFilePath := "./testdata/large_dataset.shp"reader, err := NewFileGeoReader(shpFilePath)if err != nil {t.Fatalf("创建文件读取器失败: %v", err)}gdalLayer, err := reader.ReadShapeFile()if err != nil {t.Fatalf("读取Shapefile失败: %v", err)}defer gdalLayer.Close()// 3. 打印源数据信息fmt.Println("\n========== 源数据信息 ==========")info, err := GetGDALLayerInfo(gdalLayer)if err != nil {t.Fatalf("获取图层信息失败: %v", err)}printLayerInfo(info)// 4. 批量保存到PostGIStableName := "large_test_layer"schema := "public"srid := 4326batchSize := 500 // 每批处理500条fmt.Println("\n========== 开始批量保存到PostGIS ==========")startTime := time.Now()err = SaveGDALLayerToPGBatch(db, gdalLayer, tableName, schema, srid, batchSize)if err != nil {t.Fatalf("批量保存到PostGIS失败: %v", err)}elapsed := time.Since(startTime)fmt.Printf("批量保存耗时: %v\n", elapsed)// 5. 验证数据fmt.Println("\n========== 验证写入结果 ==========")verifyResult := verifyPostGISTable(t, db, tableName, schema)if !verifyResult {t.Error("数据验证失败")}fmt.Println("\n========== 测试完成 ==========")
}// TestSaveMultipleLayersToPG 测试保存多个图层到PostGIS
func TestSaveMultipleLayersToPG(t *testing.T) {// 1. 配置数据库连接dsn := "host=localhost user=postgres password=password dbname=testdb port=5432 sslmode=disable TimeZone=Asia/Shanghai"db, err := gorm.Open(postgres.Open(dsn), &gorm.Config{})if err != nil {t.Fatalf("连接数据库失败: %v", err)}// 2. 定义要处理的文件列表testFiles := []struct {FilePath  stringTableName stringSRID      int}{{"./testdata/points.shp", "test_points", 4326},{"./testdata/lines.shp", "test_lines", 4326},{"./testdata/polygons.shp", "test_polygons", 4326},}schema := "public"// 3. 遍历处理每个文件for _, testFile := range testFiles {t.Run(testFile.TableName, func(t *testing.T) {fmt.Printf("\n========== 处理文件: %s ==========\n", testFile.FilePath)// 读取数据reader, err := NewFileGeoReader(testFile.FilePath)if err != nil {t.Logf("跳过文件 %s: %v", testFile.FilePath, err)return}gdalLayer, err := reader.ReadShapeFile()if err != nil {t.Logf("读取文件失败 %s: %v", testFile.FilePath, err)return}defer gdalLayer.Close()// 打印信息gdalLayer.PrintLayerSummary()// 保存到PostGISerr = SaveGDALLayerToPG(db, gdalLayer, testFile.TableName, schema, testFile.SRID)if err != nil {t.Errorf("保存失败 %s: %v", testFile.TableName, err)return}// 验证if !verifyPostGISTable(t, db, testFile.TableName, schema) {t.Errorf("验证失败: %s", testFile.TableName)}fmt.Printf("✓ 成功处理: %s\n", testFile.TableName)})}
}// TestSaveGDBLayerToPG 测试从GDB读取并保存到PostGIS
func TestSaveGDBLayerToPG(t *testing.T) {// 1. 配置数据库连接dsn := "host=localhost user=postgres password=password dbname=testdb port=5432 sslmode=disable TimeZone=Asia/Shanghai"db, err := gorm.Open(postgres.Open(dsn), &gorm.Config{})if err != nil {t.Fatalf("连接数据库失败: %v", err)}// 2. 读取GDB数据gdbFilePath := "./testdata/test.gdb"layerName := "test_layer" // GDB中的图层名称reader, err := NewFileGeoReader(gdbFilePath)if err != nil {t.Fatalf("创建GDB读取器失败: %v", err)}gdalLayer, err := reader.ReadGDBFile(layerName)if err != nil {t.Fatalf("读取GDB图层失败: %v", err)}defer gdalLayer.Close()// 3. 打印源数据信息fmt.Println("\n========== GDB图层信息 ==========")gdalLayer.PrintLayerInfo()// 4. 保存到PostGIStableName := "gdb_test_layer"schema := "public"srid := 4326fmt.Println("\n========== 保存GDB数据到PostGIS ==========")err = SaveGDALLayerToPG(db, gdalLayer, tableName, schema, srid)if err != nil {t.Fatalf("保存GDB数据到PostGIS失败: %v", err)}// 5. 验证fmt.Println("\n========== 验证写入结果 ==========")verifyResult := verifyPostGISTable(t, db, tableName, schema)if !verifyResult {t.Error("数据验证失败")}
}// TestSavePostGISToPostGIS 测试从PostGIS读取并保存到另一个表
func TestSavePostGISToPostGIS(t *testing.T) {// 1. 配置数据库连接dsn := "host=localhost user=postgres password=password dbname=testdb port=5432 sslmode=disable TimeZone=Asia/Shanghai"db, err := gorm.Open(postgres.Open(dsn), &gorm.Config{})if err != nil {t.Fatalf("连接数据库失败: %v", err)}// 2. 从PostGIS读取源表sourceTable := "source_table"reader := MakePGReader(sourceTable)gdalLayer, err := reader.ReadGeometryTable()if err != nil {t.Fatalf("读取PostGIS源表失败: %v", err)}defer gdalLayer.Close()// 3. 打印源表信息fmt.Println("\n========== 源表信息 ==========")gdalLayer.PrintLayerSummary()// 4. 保存到新表targetTable := "target_table"schema := "public"srid := 4326fmt.Println("\n========== 复制到新表 ==========")err = SaveGDALLayerToPG(db, gdalLayer, targetTable, schema, srid)if err != nil {t.Fatalf("保存到新表失败: %v", err)}// 5. 验证fmt.Println("\n========== 验证目标表 ==========")verifyResult := verifyPostGISTable(t, db, targetTable, schema)if !verifyResult {t.Error("数据验证失败")}
}// TestSaveWithDifferentSRID 测试不同坐标系的转换和保存
func TestSaveWithDifferentSRID(t *testing.T) {// 1. 配置数据库连接dsn := "host=localhost user=postgres password=password dbname=testdb port=5432 sslmode=disable TimeZone=Asia/Shanghai"db, err := gorm.Open(postgres.Open(dsn), &gorm.Config{})if err != nil {t.Fatalf("连接数据库失败: %v", err)}// 2. 读取源数据shpFilePath := "./testdata/test_wgs84.shp"reader, err := NewFileGeoReader(shpFilePath)if err != nil {t.Fatalf("创建文件读取器失败: %v", err)}gdalLayer, err := reader.ReadShapeFile()if err != nil {t.Fatalf("读取Shapefile失败: %v", err)}defer gdalLayer.Close()// 3. 测试不同的坐标系testCases := []struct {TableName stringSRID      intSRIDName  string}{{"test_wgs84", 4326, "WGS84"},{"test_cgcs2000", 4490, "CGCS2000"},{"test_web_mercator", 3857, "Web Mercator"},}schema := "public"for _, tc := range testCases {t.Run(tc.SRIDName, func(t *testing.T) {fmt.Printf("\n========== 测试坐标系: %s (EPSG:%d) ==========\n", tc.SRIDName, tc.SRID)err = SaveGDALLayerToPG(db, gdalLayer, tc.TableName, schema, tc.SRID)if err != nil {t.Errorf("保存失败 (SRID %d): %v", tc.SRID, err)return}// 验证坐标系if !verifySRID(t, db, tc.TableName, schema, tc.SRID) {t.Errorf("坐标系验证失败: %s", tc.SRIDName)}fmt.Printf("✓ 成功保存并验证: %s\n", tc.SRIDName)})}
}// verifyPostGISTable 验证PostGIS表中的数据
func verifyPostGISTable(t *testing.T, db *gorm.DB, tableName, schema string) bool {sqlDB, err := db.DB()if err != nil {t.Errorf("获取数据库连接失败: %v", err)return false}// 1. 检查表是否存在var exists boolquery := `SELECT EXISTS (SELECT 1 FROM information_schema.tables WHERE table_schema = $1 AND table_name = $2)`err = sqlDB.QueryRow(query, schema, tableName).Scan(&exists)if err != nil {t.Errorf("检查表存在性失败: %v", err)return false}if !exists {t.Errorf("表不存在: %s.%s", schema, tableName)return false}fmt.Printf("✓ 表存在: %s.%s\n", schema, tableName)// 2. 检查记录数var count intcountQuery := fmt.Sprintf("SELECT COUNT(*) FROM %s.%s", schema, tableName)err = sqlDB.QueryRow(countQuery).Scan(&count)if err != nil {t.Errorf("查询记录数失败: %v", err)return false}fmt.Printf("✓ 记录数: %d\n", count)// 3. 检查几何字段var geomCount intgeomQuery := fmt.Sprintf("SELECT COUNT(*) FROM %s.%s WHERE geom IS NOT NULL", schema, tableName)err = sqlDB.QueryRow(geomQuery).Scan(&geomCount)if err != nil {t.Errorf("查询几何字段失败: %v", err)return false}fmt.Printf("✓ 有效几何数: %d\n", geomCount)// 4. 检查几何类型var geomType stringtypeQuery := `SELECT type FROM geometry_columns WHERE f_table_schema = $1 AND f_table_name = $2`err = sqlDB.QueryRow(typeQuery, schema, tableName).Scan(&geomType)if err != nil {t.Logf("查询几何类型失败: %v", err)} else {fmt.Printf("✓ 几何类型: %s\n", geomType)}// 5. 检查空间索引var indexExists boolindexQuery := `SELECT EXISTS (SELECT 1 FROM pg_indexes WHERE schemaname = $1 AND tablename = $2 AND indexname LIKE '%geom%')`err = sqlDB.QueryRow(indexQuery, schema, tableName).Scan(&indexExists)if err != nil {t.Logf("检查空间索引失败: %v", err)} else {if indexExists {fmt.Println("✓ 空间索引存在")} else {fmt.Println("⚠ 空间索引不存在")}}return true
}// verifySRID 验证表的坐标系
func verifySRID(t *testing.T, db *gorm.DB, tableName, schema string, expectedSRID int) bool {sqlDB, err := db.DB()if err != nil {t.Errorf("获取数据库连接失败: %v", err)return false}var srid intquery := `SELECT srid FROM geometry_columns WHERE f_table_schema = $1 AND f_table_name = $2`err = sqlDB.QueryRow(query, schema, tableName).Scan(&srid)if err != nil {t.Errorf("查询SRID失败: %v", err)return false}if srid != expectedSRID {t.Errorf("SRID不匹配: 期望 %d, 实际 %d", expectedSRID, srid)return false}fmt.Printf("✓ SRID验证通过: %d\n", srid)return true
}// printLayerInfo 打印图层信息
func printLayerInfo(info map[string]interface{}) {fmt.Printf("图层名称: %v\n", info["layer_name"])fmt.Printf("要素数量: %v\n", info["feature_count"])fmt.Printf("几何类型: %v\n", info["geometry_type"])fmt.Printf("字段数量: %v\n", info["field_count"])if fields, ok := info["fields"].([]map[string]interface{}); ok {fmt.Println("\n字段列表:")for i, field := range fields {fmt.Printf("  %d. %v (%v)\n", i+1, field["name"], field["type"])}}if srsAuth, ok := info["srs_authority"]; ok {if srsCode, ok := info["srs_code"]; ok {fmt.Printf("坐标系: %v:%v\n", srsAuth, srsCode)}}
}// BenchmarkSaveGDALLayerToPG 性能测试
func BenchmarkSaveGDALLayerToPG(b *testing.B) {// 配置数据库连接dsn := "host=localhost user=postgres password=password dbname=testdb port=5432 sslmode=disable"db, err := gorm.Open(postgres.Open(dsn), &gorm.Config{})if err != nil {b.Fatalf("连接数据库失败: %v", err)}// 读取测试数据shpFilePath := "./testdata/benchmark.shp"reader, err := NewFileGeoReader(shpFilePath)if err != nil {b.Fatalf("创建文件读取器失败: %v", err)}gdalLayer, err := reader.ReadShapeFile()if err != nil {b.Fatalf("读取Shapefile失败: %v", err)}defer gdalLayer.Close()b.ResetTimer()for i := 0; i < b.N; i++ {tableName := fmt.Sprintf("benchmark_table_%d", i)err = SaveGDALLayerToPGBatch(db, gdalLayer, tableName, "public", 4326, 1000)if err != nil {b.Fatalf("保存失败: %v", err)}}
}

使用示例

        基于该库，我们开发了一套web端的gis空间分析模块，通过websocket进行数据传输，能实时更新分析进度，下图将加载的两个图层矢量瓦片服务进行空间分析，分析完成后，会将输出图层添加在inputlayer的目录下。

分析完成图层已经添加进目录，并且已经完成矢量瓦片服务的发布，总耗时22秒，要素数量10万，速度已经超越C端的大部分GIS产品。

小结

        因为gdal的空间分析本身就不是线程安全的，所以为了保证并发的稳定性，空间分析采用了分块并行机制，该并行机制，实现了对layer图层进行快速瓦片分割，并将分割文件保存为二进制序列化文件，再使用序列化器将临时二进制文件并发读取并序列化为图层进行并发计算，并且计算完成后按照边界索引检测，实现边界分割要素的融合。该算法调试了N个版本，现在测试完全可以平替大部分空间分析功能。当然该库目前也不算完善，最为一个合格的GIS库，还欠缺一打功能（栅格工具、其他矢量格式的兼容等），这篇文也希望能吸引到有兴趣的开发者能够加入，一起来完成该库的开发。

http://www.dtcms.com/a/525427.html

相关文章：

  • 10.24快乐！
  • 成都商报官方网站做酒业网站的要求
  • 定制网站建设公司费用中国有几大建设
  • 在 MS Excel 和 Google Sheets 中生成 3CX 可视化通话报告
  • pfc性能优化_占空比前馈优化
  • 北京大兴专业网站建设公司低价网站建设为您公司省去了什么
  • 做网站输入文本框做下拉网页页面设计叫什么
  • 【开题答辩全过程】以 “塞上江南”旅游网站为例，包含答辩的问题和答案
  • wordpress 建站教程网站建设ppt答辩
  • 特价网站建设费用拼团网站开发
  • 门户网站制作哪专业c语言网站
  • Linux小课堂: Git与版本控制之技术演进、核心原理与企业级实践
  • Spring AI Alibaba 【二】
  • eSIM上线，是全面进化还是开倒车
  • Spring AI Alibaba 10分钟快速入门
  • 做网站被罚款门户网站建设哪专业
  • 电子电气架构车载网关系列——常见网关芯片特点
  • trae ide 设置 terminal 使用 powershell , 默认加载 用户和系统环境变量
  • 有哪些摄影网站电子商务网站建设结业论文
  • 网站建设 应该考虑什么wordpress 友情链接调用
  • 嵌入式开发 | C语言 | 单精度浮点数4字节可以表示的范围计算过程
  • JMeter测试HTTP POST（附实例）
  • 网站建设 页面网站ip地址查询
  • UM681A相关参数性能介绍
  • 建企业网站多少钱官网cms
  • html制作音乐网站Erphpdown wordpress
  • 最快做网站的语言wordpress 函数api文件
  • php企业网站开发实验总结软件开发工具包sdk
  • java发送SOAP请求
  • Web单页应用（SPA）路由设计（以React为例）