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

WebGIS开发之地形土方开挖回填分析

news 来源：原创 2025/4/29 9:21:43

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

文章目录

  • 目录

    文章目录

    前言

    一、获取真实Dem的高程数据

    二、设计面高程点处理

    1.设计面边界插值

    2.三角网转换

    3.获取某点在设计面上的投影高程

    二、设计面进行微元

    三、填挖方计算

    四、实现效果展示

    总结

前言

        填挖方分析是GIS领域中非常常规的功能，不少webgis框架都有现成的功能，比如cesnium、DTS等，但是实际使用起来效果不佳，因为大部分这类算法的设计面都是统一高程，不满足实际需求。本博客将实现从底层算法开始完成这一计算过程。底层逻辑是将设计面和现状面的分割成无数能够方便计算柱状体，然后进行数据累加计算出填方和挖方。涉及技术：golang、postgis。前端通过传入高程点和几何面，后端直接计算出填挖方数据。

一、获取真实Dem的高程数据

在WebGis中，三维地形数据的发布有很多规范，本研究主要参考Mapbox的地形发布规范，即将Dem数据通过数学公式将高程值分配到RGB三个通道的数据中（因为常规图片RGB只能支持0-255之间的值，高程值通常数值较大，需要这样处理），具体计算公式如下：

height = -10000 + ((R * 256 * 256 + G * 256 + B) * 0.1)

但是获取某个点在发布瓦片中的具体高程，需要解决几个问题：

1、通过某点的坐标查询到瓦片的z，x，y序列号。

2、对应瓦片的绝对坐标转换为瓦片上的相对坐标

3、提取高程值

实现代码如下：

//瓦片转经纬度
func TileToLatLon(z, x, y int) (topLeft, topRight, bottomLeft, bottomRight [2]float64) {tileSize := 256// 计算瓦片四个角的像素坐标topLeftPx := x * tileSizetopLeftPy := y * tileSizebottomRightPx := (x + 1) * tileSizebottomRightPy := (y + 1) * tileSize// 计算并返回四个角的经纬度topLeft = pixelToLatLon(topLeftPx, topLeftPy, z)topRight = pixelToLatLon(bottomRightPx, topLeftPy, z)bottomLeft = pixelToLatLon(topLeftPx, bottomRightPy, z)bottomRight = pixelToLatLon(bottomRightPx, bottomRightPy, z)return topLeft, topRight, bottomLeft, bottomRight
}
func pixelToLatLon(px, py, z int) [2]float64 {// 将像素坐标转换为经纬度tileSize := 256mapSize := tileSize * int(math.Pow(2, float64(z)))lonDeg := float64(px)/float64(mapSize)*360.0 - 180.0latRad := math.Atan(math.Sinh(math.Pi * (1 - 2*float64(py)/float64(mapSize))))latDeg := latRad * 180.0 / math.Pivar data [2]float64data[0] = lonDegdata[1] = latDegreturn data
}
//查询瓦片号
func LonLatToTile(lon, lat float64, zoom int64) (x, y int64) {// 1. 先转换为Web墨卡托坐标const (EarthRadius = 6378137.0OriginShift = 2 * math.Pi * EarthRadius / 2.0)mercX := lon * OriginShift / 180.0mercY := math.Log(math.Tan((90+lat)*math.Pi/360.0)) * OriginShift / math.Pi// 2. 计算瓦片坐标resolution := (2 * OriginShift) / math.Exp2(float64(zoom))x = int64(math.Floor((mercX + OriginShift) / resolution))y = int64(math.Floor((OriginShift - mercY) / resolution))// 3. 处理边界情况maxTile := int64(math.Exp2(float64(zoom))) - 1if x < 0 {x = 0} else if x > maxTile {x = maxTile}if y < 0 {y = 0} else if y > maxTile {y = maxTile}return
}
func GetRGBToDEM(lon float64, lat float64, ResultTiles []models.Tile, maxZoom int64) float64 {x1, y1 := pgmvt.LonLatToTile(lon, lat, maxZoom)var realDem models.Tilefor _, item := range ResultTiles {if item.TileRow == y1 && item.TileColumn == x1 {realDem = item}}// 自动检测图片格式img, format, err := decodeImage(realDem.TileData)if err != nil {log.Printf("图片解码失败（格式：%s）：%v", format, err)return 0}// 计算lat，lon在该照片上的相对坐标x,y//计算所在瓦片的经纬度坐标topLeft, _, bottomLeft, bottomRight := pgmvt.TileToLatLon(int(maxZoom), int(realDem.TileColumn), int(realDem.TileRow))// 计算瓦片经纬度范围tileLeft := topLeft[0]tileRight := bottomRight[0]tileTop := topLeft[1]tileBottom := bottomLeft[1]lonRatio := (lon - tileLeft) / (tileRight - tileLeft)latRatio := (tileTop - lat) / (tileTop - tileBottom)// 获取图片实际尺寸bounds := img.Bounds()width := bounds.Max.X - bounds.Min.Xheight := bounds.Max.Y - bounds.Min.Y// 计算像素坐标（考虑边界溢出）x := int(math.Floor(float64(width) * lonRatio))y := int(math.Floor(float64(height) * latRatio))// 检查坐标范围bounds2 := img.Bounds()if x < bounds2.Min.X || x >= bounds2.Max.X || y < bounds2.Min.Y || y >= bounds2.Max.Y {return 0}// 获取颜色值（统一处理）r, g, b := getRGB(img.At(x, y))// Mapbox DEM高程计算公式（优化版）height2 := calculateElevation(r, g, b)return height2
}// 自动检测并解码图片
func decodeImage(data []byte) (image.Image, string, error) {// 先尝试WebP解码（Mapbox常用）if img, err := webp.Decode(bytes.NewReader(data)); err == nil {return img, "webp", nil}// 再尝试PNG解码if img, format, err := image.Decode(bytes.NewReader(data)); err == nil {return img, format, nil}return nil, "unknown", fmt.Errorf("无法识别的图片格式")
}// 统一获取RGB值
func getRGB(c color.Color) (r, g, b uint8) {switch color := c.(type) {case color.NRGBA:return color.R, color.G, color.Bcase color.RGBA:return color.R, color.G, color.Bdefault:r32, g32, b32, _ := c.RGBA()return uint8(r32 >> 8), uint8(g32 >> 8), uint8(b32 >> 8)}
}// 高程计算公式（Mapbox官方算法）
func calculateElevation(r, g, b uint8) float64 {// 公式：height = (R * 256² + G * 256 + B) * 0.1 - 10000return (float64(r)*65536+float64(g)*256+float64(b))*0.1 - 10000
}

通过以上算法，能实现输入点经纬度坐标，返回其在地形服务中的高程值

二、设计面高程点处理

前端实现传入设计面高程点和设计面本身，将其转换为三角网，这部分主要使用postgis函数实现

1.设计面边界插值

采用临近插值法，主要通过ST_InterpolatePoint函数实现查询最近点的高程

interpolated_vertices AS (SELECT ST_SetSRID(ST_MakePoint(ST_X(v.vertex),ST_Y(v.vertex),COALESCE((SELECT ST_Z(p.geom) FROM user_points p ORDER BY v.vertex <-> ST_SetSRID(ST_MakePoint(ST_X(p.geom), ST_Y(p.geom)),4326)LIMIT 1),ST_InterpolatePoint(ST_Collect(p.geom), v.vertex))),4326 )AS geomFROM surface_vertices vLEFT JOIN user_points p ON trueGROUP BY v.vertex
)

2.三角网转换

代码如下，将面分解为点，然后和高程点组合构造三角网，主要通过ST_DelaunayTriangles函数实现： 

combined_points AS (SELECT geom FROM interpolated_verticesUNION ALLSELECT geom FROM user_points
),
tin AS (SELECT ST_DelaunayTriangles(ST_Union(geom)) AS geom FROM combined_points
)

3.获取某点在设计面上的投影高程

代码如下，主要功能包括：

  1. 通过子查询point_elevations计算每个输入点的高程

  2. 使用COALESCE函数尝试两种高程计算方法：

    • 首先尝试找到包含该点的三角形平面，使用平面方程计算高程

    • 如果找不到包含点或法向量为或法向量为或法向量为0，则使用ST_3DClosestPoint计算最近点的高程

  3. 最终输出所有点的坐标数组[x,y,z]

point_elevations AS (SELECT ip.id,ip.geom,COALESCE((SELECT (ax * (ST_X(ip.geom) - bx) + ay * (ST_Y(ip.geom) - by) + az) / normal_z FROM (SELECT /*+ Parallel(t 4) */v1.x AS ax, v1.y AS ay, v1.z AS az,v2.x AS bx, v2.y AS by, v2.z AS bz,v3.x AS cx, v3.y AS cy, v3.z AS cz,( (v2.y - v1.y)*(v3.z - v1.z) - (v2.z - v1.z)*(v3.y - v1.y) ) AS normal_zFROM (SELECT (ST_Dump(geom)).path[1] AS tri_id,(ST_DumpPoints((ST_Dump(geom)).geom)).geom AS vertex FROM tin) tLEFT JOIN LATERAL (SELECT ST_X(vertex)::numeric AS x, ST_Y(vertex)::numeric AS y, ST_Z(vertex)::numeric AS z) v1 ON trueLEFT JOIN LATERAL (SELECT ST_X(vertex)::numeric AS x, ST_Y(vertex)::numeric AS y, ST_Z(vertex)::numeric AS z) v2 ON trueLEFT JOIN LATERAL (SELECT ST_X(vertex)::numeric AS x, ST_Y(vertex)::numeric AS y, ST_Z(vertex)::numeric AS z) v3 ON trueWHERE tri_id IS NOT NULL) AS planeWHERE ST_Contains(ST_MakePolygon(ST_MakeLine(ARRAY[ST_SetSRID(ST_MakePoint(ax, ay), 4326),ST_SetSRID(ST_MakePoint(bx, by), 4326),ST_SetSRID(ST_MakePoint(cx, cy), 4326)])),ip.geom)AND normal_z <> 0LIMIT 1),ST_Z(ST_3DClosestPoint(tin.geom, ip.geom))) AS zFROM tinCROSS JOIN LATERAL (SELECT /*+ Parallel(ip 4) */id, geom FROM input_points ip) ip
)
SELECT array_agg(ARRAY[ST_X(geom), ST_Y(geom), z]) AS coordinates
FROM point_elevations;

二、设计面进行微元

        说人话就是将设计面tiler，分解成无数小面，取小面的中心点来获取设计面和现状面的高程差，再计算每个小面和高程差的值，并按正负进行累加，就得到了填方和挖方。该计算方法对比方格网法会更加简单粗暴，理论小tiler的数值越小，计算值就越接近真实值，但是性能损耗会越大，所以针对面积小于7000的地块采用1米的小方格，面积大于7000的采用自适应切块。实现代码如下：

func GetTilesPoints(geo string, area float64) ([][]float64, float64) {var size float64if area >= 7000 {size2 := math.Sqrt(area / (7000 * 12321000000))size = math.Round(size2*1000000) / 1000000 // 保留两位小数} else {size = 0.00001}polygon_4490 := fmt.Sprintf(`{"type": "Polygon","coordinates": [[[104,30],[%f,%f],[%f,%f],[%f,%f],[104,30]]]}`, 104.00+size, 30.00, 104.00+size, 30.00+size, 104.00, 30.00+size)polygon_4490area := methods.CalculateAreaByStr(polygon_4490)sql := fmt.Sprintf(`WITH input_geom AS (SELECT ST_SetSRID(ST_GeomFromGeoJSON('%s'),4326) AS geom
),
tiles AS (SELECT grid.tileFROM input_geom AS i,LATERAL ST_SquareGrid(%f, ST_Expand(i.geom, %f)) AS grid(tile)WHERE ST_Intersects(grid.tile, i.geom)
),
centroids AS (SELECT ST_Centroid(tile) AS center,ST_Y(ST_Centroid(tile)) AS lat,ST_X(ST_Centroid(tile)) AS lngFROM tiles
)
SELECT ST_AsGeoJSON(ST_Union(center)) as geojson FROM centroids
`, geo, size, size)type geometryData struct {GeoJSON []byte `gorm:"column:geojson"` // 假设数据库中存储geojson的列名为"geojson"}var geomData geometryDataGeoDB := models.DBerr := GeoDB.Raw(sql).Scan(&geomData)//  如果执行数据库操作时发生错误，向客户端返回错误信息if err.Error != nil {fmt.Println(err.Error)//  提前返回，防止进一步执行}var geo2 geojson.Geometryjson.Unmarshal(geomData.GeoJSON, &geo2)var coords [][]float64switch geom := geo2.Geometry().(type) {case orb.Point:coords = append(coords, []float64{geom.Lon(), geom.Lat()})case orb.MultiPoint:for _, pt := range geom {coords = append(coords, []float64{pt.Lon(), pt.Lat()})}}return coords, polygon_4490area
}

三、填挖方计算

       将前面算法整合，实现代码如下：

func (uc *UserController) GetExcavationFillVolume(c *gin.Context) {var jsonData FillDatac.BindJSON(&jsonData) //将前端geojson转换为geo对象tile := pgmvt.Bounds(jsonData.Geojson.Features[0].Geometry)DB := models.DemDB//获取最大层级var maxZoom int64DB.Model(&models.Tile{}).Select("MAX(zoom_level)").Scan(&maxZoom)//保留最大层级的tilenewtile := pgmvt.KeepTile(maxZoom, tile)//向数据库请求需要的瓦片area := methods.CalculateArea(jsonData.Geojson.Features[0])var conditions []stringvar args []interface{}for _, t := range newtile {conditions = append(conditions, "(tile_column = ? AND tile_row = ? AND zoom_level = ?)")args = append(args, t.X, t.Y, t.Z)}// 执行批量查询var ResultTiles []models.TileDB.Where(strings.Join(conditions, " OR "), args...).Find(&ResultTiles)//获取某个点在//构造三维点字符串pointz1 := jsonData.SJMListvalues := make([]string, len(pointz1))for i, p := range pointz1 {values[i] = fmt.Sprintf(`(%f,%f,%f)`, p[0], p[1], p[2])}pointStr := strings.Join(values, ",")//构造TINgeo := Transformer.GetGeometryString(jsonData.Geojson.Features[0])//Tin := GetHeightFromTIN(geo, pointStr, 104.2969, 30.5503)//tillcoords, P_area := GetTilesPoints(geo, area)Excavation := 0.00Fill := 0.00coords2 := MakeZList(coords, geo, pointStr)var wg sync.WaitGroupresults := make(chan float64, 10)for _, coord := range coords2 {wg.Add(1)go func(x, y, z float64) {defer wg.Done()ZSBG := GetRGBToDEM(x, y, ResultTiles, maxZoom)volume := (z - ZSBG) * P_arearesults <- volume}(coord[0], coord[1], coord[2])}go func() {wg.Wait()close(results)}()for volume := range results {if volume < 0 {Excavation += volume} else {Fill += volume}}type efData struct {Excavation float64Fill       float64}var result efDataresult.Excavation = math.Abs(Excavation)result.Fill = math.Abs(Fill)c.JSON(http.StatusOK, result)}

四、实现效果展示

填挖方

总结

相关文章：

  • Vue3 + Element-Plus + 阿里云文件上传
  • SpringBoot 接口国际化i18n 多语言返回 中英文切换 全球化 语言切换
  • Mioty｜采用报文分割（Telegram Splitting）以提高抗干扰能力的无线通信技术
  • 北极花携语音唤醒、专家鉴定等新功能 亮相第十七届中国林业青年学术年会
  • 继承（c++版 非常详细版）
  • C++ CRC16校验方法详解
  • QT中的多线程
  • Leetcode算法题：字符串转换整数（atoi）
  • ship_plant船舶模型
  • 小草GrassRouter多卡聚合路由器聚合卫星、MESH网络应用解决方案
  • 低功耗蓝牙BLE之高斯频移键控（GFSK）
  • 【Vue.js】组件数据通信：基于Props 实现父组件→子组件传递数据（最基础案例）
  • 前端连接websocket服务报错 Unexpected response code: 301
  • 31、简要描述Promise.all的用途
  • 生成对抗网络（Generative Adversarial Nets，GAN)
  • WebRTC SDK是什么？
  • 机器学习中的数据转换：关键步骤与最佳实践
  • 刀客独家 | 磁力引擎副总裁康乐兼任快手电商美妆跨境行业负责人
  • uniswap getTickAtSqrtPrice 方法解析
  • C语言教程（十八）：C 语言共用体详解
  • 看展览｜建造上海：1949年以来的建筑、城市与文化
  • 君亭酒店：2024年营业收入约6.76亿元， “酒店行业传统增长模式面临巨大挑战”
  • 夜读丨怀念那个写信的年代
  • 他比李白接地气，比杜甫乐观，比白居易刚毅
  • 福建省莆田市原副市长胡国防接受审查调查
  • 四川甘孜州白玉县发生4.9级地震，震源深度10千米