WebGIS开发之编辑功能（分割、融合、捕捉、追踪）

前言

目前市面上大部分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分析、绘制功能。并通过动态瓦片技术实现了离线毫秒级加载千万级的矢量数据。如果对该产品有兴趣的同学可以在后台联系我。