前言
一、面分割
二、面融合
三、捕捉
四、追踪构面
总结
目前市面上大部分WebGIS的产品的编辑功能都很简陋，大部分都只支持简单节点编辑。稍微好一点的会支持面数据裁剪分割、融合。但是在大部分数据使用场景中，为避免出现矢量和矢量之间出现缝隙、压盖等拓扑错误，捕捉和追踪功能就非常重要了。本博客主要讲解如何通过postgis和go语言实现来实现这些功能。
面分割实现，主要是通过绘制的线来切割面数据，需要考虑的是如果线反复穿越面几何，面几何存在多环岛，且线穿过多环岛。这种逻辑如果是通过纯前端turf库是较难实现的。但是如果是通过postgis，就没有那么难实现了，废话不多说直接上代码：
// 图斑分割
func (uc *UserController) SplitGeo(c *gin.Context) {
var jsonData ZDList
c.BindJSON(&jsonData)
var pic []models.GeoPic
DB := models.DB
var Way TempWay
DB.Where("tb_id = ?", jsonData.TBID).Find(&Way.TempGeo)
if len(Way.TempGeo) == 0 {
DB.Where("tb_id = ?", jsonData.TBID).Find(&Way.TempLayer)
}
DB.Where("tb_id = ?", jsonData.TBID).Find(&pic)
var att models.TempLayerAttribute
DB.Where("tb_id = ?", jsonData.TBID).Find(&att)
//线切割面
line := Transformer.GetGeometryString(jsonData.Line.Features[0])
p := GetSingleGeo(jsonData.TBID)
polygon := Transformer.GetGeometryString(p.Features[0])
sql := fmt.Sprintf("WITH geom AS ( SELECT ST_GeomFromGeoJSON('%s') AS line, ST_GeomFromGeoJSON('%s' ) AS polygon) SELECT ST_AsGeoJSON(ST_Split(geom.polygon,geom.line)) AS geojson FROM geom", line, polygon)
var geomData Transformer.GeometryData
err := DB.Raw(sql).Scan(&geomData)
if err.Error != nil {
c.String(http.StatusBadRequest, "err")
return
}
t := p.Features[0].Type
Properties := p.Features[0].Properties
geo := PGBytesToGeojson(geomData)
var NewGeo geojson.FeatureCollection
for index, item := range geo.Geometry {
var feature struct {
Geometry map[string]interface{} `json:"geometry"`
Properties map[string]interface{} `json:"properties"`
Type string `json:"type"`
}
feature.Properties = Properties
feature.Type = t
feature.Geometry = item
data2, _ := json.Marshal(feature)
var myfeature *geojson.Feature
json.Unmarshal(data2, &myfeature)
//判断是图层还是调查
TBID := uuid.New().String()
myfeature.ID = TBID
myfeature.Properties["TBID"] = TBID
//照片ID转移
for pi, _ := range pic {
pic[pi].TBID = TBID
}
DB.Save(pic)
if len(Way.TempGeo) >= 1 {
myfeature.Properties["bsm"] = TBID
myfeature.Properties["name"] = Way.TempGeo[0].Name + strconv.Itoa(index)
newgeo := geojson.NewFeatureCollection()
newgeo.Features = append(newgeo.Features, myfeature)
geoJSONData, err := json.MarshalIndent(newgeo, "", " ")
if err == nil {
result := models.TempGeo{TBID: TBID, Geojson: geoJSONData, MAC: Way.TempGeo[0].MAC, BSM: TBID, Date: time.Now().Format("2006-01-02 15:04:05"), Name: Way.TempGeo[0].Name + strconv.Itoa(index), ZT: Way.TempGeo[0].ZT}
result_att := models.TempLayerAttribute{TBID: TBID, QKSM: att.QKSM, B: att.B, D: att.D, N: att.N, X: att.X, BZ: att.BZ, ZJR: att.ZJR, DCR: att.ZJR}
DB.Create(&result)
DB.Create(&result_att)
}
} else {
Layername := Way.TempLayer[0].Layername
bsm := Way.TempLayer[0].BSM
MAC := Way.TempLayer[0].MAC
myfeature.Properties["name"] = Way.TempLayer[0].Name + strconv.Itoa(index)
newgeo := geojson.NewFeatureCollection()
newgeo.Features = append(newgeo.Features, myfeature)
geoJSONData, err := json.MarshalIndent(newgeo, "", " ")
if err == nil {
result := models.TempLayer{Layername: Layername, MAC: MAC, BSM: bsm, TBID: TBID, Geojson: geoJSONData, ZT: Way.TempLayer[0].ZT, Name: Way.TempLayer[0].Name}
result_att := models.TempLayerAttribute{TBID: TBID, Layername: Layername, QKSM: att.QKSM, B: att.B, D: att.D, N: att.N, X: att.X, BZ: att.BZ, ZJR: att.ZJR, DCR: att.ZJR}
DB.Create(&result_att)
DB.Create(&result)
}
}
NewGeo.Features = append(NewGeo.Features, myfeature)
}
DB.Delete(&att)
//照片处理
//删除原图形
if len(Way.TempGeo) >= 1 {
DB.Delete(&Way.TempGeo)
} else {
DB.Delete(&Way.TempLayer)
}
c.JSON(http.StatusOK, NewGeo)
}分割效果展示：
分割
面融合主要逻辑是需要提前判断两个面是否符合能够融合的条件，比如是否存在共边，是否空间存在重叠，所以需要提前使用ST_Intersects和ST_Touches进行判断。前端代码则需要写好交互逻辑，比如合并需要选择一个要素作为主要属性。
// 图斑合并
type DissolverType struct {
ZD []string `json:"ZD"`
MainTBID string `json:"MainTBID"`
}
func (uc *UserController) DissolverGeo(c *gin.Context) {
var jsonData DissolverType
c.BindJSON(&jsonData)
db := models.DB
var Properties map[string]interface{}
var NewGeo geojson.FeatureCollection
var t string
var geolist []string
var zds []*geojson.FeatureCollection
for _, tt := range jsonData.ZD {
zds = append(zds, GetSingleGeo(tt))
}
var DelID []string
for _, features := range zds {
geo := Transformer.GetGeometryString(features.Features[0])
t = features.Features[0].Type
if features.Features[0].Properties["TBID"] == jsonData.MainTBID {
Properties = features.Features[0].Properties
} else {
DelID = append(DelID, features.Features[0].Properties["TBID"].(string))
}
geolist = append(geolist, geo)
}
var sqlgeo string
var geom []string
var geom1 []string
for index, item := range geolist {
sql1 := fmt.Sprintf("%s AS ( SELECT ST_GeomFromGeoJSON('%s') AS geom),", "geo"+strconv.Itoa(index), item)
sqlgeo = sqlgeo + sql1
geom = append(geom, "geo"+strconv.Itoa(index)+".geom")
geom1 = append(geom1, "geo"+strconv.Itoa(index))
}
result1 := strings.Join(geom, ",")
result2 := strings.Join(geom1, ",")
sql := fmt.Sprintf("WITH %s is_adjacent AS ( SELECT ST_Intersects(%s) AS intersects , ST_Touches(%s) AS touches FROM %s) SELECT CASE WHEN intersects OR touches THEN ST_AsGeoJSON(ST_Union(%s)) ELSE ST_AsGeoJSON(ST_GeomFromGeoJSON('{\"type\": \"GeometryCollection\", \"geometries\": []}')) END AS geojson FROM %s,is_adjacent;", sqlgeo, result1, result1, result2, result1, result2)
var geomData Transformer.GeometryData
err := db.Raw(sql).Scan(&geomData)
if err.Error != nil {
c.String(http.StatusBadRequest, "err")
return
}
var feature struct {
Geometry map[string]interface{} `json:"geometry"`
Properties map[string]interface{} `json:"properties"`
Type string `json:"type"`
}
feature.Properties = Properties
feature.Type = t
json.Unmarshal(geomData.GeoJSON, &feature.Geometry)
if feature.Geometry["geometries"] != nil {
if len(feature.Geometry["geometries"].([]interface{})) == 0 {
c.String(http.StatusBadRequest, "err")
return
}
}
data2, _ := json.Marshal(feature)
var myfeature *geojson.Feature
json.Unmarshal(data2, &myfeature)
myfeature.ID = jsonData.MainTBID
NewGeo.Features = append(NewGeo.Features, myfeature)
geoJSONData, _ := json.MarshalIndent(NewGeo, "", " ")
var templayer []models.TempLayer
db.Where("tb_id = ?", jsonData.MainTBID).Find(&templayer)
if len(templayer) >= 1 {
templayer[0].Geojson = geoJSONData
db.Save(&templayer)
} else {
var templayer2 models.TempGeo
db.Where("tb_id = ?", jsonData.MainTBID).Find(&templayer2)
templayer2.Geojson = geoJSONData
db.Save(&templayer2)
}
db.Where("tb_id IN (?)", DelID).Delete(&models.TempLayer{})
db.Where("tb_id IN (?)", DelID).Delete(&models.TempGeo{})
db.Where("tb_id IN (?)", DelID).Delete(&models.TempLayerAttribute{})
c.JSON(http.StatusOK, NewGeo)
}融合效果：
合并
捕捉的实现需要实现几个逻辑，第一个就是判断当前需要捕捉的图层，第二个是判断需要捕捉的对象的几何类型（点、线、面），第三个是需要判断捕捉的对象是矢量瓦片，还是geojson。代码如下，这里进行了6种情况的的判断并返回捕捉的点：
// 图形捕捉
func (uc *UserController) Capture(c *gin.Context) {
var jsonData CaptureType
c.BindJSON(&jsonData)
DB := models.DB
Templarers := jsonData.TempLayer
//临时图层同时存在
if len(Templarers) != 0 && jsonData.Layer != "" {
var sql string
geojsons := GetTempLayers(Templarers)
//点线面分离
var PolygonJson []*geojson.Feature
var LineJson []*geojson.Feature
var PointJson []*geojson.Feature
for _, item := range geojsons {
switch item.Geometry.GeoJSONType() {
case "Polygon":
PolygonJson = append(PolygonJson, item)
case "LineString":
LineJson = append(LineJson, item)
case "Point":
PointJson = append(PointJson, item)
}
}
//面数据捕捉
if len(PolygonJson) != 0 {
geo := Transformer.GetFeatureString(PolygonJson)
sql = fmt.Sprintf(`
WITH input_point AS (
SELECT ST_SetSRID(ST_MakePoint(%f, %f), 4326) AS geom
),
geojson_data AS (
SELECT
jsonb_array_elements('%s'::jsonb) AS feature
),
nearest_polygon AS (
SELECT
ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326) AS poly_geom
FROM
geojson_data, input_point AS input
WHERE
ST_Distance(input.geom, ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326)) < 0.0005
ORDER BY
ST_Distance(input.geom, ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326))
LIMIT 1
),
line_geom AS (
SELECT
ST_Boundary(poly_geom) AS line_geom
FROM
nearest_polygon
)
SELECT
ST_AsGeoJSON(ST_ClosestPoint(line_geom.line_geom, input.geom)) AS geojson,
ST_Distance(input.geom, ST_ClosestPoint(line_geom.line_geom, input.geom))::float AS distance
FROM
input_point AS input, line_geom;
`, jsonData.Point[0], jsonData.Point[1], geo)
}
//线数据捕捉
if len(LineJson) != 0 {
geo := Transformer.GetFeatureString(LineJson)
sql = fmt.Sprintf(`
WITH input_point AS (
SELECT ST_SetSRID(ST_MakePoint(%f, %f), 4326) AS geom
),
geojson_data AS (
SELECT
jsonb_array_elements('%s'::jsonb) AS feature
),
nearest_line AS (
SELECT
ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326) AS line_geom
FROM
geojson_data, input_point AS input
WHERE
ST_Distance(input.geom, ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326)) < 0.0005
ORDER BY
ST_Distance(input.geom, ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326))
LIMIT 1
)
SELECT
ST_AsGeoJSON(ST_ClosestPoint(line_geom, input.geom)) AS geojson,
ST_Distance(input.geom, ST_ClosestPoint(line_geom, input.geom))::float AS distance
FROM
input_point AS input, nearest_line;
`, jsonData.Point[0], jsonData.Point[1], geo)
}
//点数据捕捉
if len(PointJson) != 0 {
geo := Transformer.GetFeatureString(PointJson)
sql = fmt.Sprintf(`
WITH input_point AS (
SELECT ST_SetSRID(ST_MakePoint(%f, %f), 4326) AS geom
),
geojson_data AS (
SELECT
jsonb_array_elements('%s'::jsonb) AS feature
),
nearest_point AS (
SELECT
ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326) AS point_geom
FROM
geojson_data, input_point AS input
WHERE
ST_Distance(input.geom, ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326)) < 0.0005
ORDER BY
ST_Distance(input.geom, ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326))
LIMIT 1
)
SELECT
ST_AsGeoJSON(point_geom) AS geojson,
ST_Distance(input.geom, point_geom)::float AS distance
FROM
input_point AS input, nearest_point;
`, jsonData.Point[0], jsonData.Point[1], geo)
}
var geomData CaptureData
err := DB.Raw(sql).Scan(&geomData)
if err.Error != nil {
c.String(http.StatusBadRequest, "err")
return
}
var feature struct {
Geometry GeometryPoint `json:"geometry"`
}
json.Unmarshal(geomData.GeoJSON, &feature.Geometry)
if geomData.Distance <= 0.00015 && len(feature.Geometry.Coordinates) != 0 {
c.JSON(http.StatusOK, feature.Geometry.Coordinates)
} else {
layer := jsonData.Layer
var schema []models.MySchema
DB.Where("en = ?", layer).Find(&schema)
if len(schema) >= 1 {
switch schema[0].Type {
case "polygon":
sql = fmt.Sprintf(`
WITH input_point AS (
SELECT ST_SetSRID(ST_MakePoint(%f, %f), 4326) AS geom
),
nearest_polygon AS (
SELECT
%s.geom AS poly_geom
FROM
%s, input_point AS input
WHERE
ST_Distance(input.geom, %s.geom) < 0.0005
ORDER BY
ST_Distance(input.geom, %s.geom)
LIMIT 1
),
line_geom AS (
SELECT
ST_Boundary(poly_geom) AS line_geom
FROM
nearest_polygon
)
SELECT
ST_AsGeoJSON(ST_ClosestPoint(line_geom.line_geom, input.geom)) AS geojson,
ST_Distance(input.geom, ST_ClosestPoint(line_geom.line_geom, input.geom))::float AS distance
FROM
input_point AS input, line_geom;
`, jsonData.Point[0], jsonData.Point[1], jsonData.Layer, jsonData.Layer, jsonData.Layer, jsonData.Layer)
case "line":
sql = fmt.Sprintf(`
WITH input_point AS (
SELECT ST_SetSRID(ST_MakePoint(%f, %f), 4326) AS geom
),
closest_line AS (
SELECT
%s.geom AS line_geom
FROM
%s, input_point AS input
WHERE
ST_Distance(input.geom, %s.geom) < 0.0005
ORDER BY
ST_Distance(input.geom, %s.geom)
LIMIT 1
)
SELECT
ST_AsGeoJSON(ST_ClosestPoint(closest_line.line_geom, input.geom)) AS geojson,
ST_Distance(input.geom, ST_ClosestPoint(closest_line.line_geom, input.geom))::float AS distance
FROM
input_point AS input, closest_line;
`, jsonData.Point[0], jsonData.Point[1], jsonData.Layer, jsonData.Layer, jsonData.Layer, jsonData.Layer)
case "point":
sql = fmt.Sprintf(`
WITH input_point AS (
SELECT ST_SetSRID(ST_MakePoint(%f, %f), 4326) AS geom
),
nearest_point AS (
SELECT
%s.geom AS nearest_geom
FROM
%s, input_point AS input
ORDER BY
ST_Distance(input.geom, %s.geom)
LIMIT 1
)
SELECT
ST_AsGeoJSON(nearest_geom) AS geojson,
ST_Distance(input.geom, nearest_geom)::float AS distance
FROM
input_point AS input, nearest_point;
`, jsonData.Point[0], jsonData.Point[1], jsonData.Layer, jsonData.Layer, jsonData.Layer)
}
}
err := DB.Raw(sql).Scan(&geomData)
if err.Error != nil {
c.String(http.StatusBadRequest, "err")
return
}
json.Unmarshal(geomData.GeoJSON, &feature.Geometry)
if geomData.Distance <= 0.00015 {
c.JSON(http.StatusOK, feature.Geometry.Coordinates)
} else {
c.String(http.StatusBadRequest, "err")
}
}
} else {
var sql string
if len(Templarers) != 0 && jsonData.Layer == "" {
geojsons := GetTempLayers(Templarers)
//点线面分离
var PolygonJson []*geojson.Feature
var LineJson []*geojson.Feature
var PointJson []*geojson.Feature
for _, item := range geojsons {
switch item.Geometry.GeoJSONType() {
case "Polygon":
PolygonJson = append(PolygonJson, item)
case "LineString":
LineJson = append(LineJson, item)
case "Point":
PointJson = append(PointJson, item)
}
}
//面数据捕捉
if len(PolygonJson) != 0 {
geo := Transformer.GetFeatureString(PolygonJson)
sql = fmt.Sprintf(`
WITH input_point AS (
SELECT ST_SetSRID(ST_MakePoint(%f, %f), 4326) AS geom
),
geojson_data AS (
SELECT
jsonb_array_elements('%s'::jsonb) AS feature
),
nearest_polygon AS (
SELECT
ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326) AS poly_geom
FROM
geojson_data, input_point AS input
WHERE
ST_Distance(input.geom, ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326)) < 0.0005
ORDER BY
ST_Distance(input.geom, ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326))
LIMIT 1
),
line_geom AS (
SELECT
ST_Boundary(poly_geom) AS line_geom
FROM
nearest_polygon
)
SELECT
ST_AsGeoJSON(ST_ClosestPoint(line_geom.line_geom, input.geom)) AS geojson,
ST_Distance(input.geom, ST_ClosestPoint(line_geom.line_geom, input.geom))::float AS distance
FROM
input_point AS input, line_geom;
`, jsonData.Point[0], jsonData.Point[1], geo)
}
//线数据捕捉
if len(LineJson) != 0 {
geo := Transformer.GetFeatureString(LineJson)
sql = fmt.Sprintf(`
WITH input_point AS (
SELECT ST_SetSRID(ST_MakePoint(%f, %f), 4326) AS geom
),
geojson_data AS (
SELECT
jsonb_array_elements('%s'::jsonb) AS feature
),
nearest_line AS (
SELECT
ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326) AS line_geom
FROM
geojson_data, input_point AS input
WHERE
ST_Distance(input.geom, ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326)) < 0.0005
ORDER BY
ST_Distance(input.geom, ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326))
LIMIT 1
)
SELECT
ST_AsGeoJSON(ST_ClosestPoint(line_geom, input.geom)) AS geojson,
ST_Distance(input.geom, ST_ClosestPoint(line_geom, input.geom))::float AS distance
FROM
input_point AS input, nearest_line;
`, jsonData.Point[0], jsonData.Point[1], geo)
}
//点数据捕捉
if len(PointJson) != 0 {
geo := Transformer.GetFeatureString(PointJson)
sql = fmt.Sprintf(`
WITH input_point AS (
SELECT ST_SetSRID(ST_MakePoint(%f, %f), 4326) AS geom
),
geojson_data AS (
SELECT
jsonb_array_elements('%s'::jsonb) AS feature
),
nearest_point AS (
SELECT
ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326) AS point_geom
FROM
geojson_data, input_point AS input
WHERE
ST_Distance(input.geom, ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326)) < 0.0005
ORDER BY
ST_Distance(input.geom, ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326))
LIMIT 1
)
SELECT
ST_AsGeoJSON(point_geom) AS geojson,
ST_Distance(input.geom, point_geom)::float AS distance
FROM
input_point AS input, nearest_point;
`, jsonData.Point[0], jsonData.Point[1], geo)
}
} else if len(Templarers) == 0 && jsonData.Layer != "" {
layer := jsonData.Layer
var schema []models.MySchema
DB.Where("en = ?", layer).Find(&schema)
if len(schema) >= 1 {
switch schema[0].Type {
case "polygon":
sql = fmt.Sprintf(`
WITH input_point AS (
SELECT ST_SetSRID(ST_MakePoint(%f, %f), 4326) AS geom
),
nearest_polygon AS (
SELECT
%s.geom AS poly_geom
FROM
%s, input_point AS input
WHERE
ST_Distance(input.geom, %s.geom) < 0.0005
ORDER BY
ST_Distance(input.geom, %s.geom)
LIMIT 1
),
line_geom AS (
SELECT
ST_Boundary(poly_geom) AS line_geom
FROM
nearest_polygon
)
SELECT
ST_AsGeoJSON(ST_ClosestPoint(line_geom.line_geom, input.geom)) AS geojson,
ST_Distance(input.geom, ST_ClosestPoint(line_geom.line_geom, input.geom))::float AS distance
FROM
input_point AS input, line_geom;
`, jsonData.Point[0], jsonData.Point[1], jsonData.Layer, jsonData.Layer, jsonData.Layer, jsonData.Layer)
case "line":
sql = fmt.Sprintf(`
WITH input_point AS (
SELECT ST_SetSRID(ST_MakePoint(%f, %f), 4326) AS geom
),
closest_line AS (
SELECT
%s.geom AS line_geom
FROM
%s, input_point AS input
WHERE
ST_Distance(input.geom, %s.geom) < 0.0005
ORDER BY
ST_Distance(input.geom, %s.geom)
LIMIT 1
)
SELECT
ST_AsGeoJSON(ST_ClosestPoint(closest_line.line_geom, input.geom)) AS geojson,
ST_Distance(input.geom, ST_ClosestPoint(closest_line.line_geom, input.geom))::float AS distance
FROM
input_point AS input, closest_line;
`, jsonData.Point[0], jsonData.Point[1], jsonData.Layer, jsonData.Layer, jsonData.Layer, jsonData.Layer)
case "point":
sql = fmt.Sprintf(`
WITH input_point AS (
SELECT ST_SetSRID(ST_MakePoint(%f, %f), 4326) AS geom
),
nearest_point AS (
SELECT
%s.geom AS nearest_geom
FROM
%s, input_point AS input
ORDER BY
ST_Distance(input.geom, %s.geom)
LIMIT 1
)
SELECT
ST_AsGeoJSON(nearest_geom) AS geojson,
ST_Distance(input.geom, nearest_geom)::float AS distance
FROM
input_point AS input, nearest_point;
`, jsonData.Point[0], jsonData.Point[1], jsonData.Layer, jsonData.Layer, jsonData.Layer)
}
}
}
var geomData CaptureData
err := DB.Raw(sql).Scan(&geomData)
if err.Error != nil {
c.String(http.StatusBadRequest, "err")
return
}
var feature struct {
Geometry GeometryPoint `json:"geometry"`
}
json.Unmarshal(geomData.GeoJSON, &feature.Geometry)
if geomData.Distance <= 0.00015 {
c.JSON(http.StatusOK, feature.Geometry.Coordinates)
} else {
c.String(http.StatusBadRequest, "err")
}
}
}捕捉效果：
捕捉
追踪构面的实现逻辑可以说是最为复杂的，实现逻辑主要是根据绘制的线段和当前已经加载好了的图层进行空间分析，查询到最为合适的边界，并将该边界与绘制线进行构面。分解下来的逻辑大概是：
1、因为容差因素，需要将传入的线按照节点顺序进行两端延长，保证线会超过图层线并形成交点。
2、将延长线与几何图层进行空间分析，将与线有交集的几何提取出来
3、如果几何是面，需要转换为线要素
4、将输入线与相交线进行节点打散，并重新构造面
5、再将输入线与重新构造好的面进行线线叠加分析，找出重叠率最高的几何，并返回该几何
代码如下：
type AutoData struct {
Line geojson.FeatureCollection `json:"Line"`
Layer string
TempLayer []string
}
func (uc *UserController) AutoPolygon(c *gin.Context) {
var jsonData AutoData
c.BindJSON(&jsonData)
DB := models.DB
var sql string
line := Transformer.GetGeometryString(jsonData.Line.Features[0])
Templarers := jsonData.TempLayer
if len(Templarers) != 0 {
geojsons := GetTempLayers(Templarers)
//点线面分离
var PolygonJson []*geojson.Feature
var LineJson []*geojson.Feature
for _, item := range geojsons {
switch item.Geometry.GeoJSONType() {
case "Polygon":
PolygonJson = append(PolygonJson, item)
case "LineString":
LineJson = append(LineJson, item)
}
}
//面数据捕捉
if len(PolygonJson) != 0 {
geo := Transformer.GetFeatureString(PolygonJson)
sql = fmt.Sprintf(`
WITH input_linefrist AS (
SELECT ST_SetSRID(ST_GeomFromGeoJSON('%s'), 4326) AS geom
),
input_line AS (
SELECT
ST_LineMerge(ST_Union( ARRAY[
ST_MakeLine(ST_StartPoint(geom),ST_Project(ST_StartPoint(geom), 0.000001, ST_Azimuth(ST_PointN(geom,2),ST_StartPoint(geom)))::geometry),
geom,
ST_MakeLine(ST_EndPoint(geom),ST_Project(ST_EndPoint(geom), 0.000001, ST_Azimuth(ST_PointN(geom,ST_NumPoints(geom) - 1),ST_EndPoint(geom)))::geometry)
])) AS geom
FROM input_linefrist
),
intersecting_areas AS (
SELECT ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326) AS geom
FROM jsonb_array_elements('%s'::jsonb) AS feature
WHERE ST_Intersects(ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326), (SELECT geom FROM input_line))
),
boundary_lines AS (
SELECT ST_Boundary(geom) AS geom
FROM intersecting_areas
),
closed_geometries AS (
SELECT ST_Union(ST_Collect((SELECT geom FROM input_line), boundary_lines.geom)) AS lines
FROM boundary_lines
),
newpolygons AS (
SELECT ST_Polygonize(lines) AS polygon_geoms
FROM closed_geometries
),
boundary_lines2 AS (
SELECT (ST_Dump(ST_Boundary(polygon_geoms))).geom AS geom
FROM newpolygons
),
intersecting_lines AS (
SELECT ST_Intersection(input_line.geom, boundary_lines2.geom,0.0000001) AS geom, boundary_lines2.geom AS boundary_geom
FROM input_line, boundary_lines2
),
max_overlap AS (
SELECT
boundary_geom,
ST_Length(geom)/ST_Length(boundary_geom) AS overlap_length
FROM intersecting_lines
ORDER BY overlap_length DESC
LIMIT 1
)
SELECT ST_AsGeoJSON(ST_MakePolygon(boundary_geom)) AS geojson,
overlap_length::float AS lenth
FROM max_overlap
`, line, geo)
}
if len(LineJson) != 0 {
geo := Transformer.GetFeatureString(LineJson)
sql = fmt.Sprintf(`
WITH input_linefrist AS (
SELECT ST_SetSRID(ST_GeomFromGeoJSON('%s'), 4326) AS geom
),
input_line AS (
SELECT
ST_LineMerge(ST_Union( ARRAY[
ST_MakeLine(ST_StartPoint(geom),ST_Project(ST_StartPoint(geom), 0.000001, ST_Azimuth(ST_PointN(geom,2),ST_StartPoint(geom)))::geometry),
geom,
ST_MakeLine(ST_EndPoint(geom),ST_Project(ST_EndPoint(geom), 0.000001, ST_Azimuth(ST_PointN(geom,ST_NumPoints(geom) - 1),ST_EndPoint(geom)))::geometry)
])) AS geom
FROM input_linefrist
),
intersecting_areas AS (
SELECT ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326) AS geom
FROM jsonb_array_elements('%s'::jsonb) AS feature
WHERE ST_Intersects(ST_SetSRID(ST_GeomFromGeoJSON(feature->'geometry'::text), 4326), (SELECT geom FROM input_line))
),
boundary_lines AS (
SELECT geom AS geom
FROM intersecting_areas
),
closed_geometries AS (
SELECT ST_Union(ST_Collect((SELECT geom FROM input_line), boundary_lines.geom)) AS lines
FROM boundary_lines
),
newpolygons AS (
SELECT ST_Polygonize(lines) AS polygon_geoms
FROM closed_geometries
),
boundary_lines2 AS (
SELECT (ST_Dump(ST_Boundary(polygon_geoms))).geom AS geom
FROM newpolygons
),
intersecting_lines AS (
SELECT ST_Intersection(input_line.geom, boundary_lines2.geom,0.0000001) AS geom, boundary_lines2.geom AS boundary_geom
FROM input_line, boundary_lines2
),
max_overlap AS (
SELECT
boundary_geom,
ST_Length(geom)/ST_Length(boundary_geom) AS overlap_length
FROM intersecting_lines
ORDER BY overlap_length DESC
LIMIT 1
)
SELECT ST_AsGeoJSON(ST_MakePolygon(boundary_geom)) AS geojson,
overlap_length::float AS lenth
FROM max_overlap
`, line, geo)
}
} else {
layer := jsonData.Layer
var schema []models.MySchema
DB.Where("en = ?", layer).Find(&schema)
if len(schema) >= 1 {
switch schema[0].Type {
case "polygon":
sql = fmt.Sprintf(`
WITH input_linefrist AS (
SELECT ST_SetSRID(ST_GeomFromGeoJSON('%s'), 4326) AS geom
),
input_line AS (
SELECT
ST_LineMerge(ST_Union( ARRAY[
ST_MakeLine(ST_StartPoint(geom),ST_Project(ST_StartPoint(geom), 0.000001, ST_Azimuth(ST_PointN(geom,2),ST_StartPoint(geom)))::geometry),
geom,
ST_MakeLine(ST_EndPoint(geom),ST_Project(ST_EndPoint(geom), 0.000001, ST_Azimuth(ST_PointN(geom,ST_NumPoints(geom) - 1),ST_EndPoint(geom)))::geometry)
])) AS geom
FROM input_linefrist
),
intersecting_areas AS (
SELECT *
FROM %s
WHERE ST_Intersects(%s.geom, (SELECT geom FROM input_line))
),
boundary_lines AS (
SELECT ST_Boundary(geom) AS geom
FROM intersecting_areas
),
closed_geometries AS (
SELECT ST_Union(ST_Collect((SELECT geom FROM input_line), boundary_lines.geom)) AS lines
FROM boundary_lines
),
newpolygons AS (
SELECT ST_Polygonize(lines) AS polygon_geoms
FROM closed_geometries
),
boundary_lines2 AS (
SELECT (ST_Dump(ST_Boundary(polygon_geoms))).geom AS geom
FROM newpolygons
),
intersecting_lines AS (
SELECT ST_Intersection(input_line.geom, boundary_lines2.geom,0.0000001) AS geom, boundary_lines2.geom AS boundary_geom
FROM input_line, boundary_lines2
),
max_overlap AS (
SELECT
boundary_geom,
ST_Length(geom)/ST_Length(boundary_geom) AS overlap_length
FROM intersecting_lines
ORDER BY overlap_length DESC
LIMIT 1
)
SELECT ST_AsGeoJSON(ST_MakePolygon(boundary_geom)) AS geojson,
overlap_length::float AS lenth
FROM max_overlap
`, line, jsonData.Layer, jsonData.Layer)
case "line":
sql = fmt.Sprintf(`
WITH input_linefrist AS (
SELECT ST_SetSRID(ST_GeomFromGeoJSON('%s'), 4326) AS geom
),
input_line AS (
SELECT
ST_LineMerge(ST_Union( ARRAY[
ST_MakeLine(ST_StartPoint(geom),ST_Project(ST_StartPoint(geom), 0.000001, ST_Azimuth(ST_PointN(geom,2),ST_StartPoint(geom)))::geometry),
geom,
ST_MakeLine(ST_EndPoint(geom),ST_Project(ST_EndPoint(geom), 0.000001, ST_Azimuth(ST_PointN(geom,ST_NumPoints(geom) - 1),ST_EndPoint(geom)))::geometry)
])) AS geom
FROM input_linefrist
),
intersecting_areas AS (
SELECT *
FROM %s
WHERE ST_Intersects(%s.geom, (SELECT geom FROM input_line))
),
boundary_lines AS (
SELECT geom AS geom
FROM intersecting_areas
),
closed_geometries AS (
SELECT ST_Union(ST_Collect((SELECT geom FROM input_line), boundary_lines.geom)) AS lines
FROM boundary_lines
),
newpolygons AS (
SELECT ST_Polygonize(lines) AS polygon_geoms
FROM closed_geometries
),
boundary_lines2 AS (
SELECT (ST_Dump(ST_Boundary(polygon_geoms))).geom AS geom
FROM newpolygons
),
intersecting_lines AS (
SELECT ST_Intersection(input_line.geom, boundary_lines2.geom,0.0000001) AS geom, boundary_lines2.geom AS boundary_geom
FROM input_line, boundary_lines2
WHERE ST_Intersects(input_line.geom, boundary_lines2.geom)
),
max_overlap AS (
SELECT
boundary_geom,
ST_Length(geom)/ST_Length(boundary_geom) AS overlap_length
FROM intersecting_lines
ORDER BY overlap_length DESC
LIMIT 1
)
SELECT ST_AsGeoJSON(ST_MakePolygon(boundary_geom)) AS geojson,
overlap_length::float AS lenth
FROM max_overlap
`, line, jsonData.Layer, jsonData.Layer)
}
}
}
var geomData Transformer.GeometryData
err := DB.Raw(sql).Scan(&geomData)
if err.Error != nil {
c.String(http.StatusBadRequest, "err")
return
}
var feature struct {
Geometry map[string]interface{} `json:"geometry"`
Properties map[string]interface{} `json:"properties"`
Type string `json:"type"`
}
feature.Type = "Feature"
json.Unmarshal(geomData.GeoJSON, &feature.Geometry)
feature.Properties = make(map[string]interface{})
feature.Properties["name"] = ""
var NewGeo geojson.FeatureCollection
data2, _ := json.Marshal(feature)
var myfeature *geojson.Feature
aa := json.Unmarshal(data2, &myfeature)
if aa != nil {
fmt.Println(aa.Error())
}
NewGeo.Features = append(NewGeo.Features, myfeature)
c.JSON(http.StatusOK, NewGeo)
}
追踪构面效果如下：
追踪
这是我团队的移动端离线GIS产品《新源地图》的绘制模块的开发日志，该软件采用前后端部署一体式方案，将go语言开发的后端和postgis在安卓设备上进行编译。实现了在移动设备上能进行GIS分析、绘制功能。并通过动态瓦片技术实现了离线毫秒级加载千万级的矢量数据。如果对该产品有兴趣的同学可以在后台联系我。