1. 简介
获取特定经纬度处的高程数据需要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米。