<think>由于绘制地形图通常需要地理数据,而中越边境的地形图需要该地区的海拔数据。我们可以使用Python中的`rasterio`库来读取地理栅格数据(如DEM数字高程模型),并使用`matplotlib`进行可视化。此外,我们还可以使用`contextily`添加底图。
步骤:
1. 获取中越边境的DEM数据。我们可以使用公开的DEM数据,例如SRTM(航天飞机雷达地形测绘任务)数据,该数据提供了全球范围的高程数据,分辨率为1弧秒(约30米)。
2. 使用rasterio读取DEM数据,并裁剪出中越边境的区域(大致范围:东经102°至109°,北纬21°至24°)。
3. 使用matplotlib绘制等高线图或地形阴影图。
4. 添加颜色条和标题。
5. 使用contextily添加底图(可选)。
注意:由于DEM数据文件可能较大,这里我们假设已经下载了包含中越边境的SRTM数据文件(例如,srtm_60_06.tif,具体文件根据区域而定)。如果没有,可以从NASA EarthData等网站下载。
由于无法直接下载数据,这里我将展示一个示例代码,假设数据文件为'SRTM_NE_250m.tif',并展示如何裁剪和绘制。
如果还没有数据,可以使用`geopandas`和`rasterio`从OpenStreetMap获取边界,然后使用`elevation`库下载DEM数据。但为了简化,我们假设已有数据。
以下是一个示例代码,包括如何读取DEM、裁剪、绘制地形图,并添加底图。
注意:需要安装的库有:rasterio, matplotlib, numpy, contextily, cartopy(可选,用于地理投影)
由于中越边境区域跨越多个SRTM图块,可能需要拼接,但这里我们假设已经有一个合并的DEM文件。
我们将使用rasterio裁剪出中越边境的大致范围(102°E-109°E, 21°N-24°N),然后绘制。
代码步骤:
1. 导入库
2. 打开DEM文件
3. 定义裁剪区域(经纬度范围)
4. 将经纬度范围转换为图像的行列号
5. 读取数据
6. 绘制地形图(例如使用山体阴影)
7. 添加底图
如果使用cartopy,可以更精确地处理地理坐标,但这里我们使用简单的rasterio和matplotlib,然后叠加底图。
由于我们无法确定具体的数据文件,这里将使用虚拟数据或示例数据来展示流程。在实际应用中,请替换为正确的DEM文件路径。
以下是代码:</think>### Python绘制中越边境地形图
以下是使用Python绘制中越边境地形图的完整解决方案,使用`rasterio`读取高程数据,`matplotlib`绘制地形图,`cartopy`处理地理信息:
```python
import rasterio
from rasterio.plot import show
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.colors import LinearSegmentedColormap
# 中越边境大致范围 (东经102-109°, 北纬21-24°)
BORDER_BOX = (102, 21, 109, 24) # min_lon, min_lat, max_lon, max_lat
# 创建自定义地形色带
terrain_cmap = LinearSegmentedColormap.from_list('terrain', [
'#005000', # 深绿 (低海拔)
'#008000', # 绿
'#80C080', # 浅绿
'#F0E68C', # 卡其 (丘陵)
'#D2B48C', # 棕褐 (山地)
'#A0522D', # 深棕 (高山)
'#FFFFFF' # 白 (雪峰)
])
# 创建地图画布
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
# 设置地图范围
ax.set_extent(BORDER_BOX, crs=ccrs.PlateCarree())
# 添加地理特征
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.8)
ax.add_feature(cfeature.RIVERS, edgecolor='blue', alpha=0.5)
ax.add_feature(cfeature.LAKES, edgecolor='blue', facecolor='lightblue', alpha=0.4)
# 加载高程数据 (示例使用全球DEM数据)
# 实际使用时需替换为真实DEM文件路径
try:
# 这里使用NASA SRTM数据示例
with rasterio.open('path/to/srtm_60_06.tif') as src:
# 裁剪中越边境区域
window = src.window(*BORDER_BOX)
elevation = src.read(1, window=window)
# 创建地形阴影增强立体感
x, y = np.gradient(elevation)
slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
aspect = np.arctan2(-x, y)
azimuth_rad = np.pi/4. # 光源方向(东北)
altitude_rad = np.pi/6. # 光源高度(30°)
shaded = np.sin(altitude_rad) * np.sin(slope) + \
np.cos(altitude_rad) * np.cos(slope) * \
np.cos(azimuth_rad - aspect)
shaded = (shaded * 255).astype(np.uint8)
# 获取地理坐标范围
transform = src.window_transform(window)
extent = [transform[2],
transform[2] + transform[0] * elevation.shape[1],
transform[5] + transform[4] * elevation.shape[0],
transform[5]]
# 绘制地形图
terrain = ax.imshow(
elevation,
cmap=terrain_cmap,
extent=extent,
origin='upper',
alpha=0.9,
transform=ccrs.PlateCarree()
)
# 叠加地形阴影
ax.imshow(
shaded,
cmap='gray',
extent=extent,
origin='upper',
alpha=0.3,
transform=ccrs.PlateCarree()
)
except FileNotFoundError:
print("DEM文件未找到,使用模拟数据")
# 创建模拟地形数据
x = np.linspace(BORDER_BOX[0], BORDER_BOX[2], 500)
y = np.linspace(BORDER_BOX[1], BORDER_BOX[3], 500)
xx, yy = np.meshgrid(x, y)
elevation = 2000 * np.exp(-((xx-105.5)/1.5)**2 - ((yy-22.5)/0.8)**2) + \
1500 * np.exp(-((xx-107)/1.2)**2 - ((yy-23.2)/0.6)**2)
# 绘制模拟地形
terrain = ax.contourf(
xx, yy, elevation,
levels=20,
cmap=terrain_cmap,
alpha=0.9,
transform=ccrs.PlateCarree()
)
# 添加颜色条
cbar = plt.colorbar(terrain, ax=ax, orientation='vertical', pad=0.05, aspect=30)
cbar.set_label('海拔高度 (米)', fontsize=12)
# 添加网格和标注
gl = ax.gridlines(draw_labels=True, linestyle='--', alpha=0.7)
gl.top_labels, gl.right_labels = False, False
gl.xlabel_style, gl.ylabel_style = {'size': 10}, {'size': 10}
# 添加标题和标注
plt.title('中越边境地形图', fontsize=16, pad=20)
ax.text(105.5, 23.8, '中国', fontsize=12, ha='center',
bbox=dict(facecolor='white', alpha=0.7))
ax.text(106.5, 21.5, '越南', fontsize=12, ha='center',
bbox=dict(facecolor='white', alpha=0.7))
# 添加比例尺和指北针
ax.add_patch(plt.Rectangle((102.2, 21.2), 1, 0.1, facecolor='black'))
ax.text(102.7, 21.15, '100 km', ha='center', fontsize=9)
ax.arrow(108.5, 21.5, 0, 0.3, head_width=0.1, head_length=0.1, fc='k')
ax.text(108.6, 21.5, 'N', fontsize=12)
# 保存并显示图像
plt.savefig('中越边境地形图.jpg', dpi=300, bbox_inches='tight')
plt.show()
```
### 代码解释:
1. **数据准备**:
- 使用`rasterio`读取数字高程模型(DEM)数据
- 定义中越边境区域(东经102-109°, 北纬21-24°)
- 备选方案:当DEM文件不可用时生成模拟地形数据
2. **地形可视化技术**:
- 创建自定义地形色带`terrain_cmap`,从低海拔到高海拔渐变
- 使用`np.gradient`计算地形阴影增强立体感
- 叠加地形数据和阴影层增强三维效果
3. **地理信息处理**:
- 使用`cartopy`添加海岸线、国界、河流等地理要素
- 添加经纬度网格和标签
- 使用PlateCarree投影(等距圆柱投影)
4. **地图标注**:
- 添加国家标注和比例尺
- 包含指北针方向标识
- 添加颜色条显示海拔高度
### 数据获取建议:
1. **推荐DEM数据源**:
- NASA SRTM(航天飞机雷达地形测绘任务):30米分辨率
- ASTER GDEM:30米分辨率
- Copernicus DEM:30米分辨率
2. **免费获取途径**:
- USGS EarthExplorer(搜索SRTM)
- OpenTopography.org
- NASA Earthdata Search
3. **替代方案**:
- 使用`elevation`库下载DEM数据:
```python
import elevation
elevation.clip(bounds=BORDER_BOX, output='border_dem.tif')
```
### 注意事项:
1. 实际使用时需替换DEM文件路径
2. 高分辨率DEM文件可能较大(100MB+),处理时需足够内存
3. 可调整色带颜色和透明度优化视觉效果
4. 添加`alpha`参数使多图层叠加更自然