空间计算操作

概述

空间计算操作是对几何对象进行各种空间运算，生成新的几何对象。本章介绍常用的空间计算操作，包括缓冲区分析、叠加分析、几何简化等。

缓冲区分析

基本缓冲区

缓冲区是围绕几何对象指定距离的区域：

/**
 * 获取几何缓冲区
 * @param geom 几何对象
 * @param distance 缓冲区距离（坐标系单位）
 */
public static Geometry buffer(Geometry geom, double distance) {
    return geom.buffer(distance);
}

使用示例：

// 点缓冲区：生成圆形
Point point = factory.createPoint(new Coordinate(100, 100));
Geometry circle = buffer(point, 50);

// 线缓冲区：生成条带状区域
LineString line = factory.createLineString(coords);
Geometry corridor = buffer(line, 10);

// 面缓冲区：生成扩大的面
Polygon polygon = factory.createPolygon(shell);
Geometry expanded = buffer(polygon, 20);

负值缓冲区：

// 负值缓冲区会向内收缩
Geometry shrunk = buffer(polygon, -5);

ESRI几何引擎缓冲区

/**
 * 使用ESRI引擎计算缓冲区
 */
public static String buffer(String wkt, Integer wkid, double distance) {
    Geometry geom = EsriUtil.createGeometryByWkt(wkt);
    SpatialReference sr = SpatialReference.create(wkid);
    Geometry buffer = OperatorBuffer.local().execute(geom, sr, distance, null);
    return EsriUtil.getWktStr(buffer);
}

缓冲区注意事项

坐标系影响：

• 地理坐标系：距离单位为度，结果在两极和赤道处变形不同
• 投影坐标系：距离单位为米，结果更准确

推荐做法：

// 先转换到投影坐标系
int projWkid = CrsUtil.getProjectedWkid(geom);
Geometry projGeom = CrsUtil.transform(geom, 4490, projWkid);

// 在投影坐标系下做缓冲区
Geometry buffer = projGeom.buffer(distance);

// 如需要，转回地理坐标系
Geometry geoBuffer = CrsUtil.transform(buffer, projWkid, 4490);

叠加分析

交集（Intersection）

获取两个几何对象的公共部分：

/**
 * 获取几何交集
 */
public static Geometry intersection(Geometry a, Geometry b) {
    return a.intersection(b);
}

应用场景：

• 裁剪：按边界裁剪数据
• 统计：计算重叠面积
• 提取：提取两个图层的公共区域

并集（Union）

合并多个几何对象：

/**
 * 获取几何并集
 */
public static Geometry union(Geometry... geoms) {
    Geometry result = null;
    for (Geometry g : geoms) {
        if (result == null) {
            result = g;
        } else {
            result = result.union(g);
        }
    }
    return result;
}

应用场景：

• 合并相邻地块
• 溶解（Dissolve）操作
• 生成复合区域

差集（Difference）

从一个几何中减去另一个几何：

/**
 * 获取A与B并集擦除B的部分
 * 结果为A中不与B重叠的部分
 */
public static Geometry difference(Geometry a, Geometry b) {
    return a.difference(b);
}

应用场景：

• 擦除操作：从区域中移除指定部分
• 挖洞：在多边形中创建岛洞
• 更新：用新数据替换旧数据的部分区域

对称差集（Symmetric Difference）

获取两个几何不重叠的部分：

/**
 * 获取A与B并集减去A与B交集
 * 结果为A和B各自独有的部分
 */
public static Geometry symDifference(Geometry a, Geometry b) {
    return a.symDifference(b);
}

应用场景：

• 找出两个版本数据的差异区域
• 识别变化区域

ESRI引擎叠加操作

/**
 * 计算几何的交集
 */
public static String intersection(String awkt, String bwkt, Integer wkid) {
    Geometry a = EsriUtil.createGeometryByWkt(awkt);
    Geometry b = EsriUtil.createGeometryByWkt(bwkt);
    SpatialReference sr = SpatialReference.create(wkid);
    Geometry intersection = OperatorIntersection.local().execute(a, b, sr, null);
    return EsriUtil.getWktStr(intersection);
}

/**
 * 计算几何的并集
 */
public static String union(List<String> wkts, Integer wkid) {
    Geometry[] geoms = wkts.stream()
        .map(EsriUtil::createGeometryByWkt)
        .toArray(Geometry[]::new);
    SpatialReference sr = SpatialReference.create(wkid);
    Geometry union = GeometryEngine.union(geoms, sr);
    return EsriUtil.getWktStr(union);
}

/**
 * 获取A与B并集擦除B的部分
 */
public static String difference(String awkt, String bwkt, Integer wkid) {
    Geometry a = EsriUtil.createGeometryByWkt(awkt);
    Geometry b = EsriUtil.createGeometryByWkt(bwkt);
    SpatialReference sr = SpatialReference.create(wkid);
    Geometry difference = OperatorDifference.local().execute(a, b, sr, null);
    return EsriUtil.getWktStr(difference);
}

凸包与凹包

凸包（Convex Hull）

凸包是包围几何的最小凸多边形：

/**
 * 获取几何凸包
 */
public static Geometry convexHull(Geometry geom) {
    ConvexHull ch = new ConvexHull(geom);
    return ch.getConvexHull();
}

特点：凸包没有任何凹陷，像是用橡皮筋围住所有点。

凹包（Concave Hull）

凹包是更紧密地包围几何的边界：

/**
 * 获取几何凹包
 */
public static Geometry concaveHull(Geometry geom) {
    ConcaveHull ch = new ConcaveHull(geom);
    return ch.getHull();
}

特点：凹包会沿着点集的边缘形成凹陷。

ESRI引擎凸包：

/**
 * 计算几何的凸包
 */
public static String convexHull(String wkt) {
    Geometry geom = EsriUtil.createGeometryByWkt(wkt);
    Geometry convexHull = OperatorConvexHull.local().execute(geom, null);
    return EsriUtil.getWktStr(convexHull);
}

几何简化

Douglas-Peucker简化

减少几何的节点数量，同时保持形状特征：

/**
 * 简化几何
 * @param distance 简化容差，越大简化程度越高
 */
public static Geometry simplify(Geometry geom, double distance) {
    return DouglasPeuckerSimplifier.simplify(geom, distance);
}

应用场景：

• 数据压缩：减少数据量
• 地图显示：不同比例尺显示不同详细程度
• 性能优化：简化复杂几何提升计算速度

ESRI引擎简化

/**
 * 简化几何
 */
public static String simplify(String wkt, Integer wkid) {
    Geometry geom = EsriUtil.createGeometryByWkt(wkt);
    SpatialReference sr = SpatialReference.create(wkid);
    Geometry simplified = OperatorSimplifyOGC.local().execute(geom, sr, false, null);
    return EsriUtil.getWktStr(simplified);
}

几何增密

增加节点密度

在几何的边上添加额外的节点：

/**
 * 增加几何节点密度
 * @param distance 节点间最大距离
 */
public static Geometry densify(Geometry geom, double distance) {
    return Densifier.densify(geom, distance);
}

应用场景：

• 坐标转换前的预处理：防止长边在投影后变形
• 曲线拟合：为曲线化操作做准备

多边形化

从线生成面

将线几何转换为多边形：

/**
 * 多边形化几何
 */
public static Geometry polygonize(Geometry geom) {
    List lines = LineStringExtracter.getLines(geom);
    Polygonizer polygonizer = new Polygonizer();
    polygonizer.add(lines);
    Collection polys = polygonizer.getPolygons();
    Polygon[] polyArray = GeometryFactory.toPolygonArray(polys);
    return geom.getFactory().createGeometryCollection(polyArray);
}

应用场景：

• 从道路网络生成地块
• 从边界线生成区域

切割操作

用线切割多边形

/**
 * 按照给定线切割多边形
 */
public static Geometry splitPolygon(Geometry polygon, LineString line) {
    // 将切割线与多边形边界合并
    Geometry nodedLinework = polygon.getBoundary().union(line);
    
    // 多边形化
    Geometry polys = polygonize(nodedLinework);
    
    // 筛选在原多边形内部的部分
    List<Polygon> output = new ArrayList<>();
    for (int i = 0; i < polys.getNumGeometries(); i++) {
        Polygon candpoly = (Polygon) polys.getGeometryN(i);
        if (polygon.contains(candpoly.getInteriorPoint())) {
            output.add(candpoly);
        }
    }
    
    return polygon.getFactory()
        .createGeometryCollection(GeometryFactory.toGeometryArray(output));
}

应用场景：

• 地块分割
• 区域划分

实践案例

案例1：缓冲区分析

计算道路两侧一定范围内的影响区域：

/**
 * 道路缓冲区分析
 */
public WktLayer roadBufferAnalysis(WktLayer roadLayer, double bufferDistance) {
    WktLayer resultLayer = new WktLayer();
    resultLayer.setYwName("road_buffer");
    resultLayer.setWkid(roadLayer.getWkid());
    resultLayer.setGeometryType(GeometryType.POLYGON);
    resultLayer.setFields(roadLayer.getFields());
    
    List<WktFeature> features = new ArrayList<>();
    
    // 转换到投影坐标系
    int projWkid = CrsUtil.getProjectedWkid(38);  // 假设38带
    WktLayer projLayer = CrsUtil.reproject(roadLayer, projWkid);
    
    for (WktFeature feature : projLayer.getFeatures()) {
        Geometry geom = GeometryConverter.wkt2Geometry(feature.getWkt());
        Geometry buffer = geom.buffer(bufferDistance);
        
        WktFeature bufferFeature = new WktFeature();
        bufferFeature.setWfId(feature.getWfId());
        bufferFeature.setWkt(buffer.toText());
        bufferFeature.setFieldValues(feature.getFieldValues());
        features.add(bufferFeature);
    }
    
    resultLayer.setFeatures(features);
    
    // 转回原坐标系
    return CrsUtil.reproject(resultLayer, roadLayer.getWkid());
}

案例2：图层叠加分析

计算两个图层的重叠区域：

/**
 * 图层叠加分析
 */
public WktLayer overlayAnalysis(WktLayer layer1, WktLayer layer2) {
    WktLayer resultLayer = new WktLayer();
    resultLayer.setYwName("overlay_result");
    resultLayer.setWkid(layer1.getWkid());
    resultLayer.setGeometryType(GeometryType.POLYGON);
    
    // 合并两个图层的字段
    List<WktField> fields = new ArrayList<>();
    fields.addAll(layer1.getFields());
    for (WktField field : layer2.getFields()) {
        field.setYwName("L2_" + field.getYwName());  // 避免字段名冲突
        fields.add(field);
    }
    resultLayer.setFields(fields);
    
    List<WktFeature> features = new ArrayList<>();
    
    for (WktFeature f1 : layer1.getFeatures()) {
        Geometry g1 = GeometryConverter.wkt2Geometry(f1.getWkt());
        
        for (WktFeature f2 : layer2.getFeatures()) {
            Geometry g2 = GeometryConverter.wkt2Geometry(f2.getWkt());
            
            // 判断是否相交
            if (g1.intersects(g2)) {
                // 计算交集
                Geometry intersection = g1.intersection(g2);
                
                if (!intersection.isEmpty() && intersection.getArea() > 0) {
                    WktFeature feature = new WktFeature();
                    feature.setWfId(f1.getWfId() + "_" + f2.getWfId());
                    feature.setWkt(intersection.toText());
                    
                    // 合并属性值
                    List<WktFieldValue> values = new ArrayList<>();
                    values.addAll(f1.getFieldValues());
                    values.addAll(f2.getFieldValues());
                    feature.setFieldValues(values);
                    
                    features.add(feature);
                }
            }
        }
    }
    
    resultLayer.setFeatures(features);
    return resultLayer;
}

案例3：地块合并

合并相邻的地块：

/**
 * 合并相邻地块
 */
public String mergeAdjacentParcels(List<String> wkts) {
    List<Geometry> geoms = wkts.stream()
        .map(GeometryConverter::wkt2Geometry)
        .collect(Collectors.toList());
    
    // 使用并集操作合并
    Geometry merged = union(geoms.toArray(new Geometry[0]));
    
    return GeometryConverter.geometry2Wkt(merged);
}

案例4：区域裁剪

按行政区边界裁剪数据：

/**
 * 区域裁剪
 */
public WktLayer clipByBoundary(WktLayer dataLayer, String boundaryWkt) {
    Geometry boundary = GeometryConverter.wkt2Geometry(boundaryWkt);
    
    WktLayer resultLayer = new WktLayer();
    resultLayer.setYwName(dataLayer.getYwName() + "_clipped");
    resultLayer.setWkid(dataLayer.getWkid());
    resultLayer.setGeometryType(dataLayer.getGeometryType());
    resultLayer.setFields(dataLayer.getFields());
    
    List<WktFeature> features = new ArrayList<>();
    
    for (WktFeature feature : dataLayer.getFeatures()) {
        Geometry geom = GeometryConverter.wkt2Geometry(feature.getWkt());
        
        if (geom.intersects(boundary)) {
            Geometry clipped = geom.intersection(boundary);
            
            if (!clipped.isEmpty()) {
                WktFeature clippedFeature = new WktFeature();
                clippedFeature.setWfId(feature.getWfId());
                clippedFeature.setWkt(clipped.toText());
                clippedFeature.setFieldValues(feature.getFieldValues());
                features.add(clippedFeature);
            }
        }
    }
    
    resultLayer.setFeatures(features);
    return resultLayer;
}

案例5：多比例尺数据生成

为不同显示比例尺生成不同详细程度的数据：

/**
 * 生成多比例尺数据
 */
public Map<String, WktLayer> generateMultiScaleData(WktLayer sourceLayer) {
    Map<String, WktLayer> result = new HashMap<>();
    
    // 原始数据
    result.put("1:1000", sourceLayer);
    
    // 简化程度：根据比例尺确定容差
    double[] tolerances = {0.1, 0.5, 1.0, 5.0};
    String[] scales = {"1:5000", "1:10000", "1:50000", "1:100000"};
    
    for (int i = 0; i < tolerances.length; i++) {
        WktLayer simplifiedLayer = simplifyLayer(sourceLayer, tolerances[i]);
        simplifiedLayer.setYwName(sourceLayer.getYwName() + "_" + scales[i]);
        result.put(scales[i], simplifiedLayer);
    }
    
    return result;
}

private WktLayer simplifyLayer(WktLayer layer, double tolerance) {
    WktLayer result = new WktLayer();
    result.setYwName(layer.getYwName());
    result.setWkid(layer.getWkid());
    result.setGeometryType(layer.getGeometryType());
    result.setFields(layer.getFields());
    
    List<WktFeature> features = new ArrayList<>();
    for (WktFeature feature : layer.getFeatures()) {
        Geometry geom = GeometryConverter.wkt2Geometry(feature.getWkt());
        Geometry simplified = simplify(geom, tolerance);
        
        WktFeature simplifiedFeature = new WktFeature();
        simplifiedFeature.setWfId(feature.getWfId());
        simplifiedFeature.setWkt(simplified.toText());
        simplifiedFeature.setFieldValues(feature.getFieldValues());
        features.add(simplifiedFeature);
    }
    
    result.setFeatures(features);
    return result;
}

小结

本章介绍了空间计算操作的核心内容：

1. 缓冲区分析：生成围绕几何的指定距离区域
2. 叠加分析：交集、并集、差集、对称差集
3. 凸包与凹包：获取几何的包围边界
4. 几何简化：减少节点数量保持形状
5. 多边形化：从线生成面
6. 切割操作：用线切割多边形

下一章将介绍Shapefile文件的处理方法。