当前位置：首页>python>用Python绘制博士论文研究区概况图:以西北干旱区DEM为例

用Python绘制博士论文研究区概况图:以西北干旱区DEM为例

2026-08-20 15:36:11

讲师招募 | 免费数据资源 |最新最热

直播时间：8月14日-17日

基于AI Agent（Codex · Claude Code · Hermes）的文献计量学+Meta分析一体化融合——选题论证、证据合成、成果交付及可迁移、可复制的自动化工作流实践技术培训班

直播时间：8月21日-24日

基于Claude Code 、Codex双AI协同高水平论文撰写与质量校准：研究定位→数据分析→论文初稿→交叉审稿→投稿与返修全流程实践培训班

直播时间：8月29日-30日

基于Claude Code与Codex双AI Agent协作的WebGIS项目全链路开发与生产级部署实战高级培训班

直播时间：9月5日-6日

一、原文勘误：几个需要注意的细节

在参考同类文章时，我发现以下几个容易踩坑的地方，先帮你排雷：

问题	原文情况	修正说明
经纬度范围	73°E–108°E, 34°N–50°N	不同文献界定略有差异。多数研究采用 73°E–107°E, 35°N–50°N，本文取较宽范围以覆盖完整研究区。
ETOPO2 变量名	直接使用 `x/y/z`	ETOPO2 的 NetCDF 变量名因下载版本而异，可能是 `lon/lat/elevation` 或 `x/y/z`。建议用 `xarray` 自动识别，避免硬编码。
m_map 工具箱	未说明需额外安装	`m_map` 是 Matlab 第三方工具箱，非内置函数。Python 的 `cartopy` 同样需 `conda/pip` 单独安装。
色带文件	依赖外部 `colorsave.mat`	该文件为自定义色带，原文未提供下载方式。本文改用 Matplotlib 原生 `terrain` + 自定义 `ListedColormap`，零外部依赖。
代码索引	lat(indices_lat,:)	`lat` 为一维向量，不应使用二维索引 `(:,)`，运行会报错。

二、研究区概况：中国西北干旱区

中国西北干旱区位于欧亚大陆腹地，深居内陆，远离海洋，是亚洲中部干旱区的重要组成部分。

地理参数	描述
经纬度范围	约 73°E–108°E，34°N–50°N
行政区划	新疆全境、甘肃河西走廊、内蒙古阿拉善盟、宁夏北部
地貌格局	"三山夹两盆"——阿尔泰山、天山、昆仑山；准噶尔盆地、塔里木盆地
气候特征	温带大陆性干旱气候，年均降水不足 150 mm
核心矛盾	水资源短缺是制约社会经济发展的首要自然因素

该区域地形起伏剧烈（0–6000 m），地貌类型多样，是展示 DEM 制图技术的理想案例。

三、技术路线总览

相比 Matlab 的 m_map，Python 生态的 Cartopy 具有以下优势：

开源免费，无需商业授权
与 Matplotlib 深度集成，绘图语法一致
内置多种投影（PlateCarree、Lambert、Albers 等）
自动下载全球海岸线/国界/河流（Natural Earth 数据集）
与 xarray/netCDF4 无缝衔接，适合处理大规模栅格数据

四、环境配置

# 创建独立环境（推荐）conda create -n geoplot python=3.11conda activate geoplot# 核心依赖conda install-c conda-forge cartopy matplotlib numpy xarray netCDF4 shapely# 可选：用于读取 Shapefile 的引擎pip install pyshp geopandas

⚠️ 注意：Cartopy 依赖 PROJ、GEOS、Shapely 等 C 库，强烈建议通过 conda-forge 安装，pip 直接安装极易因库版本冲突导致投影变换失败。

五、数据准备

5.1 DEM 高程数据

推荐数据源：

数据集	分辨率	格式	下载地址
GEBCO	15 弧秒	NetCDF/GeoTIFF	https://www.gebco.net
ETOPO2	2 弧分	NetCDF	NOAA 官网
SRTM	1–3 弧秒	GeoTIFF	NASA Earthdata

本文以 GEBCO 2023（NetCDF 格式）为例，文件名为 gebco_2023_n50.0_s34.0_w73.0_e108.0.nc。

5.2 研究区边界矢量

准备研究区边界 Shapefile（如 NW_Arid_Region.shp），用于精确裁剪和描边。若无现成边界，可用以下方式近似

# 用经纬度范围+ cartopy 内置国界粗略框定# 精确制图仍需官方审图号边界数据

六、核心代码：从零绘制 DEM 地形图

以下代码分模块讲解，可直接拼接为完整脚本。

import numpy as npimport xarray as xrimport matplotlib.pyplot as pltimport matplotlib.colors as mcolorsfrom matplotlib.colors import LightSourceimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.io.shapereader import Readerimport warningswarnings.filterwarnings('ignore')

6.2 读取并裁剪 DEM 数据

# ========== 参数配置 ==========LON_MIN, LON_MAX = 73.0, 108.0LAT_MIN, LAT_MAX = 34.0, 50.0DEM_PATH = './gebco_2023_n50.0_s34.0_w73.0_e108.0.nc'BOUNDARY_PATH = './NW_Arid_Region.shp'  # 研究区边界OUTPUT_PATH = './NW_Arid_Region_DEM.png'# ========== 读取数据 ==========ds = xr.open_dataset(DEM_PATH)# 自动识别变量名（不同数据集命名不同）elevation_var = Nonefor candidate in ['elevation', 'z', 'altitude', 'band1']:    if candidate in ds.data_vars:        elevation_var = candidate        breakif elevation_var is None:    raise ValueError(f"未找到高程变量，可用变量: {list(ds.data_vars)}")dem = ds[elevation_var]# 裁剪研究区dem_clip = dem.sel(    lat=slice(LAT_MAX, LAT_MIN),   # 注意：纬度通常从北到南递减    lon=slice(LON_MIN, LON_MAX))lon = dem_clip.lon.valueslat = dem_clip.lat.valuesz = dem_clip.values# 高程截断（突出地形细节）z = np.clip(z, 0, 6000)

6.3 自定义地形色带

# ========== 构建 ArcGIS 风格地形色带 ==========terrain_colors = [    (0.20, 0.40, 0.20),   # 0-500m   深绿（河谷/绿洲）    (0.40, 0.65, 0.30),  # 500-1000m 浅绿（低山草原）    (0.75, 0.75, 0.50),  # 1000-2000m 黄绿（丘陵）    (0.90, 0.75, 0.40),  # 2000-3000m 土黄（戈壁）    (0.80, 0.55, 0.30),  # 3000-4000m 棕黄（荒漠）    (0.55, 0.35, 0.20),  # 4000-5000m 棕色（高山）    (0.40, 0.25, 0.20),  # 5000-5500m 深棕（极高山）    (0.90, 0.90, 0.95),  # >5500m    雪白色（冰川）]terrain_cmap = mcolors.LinearSegmentedColormap.from_list(    'arid_terrain', terrain_colors, N=256)

6.4 初始化地图投影

# ========== 地图投影设置 ==========# 等距圆柱投影（适合中低纬度区域概况图）proj = ccrs.PlateCarree()fig = plt.figure(figsize=(12, 8))ax = fig.add_subplot(1, 1, 1, projection=proj)# 设置显示范围ax.set_extent([LON_MIN, LON_MAX, LAT_MIN, LAT_MAX], crs=proj)# 添加地图要素ax.add_feature(cfeature.COASTLINE, linewidth=0.6, edgecolor='gray')ax.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle='--', edgecolor='gray')ax.add_feature(cfeature.RIVERS, alpha=0.4, edgecolor='blue', linewidth=0.4)ax.add_feature(cfeature.LAKES, alpha=0.5, facecolor='lightblue')

6.5 DEM 渲染 + 地形晕染

# ========== DEM 渲染 ==========im = ax.pcolormesh(    lon, lat, z,    cmap=terrain_cmap,    shading='auto',    vmin=0, vmax=6000,    transform=proj,    rasterized=True   # 栅格化，减小矢量输出体积)# ========== 地形晕染（Hillshade）增强立体感 ==========ls = LightSource(azdeg=315, altdeg=45)  # 光源：西北方向，45°仰角# 计算 hillshade（基于坡度坡向）hillshade = ls.hillshade(z, vert_exag=0.05)# 将 hillshade 作为透明度叠加层# 方法：用 RGBA 方式，把 hillshade 映射到 alpha 通道rgba = terrain_cmap((z - 0) / 6000)rgba[..., 3] = 0.9  # 基础透明度# 或者更高级的做法：用 shade_rgb 直接生成带光照的彩色地形rgb_shaded = ls.shade(    z, cmap=terrain_cmap, vert_exag=0.08,    blend_mode='soft', vmin=0, vmax=6000)ax.imshow(    rgb_shaded, extent=[LON_MIN, LON_MAX, LAT_MIN, LAT_MAX],    origin='lower', transform=proj, interpolation='bilinear')

6.6 叠加研究区边界

# ========== 叠加研究区边界 ==========# 方式1：从 Shapefile 读取try:    reader = Reader(BOUNDARY_PATH)    ax.add_geometries(        reader.geometries(), crs=proj,        facecolor='none', edgecolor='black', linewidth=1.5, zorder=10    )except Exception as e:    print(f"边界加载失败: {e}")# 方式2：用 cartopy 内置中国边界（粗略）# ax.add_feature(cfeature.BORDERS, linewidth=0.8)

6.7 添加经纬网格与标注

# ========== 经纬网格 ==========gl = ax.gridlines(    crs=proj, draw_labels=True,    linewidth=0.8, color='gray', alpha=0.5, linestyle='--',    xlocs=np.arange(73, 109, 5),   # 经度刻度    ylocs=np.arange(34, 51, 4)     # 纬度刻度)# 网格标签样式gl.top_labels = Falsegl.right_labels = Falsegl.xlabel_style = {'size': 11, 'fontname': 'Times New Roman'}gl.ylabel_style = {'size': 11, 'fontname': 'Times New Roman'}gl.xformatter = plt.matplotlib.ticker.FuncFormatter(    lambda x, _: f'{int(x)}°E')gl.yformatter = plt.matplotlib.ticker.FuncFormatter(    lambda y, _: f'{int(y)}°N')

6.8 添加 Colorbar、比例尺、指北针

# ========== Colorbar ==========cbar = plt.colorbar(    im, ax=ax, shrink=0.6, pad=0.02, aspect=25,    orientation='vertical', extend='max')cbar.set_label('Elevation (m)', fontsize=12, fontname='Times New Roman')cbar.ax.tick_params(labelsize=10)# ========== 比例尺 ==========def add_scalebar(ax, lon0, lat0, length_km, linewidth=3):    """在指定位置添加线段比例尺"""    # 1° 纬度 ≈ 111 km    length_deg = length_km / 111.0    ax.plot([lon0, lon0 + length_deg], [lat0, lat0],            'k-', linewidth=linewidth, solid_capstyle='butt')    ax.text(lon0 + length_deg/2, lat0 - 0.3, f'{length_km} km',            ha='center', va='top', fontsize=9, fontweight='bold')add_scalebar(ax, lon0=75, lat0=35.5, length_km=500)# ========== 指北针 ==========def add_north_arrow(ax, x, y, size=0.03):    """添加指北针"""    ax.annotate('N', xy=(x, y), xytext=(x, y-size*3),                arrowprops=dict(arrowstyle='->', color='black', lw=2),                fontsize=12, ha='center', va='center', fontweight='bold',                transform=ax.transAxes)  # 使用 axes 坐标系add_north_arrow(ax, x=0.92, y=0.88)# ========== 地名标注 ==========labels = [    (85.5, 46.5, '北疆'), (84.0, 39.0, '南疆'),    (100.0, 40.0, '河西走廊'), (82.0, 43.0, '天山'),    (82.0, 36.0, '昆仑山'), (88.0, 48.0, '阿尔泰山'),    (87.0, 41.0, '塔里木盆地'), (86.0, 47.0, '准噶尔盆地')]for lx, ly, txt in labels:    ax.text(lx, ly, txt, fontsize=9, fontweight='bold',            color='white' if ly < 40 else 'black',            ha='center', va='center',            path_effects=[plt.matplotlib.patheffects.withStroke(                linewidth=2, foreground='black'            )], zorder=15)# ========== 标题与输出 ==========ax.set_title('DEM Topographic Map of Northwest Arid Region, China',             fontsize=14, fontweight='bold', pad=15, fontname='Times New Roman')plt.tight_layout()plt.savefig(OUTPUT_PATH, dpi=300, bbox_inches='tight',            facecolor='white', edgecolor='none')plt.show()print(f"地图已保存至: {OUTPUT_PATH}")

七、效果预览

下图展示了 DEM 渲染 与 地形晕染（Hillshade） 两种风格的对比：

左图：标准 DEM 伪彩色渲染；右图：叠加 Hillshade 后的地形晕染效果，山体立体感显著增强，更适合论文插图。

八、进阶技巧

8.1 局部放大图（Inset Map）

博士论文常需要在主图一角添加局部放大图（如伊犁河谷、塔里木河下游）：

from mpl_toolkits.axes_grid1.inset_locator import inset_axes# 在主图内创建子轴ax_inset = inset_axes(ax, width="35%", height="35%",                      loc='lower left',                      bbox_to_anchor=(0.02, 0.02, 1, 1),                      bbox_transform=ax.transAxes)# 子图设置ax_inset.set_extent([80, 85, 42, 45], crs=proj)ax_inset.pcolormesh(lon, lat, z, cmap=terrain_cmap,                    vmin=0, vmax=6000, transform=proj)ax_inset.add_feature(cfeature.COASTLINE, linewidth=0.4)ax_inset.set_title('Ili Valley', fontsize=9)

8.2 叠加气象站点 / 采样点

# 站点坐标（示例）stations = {    'Urumqi': (87.6, 43.8),    'Kashgar': (76.0, 39.5),    'Zhangye': (100.4, 38.9)}for name, (slon, slat) in stations.items():    ax.plot(slon, slat, 'ro', markersize=5, transform=proj, zorder=15)    ax.annotate(name, (slon, slat), xytext=(5, 5),                textcoords='offset points', fontsize=8, color='darkred')

8.3 输出矢量格式（PDF/EPS）

期刊投稿常要求矢量图：

plt.savefig('./NW_Arid_Region_DEM.pdf', format='pdf',            bbox_inches='tight', dpi=300)

💡 提示：pcolormesh 默认生成栅格图像。若需矢量输出，可将 rasterized=True 去掉，但文件体积会暴增。建议底图栅格化、边界/标注矢量化的混合输出策略。

九、常见问题排查

报错信息	原因	解决方案
ModuleNotFoundError: No module named 'cartopy'	未安装或环境错误	conda install -c conda-forge cartopy
`ProjError` / 投影变换失败	PROJ 库版本冲突	用 conda 安装，避免 pip
ValueError: x and y arguments to pcolormesh	经纬度网格与数据维度不匹配	检查 `lon.shape`, `lat.shape`, `z.shape` 是否一致
地图空白/只显示网格	数据范围与 `set_extent` 不匹配	确认 `extent` 与 DEM 裁剪范围一致
色带不连续/有断层	`pcolormesh` 对 NaN 处理不当	用 `np.ma.masked_where(np.isnan(z), z)` 掩膜

十、推荐参考文献

如需在论文中引用西北干旱区背景，以下文献可供参考：

《中国西北干旱区水资源与生态环境研究报告》
The dominant warming season shifted from winter to spring in the arid region of Northwest China
Potential risks and challenges of climate change in the arid region of northwestern China
Spatiotemporal evolution and driving factors analysis of fractional vegetation coverage in the arid region of northwest China
Temporal and spatial changes of extreme precipitation and its related large-scale climate mechanisms in the arid region of Northwest China during 1961–2022
《伊犁河谷草地水分利用效率时空变化及其对气候因子的响应》