Java调用GDAL读取指定点附近指定距离方格内的栅格像元值

一、问题描述

最近遇到了需要根据指定的坐标点读取某个栅格文件对应位置附近5公里方格内的最大像元值的问题。

二、使用Java调用GDAL解决的方法

（一）实现思路

根据ReadRaster的调用方法启示：

根据坐标点和指定的距离计算出ReadRaster的调用参数，然后调用ReadRaster的方法读取。这种方法读取的区域从平面上看是“矩形”。

   public int ReadRaster(int xoff, int yoff, int xsize, int ysize, double[] array) {
       return ReadRaster(xoff, yoff, xsize, ysize, xsize, ysize, gdalconstConstants.GDT_Float64, array);
   }

 （二）实现代码

以下代码只是演示，读取的区域是个“矩形”。

package com.myself.raster;

import org.gdal.gdal.Band;
import org.gdal.gdal.Dataset;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconstConstants;

import java.util.Arrays;

public class ReadRange {
    public static void main(String[] args) {
        String filepath = "F:\\raster\\Demo.tif";
        //注册驱动
        gdal.AllRegister();
        //打开文件获取数据集
        Dataset dataset = gdal.Open(filepath,
                gdalconstConstants.GA_ReadOnly);
        if (dataset == null) {
            System.out.println("打开" + filepath + "失败" + gdal.GetLastErrorMsg());
            System.exit(1);
        }

        //构造仿射变换参数数组，并获取数据
        double[] gt = new double[6];
        dataset.GetGeoTransform(gt);

        //指定中心点经纬度
        double Longitude = 86.053;
        double Latitude = 16.529;

        //获取东西方向空间分辨率
        System.out.println("东西方向分辨率："+gt[1]+","+"南北方向分辨率："+gt[5]);

        //指定距离0.23度
        double distance = 0.23;

        //计算起点经纬度
        double startLon = Longitude - distance/2;
        double startLat = Latitude + distance/2;

        //计算终点经纬度
        double endLon = Longitude + distance/2;
        double endLat = Latitude - distance/2;

        //经纬度转换为栅格像素坐标
        int[] ColRow = Coordinates2ColRow(gt,startLon,startLat);
        int[] eColRow = Coordinates2ColRow(gt,endLon,endLat);

        //判断是否坐标超出范围
        if (ColRow[0] < 0 || ColRow[1] < 0 || ColRow[0] > dataset.getRasterXSize() || ColRow[1] > dataset.getRasterYSize()) {
            System.out.println(Arrays.toString(ColRow) + "坐标值超出栅格范围！");
            return;
        }else if (eColRow[0] < 0 || eColRow[1] < 0 || eColRow[0] > dataset.getRasterXSize() || eColRow[1] > dataset.getRasterYSize()){
            System.out.println(Arrays.toString(ColRow) + "坐标值超出栅格范围！");
            return;
        }

        //遍历波段，获取该点对应的每个波段的值并打印到屏幕
        for (int i = 0; i < dataset.getRasterCount(); i++) {
            Band band = dataset.GetRasterBand(i + 1);
            //计算xsize与ysize
            int subX = eColRow[0] - ColRow[0];
            int subY = eColRow[1] - ColRow[1];
            double[] values = new double[subX * subY];
            band.ReadRaster(ColRow[0], ColRow[1],subX,subY, values);
            System.out.println("横坐标：" + Longitude + "," + "纵坐标:" + Latitude);
            System.out.println("数组："+Arrays.toString(values));
        }
        //释放资源
        dataset.delete();
    }

    /**
     * 将地图坐标转换为栅格像素坐标
     *
     * @param gt 仿射变换参数
     * @param X  横坐标
     * @param Y  纵坐标
     * @return
     */
    public static int[] Coordinates2ColRow(double[] gt, double X, double Y) {
        int[] ints = new int[2];

        //向下取整,如果向上取整会导致计算结果偏大，从而在后面读取到邻近像元的数据
        double Yline = Math.floor(((Y - gt[3]) * gt[1] - (X - gt[0]) * gt[4]) / (gt[5] * gt[1] - gt[2] * gt[4]));
        double Xpixel = Math.floor((X - gt[0] - Yline * gt[2]) / gt[1]);
        ints[0] = new Double(Xpixel).intValue();
        ints[1] = new Double(Yline).intValue();
        return ints;
    }
}

 控制台输出

根据读取的数组再去进行求最大值就解决了问题。