提取特定经纬度范围内的高程数据

最新推荐文章于 2024-07-24 11:49:51 发布

带着地球去浪一浪

最新推荐文章于 2024-07-24 11:49:51 发布

阅读量2k

点赞数 2

分类专栏：电磁波覆盖范围分析文章标签：数据结构算法

本文链接：https://blog.csdn.net/weixin_42428226/article/details/129519237

版权

电磁波覆盖范围分析专栏收录该内容

17 篇文章 4 订阅

订阅专栏

1. 任务

以特定点，如：西藏林芝米林机场附近的一个山上（E94.36488，N29.29815），为中心，获取经度跨度为 $7.5^{\circ}$ ，纬度跨度为 $7.5^{\circ}$ 范围的高程数据，并将其存贮在 $700\times 700$ 的数组中。

2. 分析

根据任务内容可获得如下信息：

经纬度区间为 $7.5^{\circ}\times 7.5^{\circ}$ ，需要多个DEM文件来支撑相关区间的数据
数组为 $1000\times 1000$ ，数据需要由数据文件中抽取
为提高数据获取效率，需要尽量减少文件打开的次数

3. 相关参数计算

读取指定区域的高程数据的关键在于确定相应的数据文件及该文件所填充数组的行与列。

不妨设置如下参数：

$X_c$ ：中心点的经度
$Y_c$ ：中心点的纬度
$w_x$ ：显示区域的经度跨度
$w_y$ ：显示区域的纬度跨度
$N_x$ ：数组列数
$N_y$ ：数组行数
$N_{xb}$ ：特定数据文件中的有效数据在数组中的起始列号
$N_{yb}$ ：：特定数据文件中的有效数据在数组中的起始行号
$N_{xe}$ ：特定数据文件中的有效数据在数组中的终止列号
$N_{ye}$ ：特定数据文件中的有效数据在数组中的终止行号
$X_b$ ：数组中数据在特定数据文件中最左上角处数据的经度
$Y_b$ ：数组中数据在特定数据文件中最左上角处数据的纬度
$X_0$ ：数据起始点的经度
$Y_0$ ：数据起始点的纬度

对于数组中数据在特定数据文件中最左上角处数据的经纬度 $(X_b,Y_b)$ 可以利用中心点的经纬度 $(X_c,Y_c)$ 、经纬跨度 $w_x$ 和 $w_y$ 、数组行列数 $N_x$ 和 $N_y$ 等参数计算得出。数组左上角点的经纬度可方便地计算出

$\left\{\begin{matrix} X_0 = X_c - \frac{1}{2}w_x\\ Y_0 = Y_c - \frac{1}{2}w_y \end{matrix}\right.$

根据 $(X_0,Y_0)$ 可以参考如何获取特定经纬度在SRTM3中的高程值文中介绍的方法确定SRTM3文件名称SRTM_XX_YY.zip。

根据数组中行与列的经纬度间隔 $w_x/N_x$ 和 $w_y/N_y$ 确定数组中每个元素的经纬度，并依据经纬度将存在于该数据文件中数据全部读出，赋值给数组相应元素。记录填充元素的行列号。

由于显示区域是连续的，因此，当处理完成一个SRTM3文件后，只需要将其附近的相关文件打开处理即可。

遍历显示区域相关的所有文件，完成对数组数据的填充。

4. 数据提取步骤

1）获取显示区域的经纬度范围

2）获取显示区域所涉及的数据文件的经纬度编号

3）对于给定文件，首先需要确定该文件所对应的起止经纬度，以及该文件可能填充的数组的起始元素编号

4）根据与文件相关的XX与YY值，生成文件名

5）打开相应文件，将数据读入数组中

5. 实现

根据上述内容，将部分重点代码提供如下：

1）获取显示区域的经纬度范围

因为没有考虑Lambert变换，只是经纬度的直接获取，所以可以直接利用系统中心与经纬度跨度直接计算得出。

Range_LongLat range_LongLat = new Range_LongLat();
range_LongLat.min_LongitudeInImage = lon_SystemCenter - lon_range/2;
range_LongLat.max_LongitudeInImage = lon_SystemCenter + lon_range/2;
range_LongLat.min_LatitudeInImage = lat_SystemCenter - lat_range/2;
range_LongLat.max_LatitudeInImage = lat_SystemCenter + lat_range/2;

在这里定义了一个新的类Range_LongLat，用于记录区域内的最大最小经纬度值。

class Range_LongLat {
	public double max_LongitudeInImage;	//显示区域中最大的经度值
	public double min_LongitudeInImage;	//显示区域中最小的经度值
	public double max_LatitudeInImage;	//显示区域中最大的纬度值
	public double min_LatitudeInImage;	//显示区域中最小的纬度值
	
}

2）获取显示区域所涉及的数据文件的经纬度编号

/**
 * 获取显示区域所涉及的数据文件的经纬度编号
 * int no_LongitudeBegin; 	//地图数据文件中经度的起始编号
 * int no_LongitudeEnd; 	//地图数据文件中经度的终止编号
 * int no_LatitudeBegin; 	//地图数据文件中纬度的起始编号
 * int no_LatitudeEnd; 		//地图数据文件中纬度的终止编号
 */
public static No_TiffName getTiffNo(Range_LongLat myrange_LongLat) {
	No_TiffName myno_TiffName = new No_TiffName();
	myno_TiffName.no_LongitudeBegin = (int) ((myrange_LongLat.min_LongitudeInImage + 180)/5.0 + 1);
	myno_TiffName.no_LongitudeEnd   = (int) ((myrange_LongLat.max_LongitudeInImage + 180)/5.0 + 1);
	myno_TiffName.no_LatitudeBegin = (int) ((60-myrange_LongLat.max_LatitudeInImage)/5.0 + 1);
	myno_TiffName.no_LatitudeEnd   = (int) ((60-myrange_LongLat.min_LatitudeInImage)/5.0 + 1);
		
	return myno_TiffName;
}

在此，引用了一个自定义的类No_TiffName，用于说明SRTM文件名是的XX与YY的起始与终止值。

class No_TiffName {
	public int no_LongitudeBegin; 	//地图数据文件中经度的起始编号
	public int no_LongitudeEnd; 	//地图数据文件中经度的终止编号
	public int no_LatitudeBegin; 	//地图数据文件中纬度的起始编号
	public int no_LatitudeEnd; 		//地图数据文件中纬度的终止编号
	
}

根据no_TiffName变量的值，遍历符合条件的SRTM文件。

3）对于给定文件，首先需要确定该文件所对应的起止经纬度，以及该文件可能填充的数组的起始元素编号

//确定tiff文件的经纬度起止值，单位：度
double lat_High_Begin = 60-5*(no_Latitude-1);
double lat_Low_End = lat_High_Begin - 5;
double lon_Low_Begin = 5*(no_Longitude-1)-180;
double lon_High_End = lon_Low_Begin + 5;
		        
//确定tiff文件在显示区域内的位置，起止像素点编号
Point2D.Double leftUp = new Point2D.Double();
leftUp.x = (lon_Low_Begin - lon0)/delta_x;	// 文件中的最小经度相对起始经度间隔多少单元
leftUp.y = (lat0 - lat_High_Begin)/delta_y;	// 文件中的最大纬度相对起始纬度间隔多少单元
Point2D.Double rightDown = new Point2D.Double();
rightDown.x = (lon_High_End - lon0)/delta_x;	// 文件中的最小经度相对起始经度间隔多少单元
rightDown.y = (lat0 - lat_Low_End)/delta_y;	// 文件中的最大纬度相对起始纬度间隔多少单元
int wn_Begin = Math.max((int)leftUp.x, 0);
int wn_End = Math.min((int)rightDown.x, width_Image-1);
int hn_Begin = Math.max((int)leftUp.y, 0);
int hn_End = Math.min((int)rightDown.y, height_Image-1);

4）根据与文件相关的XX与YY值，生成文件名

/**
 * 生成*.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;
}

5）打开相应文件，将数据读入数组中

//将tiff文件中的高程添加到数据矩阵imageDEMData[wn][hn]中
int values[] = new int[1];
for (int wn = wn_Begin; wn <= wn_End; wn++){
	for (int hn = hn_Begin; hn <= hn_End; hn++){
		Point2D.Double lon_Lat = new Point2D.Double();
		lon_Lat.x = lon0 + wn * delta_x;
		lon_Lat.y = lat0 - hn * delta_y;
		if ((lon_Lat.x >= lon_Low_Begin)&(lon_Lat.x < lon_High_End)&
                    (lon_Lat.y > lat_Low_End)&(lon_Lat.y <= lat_High_Begin)) {
			int x = (int) ((lon_Lat.x - lon_Low_Begin)/5.0*6000.0);
			int y = (int) ((lat_High_Begin - lon_Lat.y)/5.0*6000.0);
			poBand.ReadRaster(x, y, 1, 1, values);
			imageDEMData[wn][hn] = Math.max(values[0], 0);
		}
	}
}

根据本文任务提取的高程数据所形成的地形图如下