Java开发环境中，使用GDAL进行矢量叠加，并计算面积

GDAL与Springboot的集成可参考Springboot 集成GDAL开发环境配置

Java开发环境下，GDAL的相关学习和使用的案例还是非常少的，并且部分函数的使用方式和Python、C环境下有很大的区别。已最近做的一个功能为案例，给大家分享一下Java环境下GDAL的用法。

具体流程

有一个已经做好的某种作物的种植情况遥感解译栅格影像（单波段），需要提供一个接口，入参为目标区域面的空间坐标，最后返回结果为目标区域内属于该作物的种植面积。实现流程如下：

#mermaid-svg-MOYkgoALeTv10ufs {font-family:"trebuchet ms",verdana,arial,sans-serif;font-size:16px;fill:#333;}#mermaid-svg-MOYkgoALeTv10ufs .error-icon{fill:#552222;}#mermaid-svg-MOYkgoALeTv10ufs .error-text{fill:#552222;stroke:#552222;}#mermaid-svg-MOYkgoALeTv10ufs .edge-thickness-normal{stroke-width:2px;}#mermaid-svg-MOYkgoALeTv10ufs .edge-thickness-thick{stroke-width:3.5px;}#mermaid-svg-MOYkgoALeTv10ufs .edge-pattern-solid{stroke-dasharray:0;}#mermaid-svg-MOYkgoALeTv10ufs .edge-pattern-dashed{stroke-dasharray:3;}#mermaid-svg-MOYkgoALeTv10ufs .edge-pattern-dotted{stroke-dasharray:2;}#mermaid-svg-MOYkgoALeTv10ufs .marker{fill:#333333;stroke:#333333;}#mermaid-svg-MOYkgoALeTv10ufs .marker.cross{stroke:#333333;}#mermaid-svg-MOYkgoALeTv10ufs svg{font-family:"trebuchet ms",verdana,arial,sans-serif;font-size:16px;}#mermaid-svg-MOYkgoALeTv10ufs .label{font-family:"trebuchet ms",verdana,arial,sans-serif;color:#333;}#mermaid-svg-MOYkgoALeTv10ufs .cluster-label text{fill:#333;}#mermaid-svg-MOYkgoALeTv10ufs .cluster-label span{color:#333;}#mermaid-svg-MOYkgoALeTv10ufs .label text,#mermaid-svg-MOYkgoALeTv10ufs span{fill:#333;color:#333;}#mermaid-svg-MOYkgoALeTv10ufs .node rect,#mermaid-svg-MOYkgoALeTv10ufs .node circle,#mermaid-svg-MOYkgoALeTv10ufs .node ellipse,#mermaid-svg-MOYkgoALeTv10ufs .node polygon,#mermaid-svg-MOYkgoALeTv10ufs .node path{fill:#ECECFF;stroke:#9370DB;stroke-width:1px;}#mermaid-svg-MOYkgoALeTv10ufs .node .label{text-align:center;}#mermaid-svg-MOYkgoALeTv10ufs .node.clickable{cursor:pointer;}#mermaid-svg-MOYkgoALeTv10ufs .arrowheadPath{fill:#333333;}#mermaid-svg-MOYkgoALeTv10ufs .edgePath .path{stroke:#333333;stroke-width:2.0px;}#mermaid-svg-MOYkgoALeTv10ufs .flowchart-link{stroke:#333333;fill:none;}#mermaid-svg-MOYkgoALeTv10ufs .edgeLabel{background-color:#e8e8e8;text-align:center;}#mermaid-svg-MOYkgoALeTv10ufs .edgeLabel rect{opacity:0.5;background-color:#e8e8e8;fill:#e8e8e8;}#mermaid-svg-MOYkgoALeTv10ufs .cluster rect{fill:#ffffde;stroke:#aaaa33;stroke-width:1px;}#mermaid-svg-MOYkgoALeTv10ufs .cluster text{fill:#333;}#mermaid-svg-MOYkgoALeTv10ufs .cluster span{color:#333;}#mermaid-svg-MOYkgoALeTv10ufs div.mermaidTooltip{position:absolute;text-align:center;max-width:200px;padding:2px;font-family:"trebuchet ms",verdana,arial,sans-serif;font-size:12px;background:hsl(80, 100%, 96.2745098039%);border:1px solid #aaaa33;border-radius:2px;pointer-events:none;z-index:100;}#mermaid-svg-MOYkgoALeTv10ufs :root{--mermaid-font-family:"trebuchet ms",verdana,arial,sans-serif;}

开始

接收GeoJson

保存为本地shp文件

shp与栅格影像叠加

计算有效栅格数量

计算面积

注意点

1.矢栅叠加

矢栅叠加的结果最好设置为EPSG:3857坐标系，它是以米为单位的坐标系，计算面积准确。

2.面积计算

面积计算有两种方式：

计算有效像元数量，乘以像元面积。像元面积为分辨率的平方。
叠加后的栅格转矢量，计算矢量面积。（只适合单波段）。

前面的方法相对简单，我这里采用的是第一种方法。

3.GeoJson转Shp文件

GeoJson转Shp文件时一定要设置坐标系。不然后面的叠加可能出不来结果。

全部代码

package com.example.gdal.controller;import org.gdal.gdal.*;
import org.gdal.gdalconst.gdalconstConstants;
import org.gdal.ogr.*;
import org.gdal.ogr.Driver;
import org.gdal.osr.SpatialReference;
import org.springframework.web.bind.annotation.PostMapping;
import org.springframework.web.bind.annotation.RequestBody;
import org.springframework.web.bind.annotation.RequestMapping;
import org.springframework.web.bind.annotation.RestController;import java.io.File;
import java.util.Arrays;
import java.util.Vector;/*** 调整区划代码*/
@RequestMapping("/")
@RestController
public class TestController {@PostMapping("/test")String addInsureTask(@RequestBody String geoJson) {ogr.RegisterAll();Dataset raster = gdal.Open("E:\\测试数据\\zzjg_41_16_2022_20220812\\栅格\\zzjg_41_16_2022_20220812.tif", gdalconstConstants.GA_ReadOnly);
//        //获取到六参数数组
//        double[] ori_transform = raster.GetGeoTransform();
//        //最小x
//        double minX = ori_transform[0];
//        //x方向分辨率
//        double resoutionX = ori_transform[1];
//        //x方向旋转角度
//        double angleX = ori_transform[2];
//        //最大y
//        double maxY = ori_transform[3];
//        //Y方向旋转角度
//        double angleY = ori_transform[4];
//        //y方向分辨率
//        double resoutionY = ori_transform[5];
//        //x方向的像素数
//        int pixelX = raster.getRasterXSize();
//        //y方向像素
//        int pixelY = raster.getRasterYSize();String tempShpPath="E:\\测试数据\\zzjg_41_16_2022_20220812\\栅格\\mask.shp";jsonToShp(geoJson,tempShpPath,raster.GetProjection());String resultPath="E:\\测试数据\\zzjg_41_16_2022_20220812\\栅格\\111.tiff";warp(resultPath,raster,tempShpPath);double area=getArea(resultPath);gdal.GDALDestroyDriverManager();return area+"平方米";}/*** 矢量裁剪栅格* @param targetPath* @param raster* @param shpPath*/public void warp(String targetPath,Dataset raster,String shpPath){//        SpatialReference oSRS = new SpatialReference();
//        oSRS.ImportFromEPSGA(3857);
//        System.out.println(oSRS.toString());
//        oSRS.SetWellKnownGeogCS("EPSG:3857");Vector<String> vector=new Vector<>();vector.addElement("-t_srs");vector.addElement("EPSG:3857");//结果转到3857坐标系，便于计算面积vector.addElement("-cutline");vector.addElement(shpPath);vector.addElement("-crop_to_cutline");vector.addElement("true");WarpOptions warpOptions=new WarpOptions(vector);Dataset[] sources=new Dataset[]{raster};//矢栅叠加gdal.Warp(targetPath,sources,warpOptions);}/*** geoJson转shp文件* @param geoJson* @param tempShpPath* @param srs*/public void jsonToShp(String geoJson, String tempShpPath, String srs){//从geoJson构建字符串，构建GeometryGeometry geometry=ogr.CreateGeometryFromJson(geoJson);Driver shpDriver=ogr.GetDriverByName("ESRI Shapefile");if(new File(tempShpPath).exists()){//文件存在，提前删除shpDriver.DeleteDataSource(tempShpPath);}System.out.println(srs.toString());//创建shp文件DataSource tempShp =shpDriver.CreateDataSource(tempShpPath);SpatialReference spatialReference=new SpatialReference(srs);//复制坐标系Layer layer=tempShp.CreateLayer("polygon",spatialReference);Feature feature = new Feature(layer.GetLayerDefn());feature.SetGeometry(geometry);layer.CreateFeature(feature);layer.SyncToDisk();tempShp.FlushCache();}public double getArea(String rasterPath){Dataset raster = gdal.Open(rasterPath, gdalconstConstants.GA_ReadOnly);int xSize=raster.getRasterXSize();int ySize=raster.getRasterYSize();Band band=raster.GetRasterBand(1);//读取noData值Double[] noData=new Double[1];band.GetNoDataValue(noData);double[] pixValue = new double[xSize * ySize];band.ReadRaster(0,0,xSize,ySize,pixValue);long count=0;for (int i = 0; i < pixValue.length; i++) {if(pixValue[i]!=noData[0]){count++;}}double pixel = raster.GetGeoTransform()[1];double area=pixel*pixel*count;return area;}}

参数示例

{"type": "Polygon","coordinates": [[[114.36676024162932,33.562408688929679],[114.36640332247839,33.556757469040804],[114.37135102442949,33.556444982601761],[114.37170794358042,33.562096202490636],[114.36676024162932,33.562408688929679]]]
}

影像文件示例