Java调用GDAL根据指定坐标点读取栅格数据像元值

目录

一、前言

二、实现思路 

三、实现的代码

1、程序代码

2、运行结果

3、验证

四、补充说明

一、前言

在3S技术（全球卫星定位系统、地理信息系统、遥感）中，栅格数据是一类常见的数据，栅格数据有其特定的数据结构。简单点说，栅格数据的表现的信息主要体现在每个像元对应的数值上，可以抽象为数学上的“矩阵”，这一点与数字图像类似，与其不同的是栅格数据往往带有空间信息（地理坐标、投影坐标），所以根据坐标点读取栅格数据的值是一项基础的工作。

GDAL是开源的矢量、栅格数据的处理程序，支持较多的编程语言，但网络上有关Java的帖子相对C和Python来说很少。所以在这里记录一下，已备不时之需。

二、实现思路 

根据GDAL的API文档，解析的思路应该是，首先根据文件路径获取到对应的DataSet类，然后从DataSet类获取Band类，然后调用Band类的ReadRaster获取一个一维数组，这个一维数组就是该栅格数据的某一区域的值。当ReadRaster的参数范围足够小时，就可以获取到一个像元的值，即一维数组只有一个值。

由于ReadRaster的参数需要像素坐标来定位，所以需要先使用DataSet类的GetGeoTransform方法获取像素坐标转地图坐标的参数，然后根据像素坐标转地图坐标的公式和指定的地图坐标点，解二元一次方程组，得到像素坐标。从而调用ReadRaster方法。

三、实现的代码

1、程序代码

package com.myself.raster;

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

import java.util.Arrays;

public class RasterDemo {

    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);
        }
        //获取驱动
        Driver driver = dataset.GetDriver();
        //获取驱动信息
        System.out.println("driver long name: " + driver.getLongName());
        //获取栅格数量
        int bandCount = dataset.getRasterCount();
        System.out.println("RasterCount: " + bandCount);
        //构造仿射变换参数数组，并获取数据
        double[] gt = new double[6];
        dataset.GetGeoTransform(gt);
        System.out.println("仿射变换参数"+Arrays.toString(gt));

        //指定经纬度
        double Latitude  = 86.053;
        double longitude = 16.529;

        //经纬度转换为栅格像素坐标
        int[] ColRow=Coordinates2ColRow(gt,Latitude,longitude);

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

        //遍历波段，获取该点对应的每个波段的值并打印到屏幕
        for (int i = 0;i<bandCount;i++){
            Band band = dataset.GetRasterBand(i+1);
            double[] values = new double[1];
            band.ReadRaster(ColRow[0], ColRow[1], 1, 1, values);
            double value = values[0];
            System.out.println("横坐标："+Latitude +","+"纵坐标:"+longitude);
            System.out.format("Band"+(i+1)+": %s", value);
        }
        //释放资源
        dataset.delete();

    }

    /**
     * 将地图坐标转换为栅格像素坐标
     * @param gt 仿射变换参数
     * @param X 横坐标
     * @param Y 纵坐标
     * @return
     */
    public static int[] Coordinates2ColRow(double[] gt, double X, double Y){
        int[] ints = new int[2];
//        X = gt[0] + Xpixel*gt[1] + Yline*gt[2];
//        Y = gt[3] + Xpixel*gt[4] + Yline*gt[5];
//        消元法解二元一次方程组
//        X-gt[0]=Xpixel*gt[1] + Yline*gt[2]
//        Xpixel = (X-gt[0] - Yline*gt[2])/gt[1]
//        Y - gt[3] = ((X-gt[0] - Yline*gt[2])/gt[1])gt[4] + Yline*gt[5]
//        (Y - gt[3])*gt[1] = ((X-gt[0] - Yline*gt[2]))*gt[4] + Yline*gt[5]*gt[1]
//        (Y - gt[3])*gt[1] =(X-gt[0])*gt[4] - Yline*gt[2]*gt[4] + Yline*gt[5]*gt[1]
//        (Y - gt[3])*gt[1] - (X-gt[0])*gt[4] = Yline(gt[5]*gt[1]-gt[2]*gt[4])

        //向下取整,如果向上取整会导致计算结果偏大，从而在后面读取到邻近像元的数据
        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;
    }

}
           

2、 运行结果

3、验证

这里使用QGIS来打开数据文件，并使用识别功能来验证有关点的值。

 可以发现，程序解析结果与QGIS解析结果一致。

四、补充说明

dataset.GetGeoTransform获取的gt[]含义如下：

//gt[0]  左上角x坐标

//gt[1]  东西方向分辨率

//gt[2]  旋转角度1, 0表示图像 "北方朝上"

//gt[3]  左上角y坐标

//gt[4]  旋转角度2, 0表示图像 "北方朝上"

//gt[5]  东西方向分辨率

地图坐标X =  左上角x坐标 + 列数*东西方向分辨率 + 行数*旋转角度1；

地图坐标Y =  左上角y坐标 + 列数*旋转角度2 + 行数*南北方向分辨率；

band.ReadRaster的入参含义如下：public int ReadRaster(int xoff, int yoff, int xsize, int ysize, double[] array)

public int ReadRaster(int xoff, int yoff, int xsize, int ysize, double[] array)

xoff 想要读取的部分原点位置横图像坐标

yoff 想要读取的部分原点位置纵图像坐标

xsize 指定要读取部分图像的矩形边长。 

ysize 指定要读取部分图像的矩形边长。 

array 用来存储读取到的数据的数组。

本文的数据是作者自己做的，没有任何业务意义。如有需要可以下载：Demo.tif