当前位置: 首页 > news >正文

cartopy 绘制中国地图:从基础边界到南海诸岛与十段线的完整实践

1. 为什么选择cartopy绘制中国地图

第一次接触地理数据可视化时,我尝试过Basemap、GeoPandas等多种工具,最终发现cartopy在绘制中国地图时有着不可替代的优势。这个由英国气象局开发的库,不仅继承了Basemap的易用性,还解决了投影变换时的诸多痛点。特别是在处理中国版图这种包含特殊区域的场景时,cartopy的灵活性让我印象深刻。

记得去年做气象数据可视化项目时,我需要同时展示大陆、台湾地区、南海诸岛和十段线的完整版图。当时试了三个不同的shp文件,不是缺少藏南地区,就是台湾地区的边界不完整。后来在GMT中文社区找到的CN-border-La.dat文件,配合Meteo-Python群的shp资源,才终于解决了这个问题。这也让我意识到,选择正确的数据源和工具组合有多么重要。

cartopy的核心优势在于它原生支持多种地图投影,并且能够智能处理坐标转换。比如PlateCarree投影虽然简单,但在展示跨越180度经线的区域时容易出错。这时只需要设置central_longitude参数,就能完美解决经线跳跃的问题。我在实际项目中测试过,从70°E到140°E的中国全图范围,用PlateCarree(central_longitude=180)投影效果最稳定。

2. 数据准备与常见陷阱规避

2.1 获取可靠的边界数据

绘制准确的中国地图,第一步就是找到完整可靠的边界数据。经过多次踩坑,我总结出几个可信的数据源:

  • GMT中文社区提供的CN-border-La.dat文件,包含完整的陆地边界和十段线
  • Meteo-Python技术交流群分享的china_country.shp文件
  • Natural Earth Data的1:50m行政边界数据(需要额外验证完整性)

这里有个重要提示:千万不要随便下载来路不明的shp文件。我曾经用一个国外网站下载的"China_Boundary.shp",结果发现缺少了钓鱼岛等重要区域。后来对比发现,有些数据源为了所谓的"简化",会刻意省略敏感区域,这在专业项目中是绝对不可接受的。

2.2 处理dat格式的边界数据

CN-border-La.dat文件虽然完整,但直接读取会遇到两个问题:

  1. 文件包含注释行(以#开头)
  2. 边界段用">"符号分隔

这是我处理这类文件的代码模板:

with open('CN-border-La.dat') as src: # 过滤注释行 context = ''.join([line for line in src if not line.startswith('#')]) # 按>分割边界段 blocks = [cnt for cnt in context.split('>') if len(cnt) > 0] # 转换为numpy数组 borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]

2.3 加载shp文件的正确姿势

cartopy支持两种加载shp文件的方式,经过实测都可用:

方法一:使用ShapelyFeature

shape_feature = cfeature.ShapelyFeature( shpfile.geometries(), ccrs.PlateCarree(), facecolor='teal') ax.add_feature(shape_feature)

方法二:遍历记录逐个添加

for rec in shpfile.records(): ax.add_geometries( [rec.geometry], crs=ccrs.PlateCarree(), facecolor='red')

特别注意:shp文件会覆盖cartopy自带的河流湖泊等要素,所以一定要先添加基础地理要素,最后再加载shp文件。这个顺序问题曾经让我调试了整整一个下午。

3. 完整绘图流程详解

3.1 初始化地图画布

创建画布时需要特别注意两点:

  1. 设置足够大的figsize(建议至少8×6)
  2. 选择合适的投影方式
fig = plt.figure(figsize=[8, 5.5]) ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))

3.2 添加基础地理要素

按顺序添加这些要素可以避免图层覆盖问题:

  1. 海洋和陆地背景
  2. 国界线和海岸线
  3. 河流和湖泊
  4. 自定义shp文件
  5. 十段线等特殊边界
ax.add_feature(cfeature.OCEAN.with_scale('50m')) ax.add_feature(cfeature.LAND.with_scale('50m')) ax.add_feature(cfeature.BORDERS.with_scale('50m'), lw=0.5) ax.add_feature(cfeature.COASTLINE.with_scale('50m'), lw=0.5)

3.3 绘制十段线和特殊边界

使用之前处理好的borders数据绘制完整边界:

for line in borders: ax.plot(line[0::2], line[1::2], '-', lw=1, color='k', transform=ccrs.Geodetic())

这里有个关键细节:transform参数必须设置为ccrs.Geodetic(),这样才能确保线段在不同投影下正确显示。

3.4 配置经纬度网格

网格线的配置相对复杂,但掌握这几个参数就能应对大多数场景:

gl = ax.gridlines( crs=ccrs.PlateCarree(), draw_labels=False, linewidth=1.2, color='k', alpha=0.5, linestyle='--') # 设置标签显示位置 gl.xlabels_top = False gl.ylabels_right = False # 配置格式化器 gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER # 设置网格间隔 gl.xlocator = mticker.FixedLocator(np.arange(70, 140+10, 10)) gl.ylocator = mticker.FixedLocator(np.arange(0, 55+10, 10))

4. 常见问题解决方案

4.1 跨越180度经线的显示问题

中国地图的经度范围(70°E-140°E)正好跨越了180度经线附近,这会导致地图在显示时出现断裂。解决方案是设置central_longitude参数:

# 错误的做法 - 会出现地图断裂 ax = plt.axes(projection=ccrs.PlateCarree()) # 正确的做法 ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))

4.2 台湾地区显示不完整

遇到台湾地区显示不全的情况,通常有两个原因:

  1. 使用的shp文件本身数据不全
  2. 地图范围设置不当

确保set_extent包含完整范围:

ax.set_extent([70, 140, 0, 55], crs=ccrs.PlateCarree())

4.3 南海诸岛比例失调

南海诸岛在小比例尺地图上常常显示过小,可以通过以下方式优化:

  1. 使用inset_axes添加放大的南海小图
  2. 调整主图范围,给南海留出更多空间
  3. 使用不同的投影方式,如AlbersEqualArea

4.4 输出图像分辨率不足

在保存图像时,建议:

  1. 设置dpi≥300
  2. 使用矢量格式(如PDF、SVG)
  3. 调整figsize确保细节清晰
plt.savefig('china_map.png', dpi=300, bbox_inches='tight')

5. 进阶技巧与优化建议

5.1 添加省界和城市标记

在完成国界绘制后,可以进一步添加省界和主要城市:

# 加载省界shp文件 province_shp = shapereader.Reader('china_province.shp') for rec in province_shp.records(): ax.add_geometries( [rec.geometry], crs=ccrs.PlateCarree(), edgecolor='gray', facecolor='none')

5.2 自定义样式美化

通过调整这些参数可以让地图更专业:

  • 海岸线线宽(lw参数)
  • 陆地填充色(facecolor)
  • 边界线样式(linestyle)
  • 透明度(alpha)
ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor='#f0f0f0', edgecolor='none')

5.3 添加比例尺和指北针

虽然cartopy没有内置的比例尺工具,但可以通过matplotlib实现:

# 添加比例尺 scale_bar(ax, 1000) # 需要自定义scale_bar函数 # 添加指北针 ax.annotate('N', xy=(0.9, 0.9), xycoords='axes fraction', fontsize=20, ha='center')

5.4 性能优化技巧

当处理高精度数据时,可以尝试:

  1. 使用simplify属性降低shp文件复杂度
  2. 设置合适的map_extent减少渲染范围
  3. 使用rasterize处理密集点数据
from matplotlib.path import Path from matplotlib.patches import PathPatch for rec in shpfile.records(): geom = rec.geometry.simplify(0.01) # 简化几何图形 ax.add_geometries( [geom], crs=ccrs.PlateCarree(), facecolor='red')

6. 完整代码示例与效果展示

将上述所有步骤整合,这里给出一个完整的绘图示例:

# coding:utf-8 import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER import cartopy.io.shapereader as shapereader import matplotlib.ticker as mticker # 加载dat格式边界数据 with open('CN-border-La.dat') as src: context = ''.join([line for line in src if not line.startswith('#')]) blocks = [cnt for cnt in context.split('>') if len(cnt) > 0] borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks] # 加载shp文件 shpfile = shapereader.Reader('china_country.shp') # 创建画布 fig = plt.figure(figsize=[10, 7]) ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) # 添加基础地理要素 ax.add_feature(cfeature.OCEAN.with_scale('50m')) ax.add_feature(cfeature.LAND.with_scale('50m')) ax.add_feature(cfeature.BORDERS.with_scale('50m'), lw=0.5) ax.add_feature(cfeature.COASTLINE.with_scale('50m'), lw=0.5) # 绘制shp文件 for rec in shpfile.records(): ax.add_geometries([rec.geometry], crs=ccrs.PlateCarree(), facecolor='#f5deb3', edgecolor='gray') # 绘制十段线等边界 for line in borders: ax.plot(line[0::2], line[1::2], '-', lw=1, color='k', transform=ccrs.Geodetic()) # 添加河流湖泊(必须在shp之后) ax.add_feature(cfeature.RIVERS.with_scale('50m')) ax.add_feature(cfeature.LAKES.with_scale('50m')) # 配置网格 gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--') gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER gl.xlocator = mticker.FixedLocator(np.arange(70, 140+10, 10)) gl.ylocator = mticker.FixedLocator(np.arange(0, 55+10, 10)) # 设置显示范围 ax.set_extent([70, 140, 0, 55], crs=ccrs.PlateCarree()) plt.title('中国全图(含南海诸岛和十段线)') plt.savefig('china_full_map.png', dpi=300, bbox_inches='tight') plt.show()

最终效果图中,你可以看到:

  1. 完整的中国陆地边界
  2. 清晰的台湾地区轮廓
  3. 南海诸岛和十段线的准确位置
  4. 协调的色彩搭配和专业的标注

在实际项目中,这套代码框架已经成功应用于气象数据可视化、地理信息系统开发等多个领域。根据具体需求,你可以灵活调整颜色方案、添加数据叠加层,或者集成到更复杂的可视化流程中。

http://www.cnnetsun.cn/news/2587225.html

相关文章:

  • 5分钟学会B站缓存视频转换:永久保存你收藏的珍贵内容
  • Linux---进程(概念,PCB,进程属性,标示符,fork)
  • RAG 高级技术与调优实战手册
  • 自治系统失控:从故障模式到抗错设计的工程实践
  • 构建稳健AI应用:隔离、容错与可观测性架构设计实践
  • pypto:用Python直接写NPU算子,门槛有多低?
  • 保姆级教程:用RDPWrap解锁Win10/11家庭版远程桌面,还能多人同时登录
  • 告别混乱状态机!用UE4行为树+黑板实现智能敌人AI(实战案例解析)
  • Unity 2022.3.3 LTS + Visual Studio 2022:手把手教你复刻《吸血鬼幸存者》核心战斗(附完整源码)
  • Taotoken模型广场首发更新Qwen与Gemini等旗舰模型体验
  • 模型评测为什么一上对抗攻击测试就开始高分低防御:从 Adversarial Prompt 到 Robustness Budget 的工程实战
  • 淘宝任务自动化终极指南:5分钟解放双手的免费淘金币脚本
  • “襄阳造”打磨车出口毛里塔尼亚
  • 贝叶斯双重机器学习:高维因果推断的去偏与不确定性量化
  • Claude Code VS Code扩展:AI编程代理的工程化实践
  • TikTok 短视频生成工具哪家好?爆款视频复刻工具实用推荐
  • Godot PCK文件结构解析与安全解包实战指南
  • sqlmap原理深度解析:从DVWA靶场看SQL注入本质
  • 机器学习辅助高通量筛选:uMLIP与迁移学习加速功能材料发现
  • GBase 8s数据库常见问题排查及解决方法简述
  • 机器学习与模拟退火优化布尔特征集变量排序,加速密码分析计算
  • Unity Hub安装Android组件失败的真相与三步修复法
  • 大厂级AI服务对接实战(OpenAI/Anthropic/Claude全栈集成手册)
  • Unity工控机HMI开发实战:从协议接入到工业级部署
  • 开源免费!这款 AI 语音工作室让 ElevenLabs 都感到压力
  • 模拟实现:glibc_1.0-文件操作函数fopen fclose fwrite fflush实现。
  • 零样本与开放词汇目标检测:从语义对齐到开放世界感知的技术演进与实践
  • 别再手动折腾了!用Docker Compose一键部署Yapi接口管理平台(附完整配置文件)
  • AR物体识别抖动原理与四层实战优化方案
  • Unity Shader Graph溶解特效的物理建模与多尺度实现