在Java GDAL环境中提取指定经纬度的高程数据

1. 简介

获取特定经纬度处的高程数据需要3个步骤：

1. 根据经纬度确定数据文件
2. 确定该经纬度点所对于的栅格数据的行号与列号
3. 读取相应数据文件中对应行与列处的高程值

对于第1步，可以采用“”博文中描述的方法来获得。对于第2步，对于位置精度要求不高的场合，也可以采用“如何获取特定经纬度在SRTM3中的高程值”博文提供的简单方式获得。最后的第3步，可以采用“在Java GDAL环境中读取Tiff文件中的DEM数据”博文中介绍的方法来获取。

然而，Tiff文件中的数据通常不会如“如何获取特定经纬度在SRTM3中的高程值”博文中所描述得那样理想，能够在栅格的左上角处取得整数经纬度。在实际的数据中，栅格点所对应的经纬度是需要利用仿射变换参数，经过仿射变换后得到，如“在Java GDAL环境中读取Tiff文件中的DEM数据”博文所述。

因此，在Java GDAL环境中，利用仿射变换参数，可以获得指定经纬度点的更为准确的高程数据。

2. 经仿射变换获得栅格数据中的行与列

在Java GDAL环境中，可以利用poDataset.GetGeoTransform(gt)命令获得栅格数据的仿真变换参数gt，其中在“上北下南”的图像中，系数gt(2)与gt(4) 为0， 系数gt(1) 是像素宽度，gt(5)是像素高度，(gt(0),gt(3))位置是栅格左上角像素的经纬度。

栅格坐标和地理坐标之间的仿射变换关系为

表示实际的地址坐标，即经纬度；为栅格数据中的坐标。

根据给定点的经纬度来获得栅格数据的行与列，则需要对上面的仿射变换进行逆运算，即

计算出的进行四舍五入为整数，即为相应的行与列。

3. 示例

给定点的经纬度为（E116.75, N37.112），由SRTM3数据中提取出相应的高程信息。

首先需要由经纬度（E116.75, N37.112）确定该点所在文件的名称为"srtm_60_05.tif"，具体方法详见“如何获取特定经纬度在SRTM3中的高程值”，本文不再赘述。

生成文件名的代码如下：

/**
 * 生成*.tif文件名称
 * @param longitude：文件描述区域中心点的经度，单位：度
 * @param latitude：文件描述区域中心点的纬度，单位：度
 * --------------------------------------------
 * @return fileName：文件名称，如："srtm_" + XX + "_" + YY +".tif"
 */
public static String GenerateFileName(double longitude, double latitude) {
	String fileName = "";
	String XX, YY;
	
	int longX = (int) Math.round((longitude + 180)/5+0.5);
	int latY = (int) Math.round((60 - latitude)/5+0.5);
	if (longX < 10) {
		XX = "0" + String.valueOf(longX);
	} else {
		XX = String.valueOf(longX);
	}
	if (latY < 10) {
		YY = "0" + String.valueOf(latY);
	} else {
		YY = String.valueOf(latY);
	}
	
	fileName = "srtm_" + XX + "_" + YY +".tif";
	
	return fileName;
}

根据文件名srtm_60_05.tif，打开相关数据文件，提取出对应的仿射变换参数。由文件“”获得的仿射变换参数为：

• gt(0) = 114.99958393868292,
• gt(1) = 8.333333333333334E-4,
• gt(2) = 0.0,
• gt(3) = 40.000416763514124,
• gt(4) = 0.0,
• gt(5) = -8.333333333333334E-4

因此，可以利用上述仿射变换参数，根据逆变换公式计算出栅格数据文件中所对应的行与列，相关代码为

/**
 * 将地理坐标转换为栅格像素坐标
 * @param gt 仿射变换参数
 * @param J 经度，单位：度
 * @param W 纬度，单位：度
 * @return
 */
public static int[] JW2ColRow(double[] gt, double J, double W){
    int[] colrow = new int[2];

    // 四舍五入，为了选择邻近像元的数据
    double  Xp = Math.round(((J - gt[0])*gt[5] - (W-gt[3])*gt[2]) / (gt[1]*gt[5]-gt[2]*gt[4]));
    double  Yp = Math.round(((W - gt[3])*gt[1] - (J-gt[0])*gt[4]) / (gt[1]*gt[5]-gt[2]*gt[4]));

    colrow[0] = new Double(Xp).intValue();
    colrow[1] = new Double(Yp).intValue();

    return colrow;
}

由上述代码计算出的行与列分别为3466和2100。据此可以读出其对应的高程为22米。