目录
一、前言
在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(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);
}
}
/**
* 将地图坐标转换为栅格像素坐标
* @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?