GDAL+Java实现获取对应栅格影像经纬度对应的像素值

从前面的GDAL系列博文中，可以指导GDAL可以将栅格影像文件读出为对应的多维数组，可以读出每一个像素格对应的像素值。但如何根据经纬度直接读取像素值呢？博主从查阅了网上的相关文档，发现有个人写的计算公式是错误的，用代码跑出来的结果都是错误的。于是自己查阅了相关文档，自己实现了一遍，跟大家分享一下。

图像坐标空间到地理参考坐标空间的转换

可以用GDAL读取出每个影像文件的空间地理变换，地理变换是从图像坐标空间（行、列），也称为（像素、线）到地理参考坐标空间（投影或地理坐标）的仿射变换。可以用GetGeoTransform读取出来，由以下六个参数组成

• GT(0) 左上像素左上角的x坐标。
• GT(1) w-e像素分辨率/像素宽度。
• GT(2) 行旋转（通常为零）。
• GT(3) 左上像素左上角的y坐标。
• GT(4) 列旋转（通常为零）。
• GT(5) n-s像素分辨率/像素高度（北上图像为负值）。

从图像坐标空间到地理参考坐标空间的转换关系如下：

X_geo = GT(0) + X_pixel * GT(1) + Y_line * GT(2)
Y_geo = GT(3) + X_pixel * GT(4) + Y_line * GT(5)
1
2

X_geo为实际X坐标，Y_geo为实际Y坐标，X_pixel代表栅格数据的第n列，Y_line代表栅格数据的第n行。
注意，上面的像素/线坐标是从左上角的（0.0,0.0）到右下角的（宽度单位为像素，高度单位为像素）。因此，左上角像素中心的像素/线位置为（0.5,0.5）。
如果是北向图像：
GT(2) ， GT(4) 系数为零。
GT(1) ， GT(5) 是像素大小。
GT(0) ， GT(3) 位置是栅格左上角像素的左上角。

计算公式

从前面的关系中，可以倒退出对应的计算关系，和计算二元一次方程是一样的。

X_pixel =(GT(5)*X_geo-GT(2)*Y_geo-GT(5)*GT(0)+GT(2)*GT(3))/(GT(5)*GT(1)-GT(2)*GT(4))
Y_line =(Y_geo-GT(3)-xPixel*GT(4))/GT(5)
1
2

注意：
GT(2) ， GT(4)不能作为被除数，这可能导致计算错误

代码实现

package com.meteorology.clock.common;

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 RasterTool {
    public static void main(String[] args) {
        //指定经纬度
        double Latitude  =  28.850445;
        double longitude = 119.100341;

        String filepath = "E:\\Z_NAFP_C_BABJ_20230712180520_P_HRCLDAS_RT_BEZZ_0P01_HOR-WIND-2023071218.img";
        //注册
        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));



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

        //判断是否坐标超出范围
        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){
//        GT(0) 左上像素左上角的x坐标。
//        GT(1) w-e像素分辨率/像素宽度。
//        GT(2) 行旋转（通常为零）。
//        GT(3) 左上像素左上角的y坐标。
//        GT(4) 列旋转（通常为零）。
//        GT(5) n-s像素分辨率/像素高度（北上图像为负值）

//        如果是北向图像：
//        GT(2) ， GT(4) 系数为零。
//        GT(1) ， GT(5) 是像素大小。
//        GT(0) ， GT(3) 位置是栅格左上角像素的左上角。

//        X_geo = GT(0) + X_pixel * GT(1) + Y_line * GT(2)
//        Y_geo = GT(3) + X_pixel * GT(4) + Y_line * GT(5)
        int[] ints = new int[2];
        //向下取整,如果向上取整会导致计算结果偏大，从而在后面读取到邻近像元的数据
        double xPixel=(gt[5]*X-gt[2]*Y-gt[5]*gt[0]+gt[2]*gt[3])/(gt[5]*gt[1]-gt[2]*gt[4]);
        double yLine=(Y-gt[3]-xPixel*gt[4])/gt[5];
        ints[0]= (int) Math.floor(xPixel);
        ints[1]=(int) Math.floor(yLine);
        System.out.println(ints[0]);
        System.out.println(ints[1]);
        return ints;
    }
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

结果验证

使用ArcGIS工具，打开栅格影像，并查询对应的像素值

程序结果：

从上面可以看出，结果一致，而且小数位比ArcGIS更多。

注意

我这里使用的测试栅格影像是WGS84坐标系的，所以这样的计算没有问题，要是其他坐标系的，需要进行坐标转换。