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

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

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

具体流程

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

注意点

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构建字符串，构建Geometry
        Geometry 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
            ]
        ]
    ]
}

影像文件示例