Java调用GDAL实现融合有相同字段属性的多边形矢量数据的一种方法

目录

一、写在文章前

二、实现思路

三、实现过程

1.打开矢量数据

 2.生成融合结果

 3.不融合不相邻的几何图形

4.融合后的几何图形写入到新的图层

5.融合的效果

四、完整的示例 

---

一、写在文章前

如下图所示，融合（dissolve）具有相同字段属性的多边形矢量数据是日常的GIS工作中经常会用到功能，它在数据分析、制图工作中都有重要的作用，目前的GIS软件中均有此功能。个人认为在网络地理信息系统的开发中，使用GDAL是一种性价比较高的方式。据作者所知，GDAL暂时还没有提供融合（dissolve）功能的API，所以作者自己编写了一个用Java调用GDAL实现融合功能的示例程序，在这里共享出来供大家参考、指正。

二、实现思路

一种思路是：

1. 打开矢量数据，将矢量图层所有的要素存放到一个集合中。
2. 两次遍历集合将属性值一致的的要素的几何图形（geometry）进行合并（union）操作，并存入融合结果的集合。
3. 若不希望将不相邻的几何图形（geometry）进行融合，可以遍历融合结果的集合，对于不是环（ring）的几何图形部件逐个将其存入到拆分多部件的集合中，对于是环的几何图形部件，不处理，将其对应的父几何图形直接存入到拆分多部件的集合中。
4. 将拆分多部件的集合中的要素逐个写入到新的图层，将新的图层存入新的数据中，保存。

这样就得到了融合后的矢量数据。

还有一种方法是对输入的矢量数据执行分组查询的SQL语句，对同一组的数据执行union操作，这样可以替代步骤2，从而减少代码量。GDAL支持的SQL语句语法与SQL99是基本一致的，但需要了解GDAL的空间操作的函数。

本文的实现使用了后者。

三、实现过程

本节主要展示效果及关键代码，完整实现过程在下文。

1.打开矢量数据

        ogr.RegisterAll();
        String dataSet = "E:\\grids.shp";
        String outPath = "E:\\grids_merge.shp";

注册了所有的驱动，指定了输入和输出文件。

输入的数据是一个用QGIS生成的网格，存在一个字段属性”grade“。几何图形如下图所示。

        DataSource ds = ogr.Open(dataSet);
        Layer layer = ds.GetLayerByName(dataName);

这是打开数据源和图层的代码。 

 2.生成融合结果

        //拼接SQL语句，严格一点的话，可能需要用JDBC先预编译SQL语句，这里简化了。
        String sql = "select ST_Union(GEOMETRY) AS GEOMETRY,grade,count(*) as count from " + dataName + " group by grade";
        // 执行SQL语句
        Layer layerSQL = ds.ExecuteSQL(sql, null, "SQLITE");
        //通过ExecuteSQL获取的Layer必须在销毁数据源之前使用ReleaseResultSet销毁。
        ds.ReleaseResultSet(layerSQL);

根据GDAL的API文档要求对于执行SQL获取的图层，必须要执行ReleaseResultSet销毁。销毁代码需要放在销毁数据源的代码之前。

 3.不融合不相邻的几何图形

如果需要不融合不相邻的几何图形，还需要下面的代码。

        Feature feature;
        Map<Geometry, Feature> map = new HashMap<>();
        //将执行SQL语句后的图层要素中的”多部件“几何图形分割为单部件
        while ((feature = layerSQL.GetNextFeature()) != null) {
            Geometry geometry = feature.GetGeometryRef();
            int count = geometry.GetGeometryCount();
            for (int i = 0; i < count; i++) {
                Geometry part = geometry.GetGeometryRef(i);
                boolean isRing = part.GetGeometryType() == ogrConstants.wkbLineString || part.GetGeometryType() == ogrConstants.wkbMultiLineString;

                if (isRing) {
                    map.put(geometry, feature);
                    break;
                } else {
                    map.put(part, feature);
                }
            }
        }

4.融合后的几何图形写入到新的图层

将融合后的几何图形写入新的图层，对于不需要不融合不相邻的几何图形的情况，只需要复制SQL执行结果的图层到新数据源即可。

DataSource outDataSource = driver.CreateDataSource(outPath);
outDataSource.CopyLayer(layerSQL,name,options);

对于需要不融合不相邻的几何图形，则需要执行下列代码。

        for (int i = 0; i < layerSQL.GetLayerDefn().GetFieldCount(); i++) {
            FieldDefn fieldDefn = layerSQL.GetLayerDefn().GetFieldDefn(i);
            layer1.CreateField(fieldDefn);
        }

        for (Map.Entry<Geometry, Feature> entry : map.entrySet()) {
            Feature feature1 = entry.getValue();
            Feature feature2 = new Feature(feature1.GetDefnRef());

            for (int j = 0; j < feature1.GetFieldCount(); j++) {
                //复制属性表中的数据，这里简写了，只考虑的短整形的一种情况
                //对shapefile而言，实际还有长整形、浮点型、双精度、文本、日期等数据类型
                feature2.SetField(j, feature1.GetFieldAsInteger(j));
            }
            Geometry geometry2 = entry.getKey();

            feature2.SetGeometry(geometry2);
            layer1.CreateFeature(feature2);
        }

保存数据的代码如下 

        ds.ReleaseResultSet(layerSQL);
        outDataSource.SyncToDisk();

        ds.delete();
        outDataSource.delete();

5.融合的效果

融合数据的几何图形如下图所示：

融合数据的几何图形属性表如下： 

四、完整的示例 

需要说明的是，本例子也只是简单的实现，距离在生产环境使用还有一定的差距，仅供学习使用。

import org.gdal.gdal.gdal;
import org.gdal.ogr.*;

import java.io.File;
import java.util.HashMap;
import java.util.Map;
import java.util.Vector;

public class GDALMergeByAttribute {
    public static void main(String[] args) {
        ogr.RegisterAll();

        String dataSet = "E:\\grids.shp";
        String outPath = "E:\\grids_merge.shp";

        System.out.printf("开始融合数据！输入文件：%s\n", dataSet);
        // 打开数据源
        DataSource ds = ogr.Open(dataSet);

        if (ds == null) {
            System.out.println("打开数据源失败");
            return;
        }
        //一般GIS软件生成shapefile的”图层名“和文件名一般是一致的
        File dataSetFile = new File(dataSet);
        String dataName = dataSetFile.getName();
        dataName = dataName.substring(0, dataName.lastIndexOf("."));

        Layer layer = ds.GetLayerByName(dataName);
        if (layer == null) {
            System.out.println("打开图层失败");
            return;
        }


        //拼接SQL语句，严格一点的话，可能需要用JDBC先预编译SQL语句，这里简化了。
        String sql = "select ST_Union(GEOMETRY) AS GEOMETRY,grade,count(*) as count from " + dataName + " group by grade";
        System.out.printf("SQL语句的内容为：%s\n", sql);
        // 执行SQL语句
        Layer layerSQL = ds.ExecuteSQL(sql, null, "SQLITE");
        if (0 != gdal.GetLastErrorNo()) {
            System.out.println("执行SQL语句出错！");
            return;
        }
        System.out.printf("执行SQL语句成功！返回数据行数：%s\n", layerSQL.GetFeatureCount());
        Driver driver = ds.GetDriver();
        Vector<String> options = new Vector<>();
        //防止属性表中文乱码
        options.add("ENCODING=UTF-8");


        File file = new File(outPath);
        String name = file.getName();
        name = name.substring(0, name.lastIndexOf("."));
        Feature feature;
        Map<Geometry, Feature> map = new HashMap<>();
        //将执行SQL语句后的图层要素中的”多部件“几何图形分割为单部件
        while ((feature = layerSQL.GetNextFeature()) != null) {
            Geometry geometry = feature.GetGeometryRef();
            int count = geometry.GetGeometryCount();
            for (int i = 0; i < count; i++) {
                Geometry part = geometry.GetGeometryRef(i);
                boolean isRing = part.GetGeometryType() == ogrConstants.wkbLineString || part.GetGeometryType() == ogrConstants.wkbMultiLineString;
                if (isRing) {
                    map.put(geometry, feature);
                    break;
                } else {
                    map.put(part, feature);
                }
            }
        }

        System.out.printf("融合后的要素：%s个，分割不相邻要素后要素有%s个\n", layerSQL.GetFeatureCount(), map.size());
        DataSource outDataSource = driver.CreateDataSource(outPath);
        Layer layer1 = outDataSource.CreateLayer(name, layerSQL.GetSpatialRef(), layer.GetGeomType(), options);


        for (int i = 0; i < layerSQL.GetLayerDefn().GetFieldCount(); i++) {
            FieldDefn fieldDefn = layerSQL.GetLayerDefn().GetFieldDefn(i);
            layer1.CreateField(fieldDefn);
        }
        System.out.printf("创建属性字段成功！字段数：%s\n", layer1.GetLayerDefn().GetFieldCount());

        for (Map.Entry<Geometry, Feature> entry : map.entrySet()) {
            Feature feature1 = entry.getValue();
            Feature feature2 = new Feature(feature1.GetDefnRef());

            for (int j = 0; j < feature1.GetFieldCount(); j++) {
                //复制属性表中的数据，这里简写了，只考虑的短整形的一种情况
                //对shapefile而言，实际还有长整形、浮点型、双精度、文本、日期等数据类型
                feature2.SetField(j, feature1.GetFieldAsInteger(j));
            }
            Geometry geometry2 = entry.getKey();

            feature2.SetGeometry(geometry2);
            layer1.CreateFeature(feature2);
        }
        System.out.printf("新图层写入要素：%s\n", layer1.GetFeatureCount());
        //通过ExecuteSQL获取的Layer必须在销毁数据源之前使用ReleaseResultSet销毁。
        ds.ReleaseResultSet(layerSQL);
        outDataSource.SyncToDisk();
        System.out.printf("融合数据成功！文件名：%s\n", outDataSource.GetName());

        ds.delete();
        outDataSource.delete();

    }

}

 正常情况下的输出：

开始融合数据！输入文件：E:\grids.shp
SQL语句的内容为：SELECT ST_Union(GEOMETRY) AS GEOMETRY,grade,count(*) as count FROM grids group by grade
执行SQL语句成功！返回数据行数：3
融合后的要素：3个，分割不相邻要素后要素有4个
创建属性字段成功！字段数：2
新图层写入要素：4
融合数据成功！文件名：E:\grids_merge.shp