GDAL矢量图层裁剪为矢量网格（java）

1、矢量图层

打开矢量图层，并显示基本信息。

图层数，驱动名，图层网格范围，属性表结构信息，矢量元素，及矢量元素的属性

public static void info(String srcFile) {
     // 读取原图
     DataSource srcDataSource = ogr.Open(srcFile, 0);
     if (srcDataSource == null) {
         System.err.println("GDALOpen failed - " + gdal.GetLastErrorNo());
         System.err.println(gdal.GetLastErrorMsg());
         return;
     }

     // 获取该数据源中的图层个数，一般shp数据图层只有一个，如果是mdb、dxf等图层就会有多个
     int layerCount = srcDataSource.GetLayerCount();
     System.out.println("图层数 = " + layerCount);
     if (layerCount < 1) {
         System.err.println("异常：没有图层。");
         return;
     }
     System.out.println(srcDataSource.GetDriver().GetName());

     Layer layer = srcDataSource.GetLayer(0);
     // 对图层进行初始化，如果对图层进行了过滤操作，执行这句后，之前的过滤全部清空
     layer.ResetReading();
     // 通过属性表的SQL语句对图层中的要素进行筛选，这部分详细参考SQL查询章节内容
     // layer.SetAttributeFilter("\"NAME99\"LIKE \"北京市市辖区\"");
     // 通过指定的几何对象对图层中的要素进行筛选
     // layer.SetSpatialFilter();
     // 通过指定的四值范围对图层中的要素进行筛选
     // layer.SetSpatialFilterRect();

     double[] extent = layer.GetExtent(true);
     System.out.println("extent is " + extent[0] + ", " + extent[1] + ", " + extent[2] + ", " + extent[3]);

     // 获取图层中的属性表表头并输出
     System.out.println("属性表结构信息：");
     FeatureDefn featureDefn = layer.GetLayerDefn();
     int fieldCount = featureDefn.GetFieldCount();
     for (int index = 0; index < fieldCount; index++) {
         FieldDefn fieldDefn = featureDefn.GetFieldDefn(index);
         System.out.println(fieldDefn.GetNameRef() + ":" +
                 fieldDefn.GetFieldTypeName(fieldDefn.GetFieldType()) + "(" +
                 fieldDefn.GetWidth() + "." +
                 fieldDefn.GetPrecision() + ")");
     }
     // 输出图层中的要素个数
     System.out.println("要素个数 = " + layer.GetFeatureCount(0));

     StringBuffer sb = new StringBuffer();
     Feature feature = layer.GetNextFeature();
     while (feature != null) {
         System.out.println("当前处理第：" + feature.GetFID());
         // 获取要素中的属性表内容
         for (int index = 0; index < fieldCount; index++) {
             FieldDefn fieldDefn = featureDefn.GetFieldDefn(index);
             sb.append(fieldDefn.GetNameRef());
             sb.append("(");
             sb.append(ogr.GetFieldTypeName(fieldDefn.GetFieldType()));
             sb.append(") = ");
             if (feature.IsFieldSet(index)) {
                 sb.append(feature.GetFieldAsString(index));
             } else {
                 sb.append("(null)");
             }
             System.out.println(sb.toString());
         }
         // 获取要素中的几何体
         Geometry geometry = feature.GetGeometryRef();
         feature = layer.GetNextFeature();
     }
     srcDataSource.delete();
 }

2、根据图层的网格范围（extent）创建网格

使用【ESRI Shapefile】矢量驱动，
建立指定投影【4326】的图层，
根据图层的范围，及水平间隔、垂直间隔 创建矩形矢量图层。

注意：矩形的wkt为 POLYGON ((30 0,60 0,60 30,30 30,30 0)) 五个点

private static void writeMesh(String destFile, double startX, double endX, double startY, double endY, float interX, float interY) {

     // 创建数据，这里以创建ESRI的shp文件为例
     Driver driver = ogr.GetDriverByName("ESRI Shapefile");
     if (driver == null) {
         System.err.println("驱动不可用！");
         return;
     }

     //创建数据源
     DataSource dataSource = driver.CreateDataSource(destFile);
     if (dataSource == null) {
         System.err.println("创建文件【" + dataSource + "】失败！");
         return;
     }

     SpatialReference sr = new SpatialReference();
     sr.ImportFromEPSG(4326);

     // 创建图层，创建一个多边形图层，这里没有指定空间参考，如果需要的话，需要在这里进行指定
     Layer layer = dataSource.CreateLayer("TestPolygon", sr);
     if (layer == null) {
         System.err.println("层创建失败！");
         return;
     }

     // 下面创建属性表
     // 先创建一个叫FieldID的整型属性
     FieldDefn oFieldID = new FieldDefn("id", ogr.OFTInteger);
     layer.CreateField(oFieldID, 1);
     layer.CreateField(new FieldDefn("left", ogr.OFTString), 1);
     layer.CreateField(new FieldDefn("top", ogr.OFTString), 1);
     layer.CreateField(new FieldDefn("right", ogr.OFTString), 1);
     layer.CreateField(new FieldDefn("bottom", ogr.OFTString), 1);

     FeatureDefn featureDefn = layer.GetLayerDefn();

     int index = 1;
     // 创建矩形要素
     for (double x = startX; x < endX; x = x + interX) {
         for (double y = startY; y < endY; y = y + interY) {
             StringBuffer sb = new StringBuffer();
             sb.setLength(0);
             sb.append("POLYGON (");

             sb.append("(");
             sb.append(x).append(" ").append(y).append(",");
             sb.append(x + interX).append(" ").append(y).append(",");
             sb.append(x + interX).append(" ").append(y + interY).append(",");
             sb.append(x).append(" ").append(y + interY).append(",");
             sb.append(x).append(" ").append(y);
             sb.append("))");

             System.out.println("wkt = " + sb.toString());
             Geometry geomRectangle = ogr.CreateGeometryFromWkt(sb.toString(), sr);
             Feature feature = new Feature(featureDefn);
             feature.SetField(0, index++);
             feature.SetField(1, String.valueOf(x));
             feature.SetField(2, String.valueOf(y));
             feature.SetField(3, String.valueOf(x + interX));
             feature.SetField(4, String.valueOf(y + interY));

             feature.SetGeometry(geomRectangle);
             layer.CreateFeature(feature);
         }
     }
     dataSource.delete();
 }

3、使用网格和最初矢量图层求交集

核心是直接调用layer的Intersection方法。
类似还有，裁剪，差分，联合等方法。

public static void Intersection(String srcFile, String targetFile, String dstFile) {
     DataSource srcDataSource = ogr.Open(srcFile, 0);
     if (srcDataSource == null) {
         System.err.println("GDALOpen failed - " + gdal.GetLastErrorNo());
         System.err.println(gdal.GetLastErrorMsg());
         return;
     }

     DataSource tgtDataSource = ogr.Open(targetFile, 0);
     if (tgtDataSource == null) {
         System.err.println("GDALOpen failed - " + gdal.GetLastErrorNo());
         System.err.println(gdal.GetLastErrorMsg());
         return;
     }

     Driver dstDriver = ogr.GetDriverByName("ESRI Shapefile");
     DataSource dstDataSource = dstDriver.CreateDataSource(dstFile);
     if (dstDataSource == null) {
         System.err.println("创建文件 失败！");
         return;
     }
     Layer srcLayer = srcDataSource.GetLayer(0);
     Layer tgtLayer = tgtDataSource.GetLayer(0);

     SpatialReference sr = new SpatialReference();
     sr.ImportFromEPSG(4326);
     Layer dstLayer = dstDataSource.CreateLayer("Intersection", sr);

     srcLayer.Intersection(tgtLayer, dstLayer);

     tgtDataSource.delete();
     dstDataSource.delete();
     srcDataSource.delete();
 }

4、最终效果

1. 原图
   
2. 创建网格
   
3. 矢量图像相交