<think>我们正在处理DEM高程数据的拼接问题。根据引用[1]和[3]，我们可以知道DEM数据通常以栅格文件（如.tif）存储，并且拼接过程包括查找文件、读取文件、合并文件等步骤。引用[3]提供了使用Python查找符合条件的DEM文件并写入列表的代码。引用[1]提到ASTGTM数据是按经纬度分块的，因此拼接需要将这些分块合并成一个大的栅格。
接边处理：在拼接相邻的DEM数据时，边缘部分可能会出现不连续的情况，因此需要进行接边处理，使相邻图块之间的过渡自然。常见的接边处理方法包括平滑过渡（如使用加权平均）或羽化（feathering）技术。
根据引用[2]，我们可以使用数据平滑的方法来处理接边，比如移动平均法、高斯滤波等。在接缝处，我们可以对重叠区域进行加权平均，权重根据距离接缝的远近而变化。
步骤：
1. 读取所有需要拼接的DEM文件（通常为GeoTIFF格式）。
2. 确定每个文件的地理位置（根据文件名或文件中的地理信息）。
3. 按照地理位置排列这些文件（例如，按经度、纬度排序）。
4. 确定最终拼接后的大图范围，并创建一个空的大矩阵（或使用栅格处理库创建新栅格）。
5. 将每个图块的数据放入大矩阵的相应位置。
6. 对于相邻图块的重叠区域，进行接边处理（如加权平均）。
但是，在实际操作中，我们通常使用现成的库来处理栅格数据的拼接，因为手动处理地理坐标变换和重采样等非常复杂。常用的库有rasterio（基于GDAL）和gdal。
下面是一个使用rasterio和gdal进行DEM拼接的示例步骤（注意：这里我们以rasterio为例，因为其API更友好，但gdal也有拼接功能）。
然而，需要注意的是，引用[4]提到ArcGIS和Erdas的拼接工具，但我们现在使用Python。
使用rasterio.merge模块的merge函数可以方便地拼接多个栅格。该函数会自动处理重叠区域（默认取最大值，但我们可以通过方法参数指定，比如取平均值）。但是，对于接边处理，可能需要自定义一个函数来对重叠区域进行平滑。
以下是一个使用rasterio进行拼接的基本流程：
步骤1：安装必要的库（rasterio和numpy）
步骤2：读取所有DEM文件
步骤3：使用rasterio.merge.merge函数进行拼接
但是，rasterio.merge.merge默认的拼接方法有：
- 第一个像素（覆盖）
- 最后一个像素（覆盖）
- 最小值、最大值、平均值等
如果我们想要进行接边平滑，可能需要分两步：
a. 先拼接，不处理重叠区域（比如使用第一个像素，即后放入的覆盖先放入的）
b. 然后找到所有重叠区域（可以通过记录每个图块放置的位置，然后检查哪些像素位置有多个图块覆盖）
c. 对重叠区域进行加权平均（例如，根据每个像素到各自图块边界的距离设置权重）
然而，手动实现接边处理比较复杂。另一种思路是：在拼接时，每个图块在重叠区域按照一定的权重进行融合。我们可以使用rasterio.merge.merge函数的自定义回调函数（通过`method`参数传递一个函数）来实现。
但是，rasterio.merge.merge函数目前内置的方法只有几种，并不支持自定义的加权平均。因此，我们可以考虑以下方法：
方法1：先使用rasterio.merge.merge拼接，然后对重叠区域进行后处理（需要记录每个位置有哪些图块贡献，然后重新计算平均值，但这样需要保存每个图块在最终图像中的贡献次数和值，比较复杂）。
方法2：分块处理，每次处理两个图块的拼接，并逐步合并。在每次合并两个图块时，对重叠区域进行加权平均。这样，我们可以控制每个接缝的处理。
由于方法2更可控，但效率较低（需要多次合并），而方法1需要额外的内存记录贡献次数，我们选择一种折中的方式：使用rasterio.merge.merge先拼接，然后识别重叠区域（通过非0像素的叠加区域）并应用平滑滤波（如高斯滤波）到整个图像的接缝处？但这样可能会影响非接缝区域。
因此，我们可能需要一个更直接的方法：在合并时，对于重叠区域，我们计算每个像素来自多个图块，然后取平均值（或加权平均）。这可以通过以下步骤：
1. 创建一个空的目标数组和一个计数数组（记录每个像素被覆盖的次数）。
2. 遍历每个图块，将图块的值加到目标数组的对应位置，同时计数数组对应位置加1。
3. 最后，将目标数组除以计数数组（即取平均值）。
但是，这种方法没有考虑图块边缘的权重（即距离图块中心越近的像素权重越高，而边缘则低）。我们可以改进为：
1. 为每个图块创建一个权重数组（例如，图块中心权重为1，边缘为0，中间渐变）。可以使用距离图块边界的距离来生成权重（如线性权重或高斯权重）。
2. 然后，创建一个空的目标数组和一个权重累积数组（全零，与目标数组同形状）。
3. 对于每个图块，将图块数据乘以权重数组后加到目标数组，同时将权重数组加到权重累积数组。
4. 最后，目标数组除以权重累积数组（避免除零错误）。
这种方法称为羽化（feathering）或混合（blending）。但是，它需要每个图块都有权重数组，并且要知道每个图块在目标图像中的位置。
由于DEM图块通常是规则的网格（相邻图块有相同的行数或列数），我们可以假设每个图块的大小相同，并且相邻图块有重叠（但实际DEM数据可能是无缝相邻的，即没有重叠，但边缘是连续的？）。实际上，ASTER GDEM数据是分块提供的，相邻块之间有一列重叠（具体需要看数据说明）。如果没有重叠，我们就需要创建重叠区域（通过扩展图块边缘），然后进行混合。
考虑到复杂性，我们首先尝试使用简单的平均值方法（即上述的计数平均），然后再考虑更复杂的权重方法。
实现步骤（使用rasterio和numpy）：
1. 读取所有DEM文件，并获取它们的地理范围（边界）。
2. 计算整个拼接区域的范围（最小最大经纬度）。
3. 根据整个范围和分辨率计算目标图像的尺寸（行数和列数）。
4. 创建目标数组（全零或NaN）和权重累积数组（全零）。
5. 对于每个DEM文件：
a. 读取数据数组和变换矩阵。
b. 计算该图块在目标图像中的位置（起始行、起始列，结束行、结束列）。
c. 创建一个权重数组（与图块数组同形状）。例如，我们可以使用一个简单的线性权重：从图块边缘向内，权重从0线性增加到1（或者使用常数1，这样就变成简单平均）。
d. 将图块数据数组乘以权重数组，然后加到目标数组的对应位置。
e. 将权重数组加到权重累积数组的对应位置。
6. 将目标数组除以权重累积数组（注意避免除0，将0的位置设为NaN或0）。
7. 写入输出文件。
注意：如果图块之间没有重叠，那么权重累积数组在非重叠区域为1，重叠区域为重叠的图块数（如果使用常数权重）。如果使用边缘权重，则重叠区域的权重是多个图块权重的叠加。
下面是一个简化的代码框架（使用常数权重1，即简单平均）：
但是，简单平均在接缝处可能会产生模糊，特别是当相邻图块有较大差异时。因此，我们更倾向于使用边缘权重（羽化），使得每个图块在边缘的贡献小，而中心贡献大。
如何生成权重数组？
我们可以为每个图块生成一个权重数组，该数组在中心为1，在边缘（比如1像素宽）为0.5，然后渐变？或者使用距离函数：
权重 = 1 - 距离边缘的归一化距离（在边缘处为0，中心为1） 或者使用其他函数。
例如，我们可以这样计算权重数组（以图块宽度和高度为基准）：
权重 = (distance_from_left_boundary / width) * (distance_from_right_boundary / width) * ... 这样并不好。
一个更简单的方法是：分别计算每个像素到上、下、左、右边界的距离，然后取最小值，再除以到中心的最大距离（即宽度/2或高度/2）？这样中心为1，边缘为0。
但是，我们也可以使用二维的线性渐变。例如，对于图块的左边缘，从左到右权重从0线性增加到1（在左边界处权重0，右边界处权重1）？这样不行，因为右边缘可能和下一个图块重叠。所以，我们可以在两个方向（上下左右）都做渐变，然后取最小值？或者分别对左右和上下做渐变，然后相乘？
实际上，我们可以生成一个权重数组，使得图块的四个边缘的权重为0，然后向中心线性增加到1。具体：
1. 创建一个与图块相同形状的权重数组，初始值全1。
2. 对于每个边缘（左、右、上、下），分别创建一个线性渐变权重（从0到1），然后取这四个权重的乘积（或最小值）作为最终的权重数组。
例如，对于左边缘，我们创建一个从左到右的渐变：从0（最左）到1（最右），但注意，我们只希望边缘部分渐变，比如渐变范围是边缘的N个像素（如10个像素）。那么：
left_weight = np.linspace(0, 1, N) # 前N个像素渐变
然后，将left_weight扩展到整个高度（重复N次），形成一个二维数组（宽度为N，高度为图块高度），然后图块左侧N列用这个数组，右侧N列用类似的从1到0的渐变（右边缘），上下同理。
然后，最终的权重数组是这四个权重的乘积（这样四个角就会是0，而中心是1）。
但是，这样处理会使四个角为0，中心为1，而边缘部分渐变。在重叠区域，多个图块的权重数组会叠加，然后我们按照加权平均进行融合。
考虑到复杂度，我们可以先尝试在重叠区域使用简单的线性渐变（只考虑左右或上下，因为相邻图块通常是左右或上下相邻，不会对角相邻）。所以，我们可以只考虑左右边缘（对于左右相邻的图块）和上下边缘（对于上下相邻的图块）。但是，在多个图块重叠的角落（四个图块重叠），我们需要四个方向都考虑。
因此，我们为每个图块生成一个权重数组，该数组在四个边缘都进行渐变（渐变宽度为overlap_width，比如10个像素）。渐变方式：从边缘向内，权重从0线性增加到1（在边缘处为0，在向内overlap_width处为1，然后保持1直到另一边缘）。
具体步骤：
1. 创建一个全1的权重数组（形状为图块的高度和宽度）。
2. 左边缘渐变：从左边界开始，向右overlap_width个像素，权重从0线性增加到1。
weight_array[:, 0:overlap_width] = np.linspace(0, 1, overlap_width)[np.newaxis, :]
3. 右边缘渐变：从右边界向左overlap_width个像素，权重从1线性减小到0。
weight_array[:, -overlap_width:] = np.linspace(1, 0, overlap_width)[np.newaxis, :]
4. 上边缘渐变：从上边界向下overlap_width个像素，权重从0线性增加到1。
weight_array[0:overlap_width, :] = np.linspace(0, 1, overlap_width)[:, np.newaxis]
5. 下边缘渐变：从下边界向上overlap_width个像素，权重从1线性减小到0。
weight_array[-overlap_width:, :] = np.linspace(1, 0, overlap_width)[:, np.newaxis]
6. 然后，对于四个角，我们取对应边缘权重的最大值或最小值？这里我们取乘积，但这样四个角会变成0（因为两个0相乘），所以我们可以取每个位置的最大值（这样角落的权重就是两个边缘权重的最大值，即角落的权重是0，因为两个边缘都是0）。或者，我们可以取每个位置的最小值？这样也不行，因为角落同时受到两个边缘的影响，我们希望角落的权重是0。所以，乘积也可以，因为只要有一个边缘权重为0，乘积就是0。
但是，这样生成的权重数组在四个角为0，在中心为1，在边缘渐变区域为渐变值。
然后，我们按照这个权重数组进行加权平均拼接。
注意：重叠区域（多个图块重叠）的权重累积数组的值会大于1（因为多个图块都有权重），然后我们进行归一化（加权平均）。
代码实现步骤：
1. 定义重叠像素宽度（overlap_width，比如10）。
2. 计算整个拼接区域的范围（最小经度、最大经度、最小纬度、最大纬度）和分辨率（假设所有图块分辨率相同）。
3. 根据总范围和分辨率计算目标图像的行数和列数。
4. 创建目标数组（全0或NaN）和权重累积数组（全0），形状为（行数，列数）。
5. 遍历每个图块：
a. 读取图块数据（数组arr）和变换矩阵（transform）。
b. 根据图块的变换矩阵和尺寸计算其在目标图像中的位置（行起始、行结束，列起始、列结束）。
c. 创建该图块的权重数组（如上所述，形状与arr相同）。
d. 将arr乘以权重数组，然后加到目标数组的对应位置。
e. 将权重数组加到权重累积数组的对应位置。
6. 将目标数组除以权重累积数组（注意：分母为0的地方设为NaN或0）。
7. 创建目标图像的变换矩阵（左上角坐标和分辨率）。
8. 写入输出GeoTIFF文件。
注意：读取DEM文件时，我们使用rasterio。另外，需要处理NoData值（通常DEM中会用某个特定值表示无效，如-9999）。在计算加权平均时，应该忽略NoData值。因此，我们需要在读取时设置mask，并在计算时跳过NoData。
改进：在步骤5中，对于每个像素，如果图块中的像素是NoData，则权重数组对应位置设为0（这样就不会参与计算），并且该像素不参与累加。
因此，我们在读取图块时，获取其NoData值（如果没有，则默认为None），然后创建一个掩码（mask）数组（有效数据为True，NoData为False）。然后，在权重数组中，将NoData位置设为0。
步骤修改：
a. 读取图块数据（arr）和NoData值（nodata），以及transform。
b. 创建掩码数组：mask = (arr != nodata) 注意：如果nodata为None，则跳过（即全部有效）。
c. 创建权重数组（全1，然后设置边缘渐变），然后将权重数组乘以mask（这样NoData位置的权重为0）。
d. 然后继续。
这样，NoData值就不会参与加权平均。
下面我们给出代码示例（注意：由于代码较长，我们只写关键部分，并假设所有图块的分辨率相同，且坐标系相同）。
代码示例（使用rasterio）：</think>### Python实现DEM高程数据接边处理与拼接
在Python中实现DEM高程数据的接边处理和拼接，主要使用`rasterio`和`gdal`库处理栅格数据，结合`numpy`进行数值计算。以下是关键步骤和代码实现：
---
### **1. 数据准备与文件查找**
```python
import os
import rasterio
from rasterio.merge import merge
from rasterio.warp import calculate_default_transform, reproject, Resampling
import numpy as np
# 查找所有DEM文件（假设为.tif格式）
def find_dem_files(directory, suffix='.tif'):
dem_files = []
for root, _, files in os.walk(directory):
for file in files:
if file.endswith(suffix) and 'dem' in file.lower():
dem_files.append(os.path.join(root, file))
return dem_files
dem_dir = "path/to/DEMdata"
dem_files = find_dem_files(dem_dir)
```
---
### **2. 接边处理（羽化混合）**
相邻DEM图块边缘需通过**加权平均**消除接缝，权重由像素到边界的距离决定：
```python
def feather_edges(arr, feather_width=50):
"""对数组边缘进行羽化处理"""
height, width = arr.shape
weights = np.ones_like(arr, dtype=float)
# 生成左右边缘权重（线性渐变）
for i in range(feather_width):
weights[:, i] = i / feather_width
weights[:, width - 1 - i] = i / feather_width
# 生成上下边缘权重
for i in range(feather_width):
weights[i, :] = np.minimum(weights[i, :], i / feather_width)
weights[height - 1 - i, :] = np.minimum(weights[height - 1 - i, :], i / feather_width)
return arr * weights
# 示例：处理单个DEM边缘
with rasterio.open(dem_files[0]) as src:
dem_data = src.read(1)
processed_data = feather_edges(dem_data)
```
---
### **3. 多图块拼接**
使用`rasterio.merge`实现无缝拼接：
```python
def merge_dems(dem_files, output_path):
"""拼接DEM并保存"""
src_files = [rasterio.open(f) for f in dem_files]
mosaic, transform = merge(src_files, method='first') # 优先取非空值
# 获取元数据
meta = src_files[0].meta.copy()
meta.update({
"height": mosaic.shape[1],
"width": mosaic.shape[2],
"transform": transform
})
# 写入结果
with rasterio.open(output_path, 'w', **meta) as dest:
dest.write(mosaic)
# 关闭所有文件
for src in src_files:
src.close()
output_path = "merged_dem.tif"
merge_dems(dem_files, output_path)
```
---
### **4. 接缝优化（后处理）**
拼接后对重叠区进行**高斯滤波**平滑：
```python
from scipy.ndimage import gaussian_filter
def smooth_seams(input_path, output_path, sigma=2):
"""高斯平滑接缝区域"""
with rasterio.open(input_path) as src:
data = src.read(1)
meta = src.meta
# 识别重叠区（值为0表示无数据）
overlap_mask = (data == 0)
# 高斯滤波（仅作用于非空区域）
smoothed = gaussian_filter(data, sigma=sigma)
data[~overlap_mask] = smoothed[~overlap_mask]
# 保存结果
with rasterio.open(output_path, 'w', **meta) as dest:
dest.write(data, 1)
smooth_seams("merged_dem.tif", "smoothed_dem.tif")
```
---
### **关键注意事项**
1. **坐标系统一**
确保所有DEM使用相同的坐标系（如WGS84），否则需重投影：
```python
def reproject_dem(input_path, output_path, target_crs='EPSG:4326'):
with rasterio.open(input_path) as src:
transform, width, height = calculate_default_transform(
src.crs, target_crs, src.width, src.height, *src.bounds)
meta = src.meta.copy()
meta.update({
'crs': target_crs,
'transform': transform,
'width': width,
'height': height
})
with rasterio.open(output_path, 'w', **meta) as dest:
reproject(
source=rasterio.band(src, 1),
destination=rasterio.band(dest, 1),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=transform,
dst_crs=target_crs,
resampling=Resampling.bilinear)
```
2. **处理NoData值**
拼接时忽略无效值（如-9999）：
```python
mosaic, transform = merge(src_files, nodata=-9999)
```
3. **性能优化**
- 分块处理：使用`rasterio.windows.Window`分块读写大文件
- 并行计算：用`concurrent.futures`多线程处理图块
---
### **应用场景**
- 地形分析（坡度、流域提取）
- 三维可视化
- 遥感影像正射校正
- 洪水模拟[^1]
本文介绍如何使用GDAL库从DEM数据中读取特定经纬度的高程信息，包括库的使用方法、代码实现及结果验证。
扫一扫
601
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
