gdal统计影像矢量切割后的土地利用面积

1. 本次使用的影像是某个地区的土地利用影像，该影像已经处理过，7个像素分别代表7种土地利用。
   分两步，第一步使用某县乡的shp文件切割影像，获取该县乡的土地利用影像（该县乡shp是影像区域的子集）；
   第二步统计切割后的影像，各土地种类像元个数和每个像元的面积，进而获取各土地利用的面积。
2. 程序代码如下

package com.aerors.etl.utils;

import java.util.HashMap;
import java.util.Iterator;
import java.util.Map;
import java.util.Map.Entry;
import java.util.Vector;

import org.gdal.gdal.Dataset;
import org.gdal.gdal.WarpOptions;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconst;
import org.gdal.ogr.ogr;

public class RasterClip {

    public static void main(String[] args) {
        gdal.AllRegister();
        ogr.RegisterAll();

        String inTiff = "E:\\data_sjy\\sjy_landuse\\qinghai2017";
        String inShp = "E:\\data_sjy\\xiang_shp\\chengduoxianxiang.shp";
        Dataset src = gdal.Open(inTiff);
        //单个像元宽度，高度和宽度一致
        double[] geotran = src.GetGeoTransform();
        System.out.println("geotran[1]:"+geotran[1]);

        Dataset[] src_array = { src };

        Vector<String> options = new Vector<>();
        options.add("-of");
        options.add("GTiff");
        options.add("-cutline");
        options.add(inShp);
        options.add("-crop_to_cutline");
        options.add("-dstalpha");
        WarpOptions warpOptions = new WarpOptions(options);

        Dataset warp = gdal.Warp("warp.vrt", src_array, warpOptions);
        
        long startTime = System.currentTimeMillis();

        //统计
        int bandCount = warp.getRasterCount();
        System.out.println("bandCount:"+bandCount);
        int width = warp.GetRasterBand(1).getXSize();
        int height = warp.GetRasterBand(1).getYSize();
        
        int blockX = warp.GetRasterBand(1).GetBlockXSize();
        int blockY = warp.GetRasterBand(1).GetBlockYSize();
        
        System.out.println("blockX:" + blockX);
        System.out.println("blockY:" + blockY);
                //map存放土地类别和像元总个数
        Map<Integer, Integer> map = new HashMap<Integer, Integer>();
        //逐行扫描统计
        for(int i=0;i<height;i++) {
            int[] lineBuffer = new int[width * bandCount];
            warp.ReadRaster(0, 0, width, 1, width, 1,gdalconst.GDT_Int32, lineBuffer, new int[] {1});
            
            //System.out.println(lineBuffer.length);
            int lbefore = lineBuffer[0];
            for(int j=0; j<lineBuffer.length;j++) {
                if(lineBuffer[j] != lbefore) {
                    if(lineBuffer[j] != 0) {
                        if(map.containsKey(lineBuffer[j])){
                            map.put(lineBuffer[j], map.get(lineBuffer[j])+1);
                        }else {
                            map.put(lineBuffer[j], 1);
                        }
                    }
                }
                lbefore = lineBuffer[j];
            }
        }
                //map存放土地类别和土地面积
        Map<Integer, Double> areaMap = new HashMap<Integer, Double>();
        Iterator<Entry<Integer, Integer>> mi = map.entrySet().iterator();
        while(mi.hasNext()) {
            Entry<Integer, Integer> element = mi.next();
            double area = element.getValue()*geotran[1]*geotran[1];
            areaMap.put(element.getKey(), area);
            System.out.println("reulst area: "+element.getKey() +": "+area);
        }
        long endTime = System.currentTimeMillis();
        System.out.println("dur:"+(endTime-startTime));

//      Driver tiffDriver = gdal.GetDriverByName("GTiff");
//      Dataset warp_out = tiffDriver.CreateCopy(outTiff, warp);

        src.delete();
        warp.delete();
//      warp_out.delete();
    }
}