画北极区域的图说难不难说简单也真不简单。最近有个朋友问我用python的cartopy怎么画“北极部分区域”的扇形图他想把特定经度范围内的极区海冰数据展示出来而不是每次都出全场南半球那种大圆饼。这个问题特别典型我干脆把整个思路、代码、以及我踩过的坑都整理出来希望能帮到想画极区扇形图的人。先说结论用cartopy做扇形区域图核心就三件事——选对投影、控制范围、用set_boundary做精确裁剪。实测下来这一套在北极海冰分析、航道研究、气象数据可视化里都非常实用而且不需要多复杂的数学基础只要理解几个坐标系概念半小时就能跑出自己的第一张扇形图。1. 扇形区域图到底难在哪先想清楚再动手1.1 用普通地图投影画北极会发生什么我自己第一次画北极的时候直接用了最常见的PlateCarree等距圆柱投影就是那种把地球直接摊开的平面经纬度地图。结果可想而知高纬度的格陵兰岛和加拿大北极群岛被横向拉伸得又扁又宽北极点附近的数据全部挤成一团根本看不出空间分布规律。这种现象叫“极区畸变”本质是球面坐标强行映射到平面时经线在高纬度区域的间距被严重压缩视觉上给人一种“北冰洋比实际大好几倍”的错觉。所以在画极区图之前必须搞清楚一个基本逻辑不是选“看起来好看”的投影而是选“在该区域变形最小”的投影。北极图最常用的方案是极射赤面投影NorthPolarStereo它的特点是投影中心在北极点变形主要发生在靠近赤道的区域而极区内部的角度和面积关系保持得比较好。这也是为什么气象、海洋领域做极区分析默认都是这类投影。1.2 “扇形”的本质画得好不好全看投影与裁剪理解了投影问题后就要说“扇形”这个词。很多人口中的扇形区域其实包含了两种需求一种是只画一个经纬度范围的矩形区域比如东经30度到150度、北纬60度到90度另一种是真正意义上“带角度扇面”的图也就是从北极点向外辐射的两条经线加上一条纬度圆弧围出来的饼状区域。这两种需求在cartopy里的实现方式完全不一样。第一种非常容易用set_extent指定范围就行第二种就需要自己构造一条闭合的扇形边界路径再通过set_boundary做裁剪。更麻烦的是如果你直接在北半球极射投影下用经纬度矩形去做set_extent出来的图边缘是直线切割的它看起来像一块“被切过的圆饼”而不是标准的“扇面”。具体要哪种得在动手前想清楚是数据范围从视觉上本来就是扇形还是你希望出图后裁成扇形。2. 投影选型与坐标系理解别在这步犯错2.1 三种常见投影方案对比cartopy支持的投影非常多但画北极部分区域常用的其实就那几款我自己实际用下来做了一个对比。投影名称cartopy类名适合场景注意事项北半球极射赤面投影NorthPolarStereo北极区域海冰、气象要素图纬度越低变形越大一般只画50度以北兰伯特等积方位投影LambertAzimuthalEqualArea极区面积对比、统计类分析面积关系精准但形状在高纬边缘有拉伸等距圆柱投影PlateCarree全图概览低纬度区域极区严重变形不适合北极局部图我个人的建议是如果你只是展示北冰洋附近的海冰范围或温度分布NorthPolarStereo是最稳妥的选择如果你要做面积统计或者强调“哪个区域海冰少了多少”那LambertAzimuthalEqualArea在等积性上更专业。至于PlateCarree除非你只是快速预览数据否则尽量别用在北极上。2.2 理解cartopy里的crs与transform很多刚接触cartopy的人会被“crs”和“transform”这两个概念搞晕尤其当代码里一会儿出现ccrs.NorthPolarStereo一会儿出现ccrs.PlateCarree还有人跟你说“数据要指定transformccrs.PlateCarree()”这时候很容易蒙。我打个比方可以把投影坐标系crs理解成“画布”的透视方式而transform是“画笔”在落笔前用的坐标系统。画布是什么投影决定了图底怎么画但你的数据是经纬度网格在画布上落笔前必须告诉cartopy这些经纬度属于哪个坐标系它才能换算到画布坐标里。实战中最常见的数据是经纬度网格而画布经常是NorthPolarStereo所以代码就写成这样ax fig.add_subplot(projectionccrs.NorthPolarStereo()) mesh ax.pcolormesh(lon, lat, data, transformccrs.PlateCarree())记住一个口诀画布用目标投影数据源坐标用transform指定。你只要保证数据的transform和画布投影两者定义正确cartopy会自动完成换算不需要自己手动做坐标变换。3. 从矩形范围到真扇形边界构造与裁剪3.1 先用set_extent搞定“部分区域”如果你只需要画一个经纬度范围内的“北极局部图”那就不用折腾扇形边界了一个set_extent就能解决。以我常用的北极海冰分析为例我想看北欧海和巴伦支海一带范围大概是西经20度向东到东经80度、北纬65度到90度。代码很简单import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature proj ccrs.NorthPolarStereo(central_longitude30) fig, ax plt.subplots(figsize(8, 8), subplot_kw{projection: proj}) ax.set_extent([-45, 80, 60, 90], crsccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linestyle:, linewidth0.5) ax.gridlines(draw_labelsFalse, linestyle--, alpha0.6) plt.show()这里有两个容易踩坑的点。第一set_extent的前四个参数是[min_lon, max_lon, min_lat, max_lat]但经纬度本身要指定为PlateCarree坐标不然会报错或者裁剪位置完全不对。第二central_longitude建议设置在你要展示区域的中心附近这样图的展开范围更自然不会把你关心的区域劈成两半。实测中如果把central_longitude设成0而你想看东经100到160的区域图会被整体旋转一个角度视觉上不够直观。3.2 用set_boundary构造真正的扇形如果是想要真正的扇形区域set_extent就不够了因为它只能做矩形或椭圆的裁剪窗口没办法画一个“弧底两条半径”的扇面。这个时候需要用matplotlib.path下的Path来构造一条闭合边界然后用GeoAxes.set_boundary把它裁出来。我自己最常用的一种方法是在经纬度平面里定义扇形边界上的采样点然后转成投影坐标。这里的关键是要把扇形的顶点放在北极点沿两条经线向外扩展到底部纬度最后沿底边纬线画一段圆弧。代码如下import numpy as np import matplotlib.pyplot as plt import matplotlib.path as mpath import cartopy.crs as ccrs import cartopy.feature as cfeature proj ccrs.NorthPolarStereo(central_longitude45) pc ccrs.PlateCarree() # 扇形参数 lon_min, lon_max 0, 90 lat_min 60 # 在经纬度坐标内构造扇形边界 arc_lons np.linspace(lon_min, lon_max, 200) arc_lats np.full_like(arc_lons, lat_min) lon_seq [lon_min, lon_min] list(arc_lons) [lon_max, lon_max] lat_seq [90, lat_min] list(arc_lats) [lat_min, 90] # 转成NorthPolarStereo投影坐标 verts proj.transform_points(pc, np.array(lon_seq), np.array(lat_seq)) sector_path mpath.Path(verts[:, :2]) fig, ax plt.subplots(figsize(8, 8), subplot_kw{projection: proj}) ax.set_boundary(sector_path) ax.set_extent([-10, 100, 60, 90], crspc) ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linestyle:, linewidth0.5) ax.gridlines(linestyle--, alpha0.6, draw_labelsFalse) plt.show()重点解释一下这段代码的逻辑。lon_seq和lat_seq里我先把起点放在北极点纬度90然后让它沿着lon_min这条经线一路往下到底部纬度lat_min再沿着纬线lat_min从lon_min平滑地走到lon_max最后沿着lon_max经线回到北极点。这样一整圈围出来的就是一个标准的扇面。之所以要采样200个点来做底部圆弧是因为如果只有3个点投影到NorthPolarStereo之后底部就是一条折线视觉上很生硬。采样点越多弧度越圆润。这里选择在经纬度坐标内采样再由transform_points投影到目标坐标系是最不容易出错的方式因为它保证了每两个相邻采样点之间的连线在地理意义上是合理的。3.3 边界裁剪的常见误区用set_boundary的时候我最开始犯过一个错误以为传进去的路径必须在画布投影坐标下定义所以就自己噼里啪啦算了一堆立体几何公式结果不是缝隙就是变形。后来才意识到最简单稳妥的办法就是先在经纬度坐标里构造好路径再交给transform_points转换这个步骤属于常见实践强烈建议直接抄作业。另一个常见误区是set_boundary和set_extent同时使用时的顺序问题。我习惯先set_boundary再set_extent这样既能裁剪扇面又能保证地图的可视范围不会无限延伸到北冰洋以外的区域。如果只set_boundary不set_extent图的范围可能会被自动扩展到很远的地方视觉上扇面会很小。如果只set_extent不set_boundary那就还是一个矩形的裁剪窗口不是扇形。4. 完整实操从底图到数据叠加一步不落4.1 环境准备与基础底图开始跑代码之前先确认环境里装了这几个库python本体之外还需要numpy、matplotlib、cartopy。cartopy的安装偶尔会让你头疼如果使用conda环境会顺畅一些。以我自己的经验到了新电脑上我最常用下面这种方式创建环境conda create -n cartopy_env python3.10 conda activate cartopy_env conda install -c conda-forge cartopy matplotlib numpy装好之后先跑一个最基础的底图脚本确认coastline能正常加载。如果运行时提示“DownloadError”或者“OError: [Errno -2] Name or service not known”多半是海岸线数据下载失败。解决办法有两个一是检查网络二是用cfeature.NaturalEarthFeature的cache目录手动放置下载好的矢量数据。我建议新人在离线或网络不好的环境下直接改用ax.add_feature(cfeature.OCEAN, colorlightblue) ax.add_feature(cfeature.LAND, colorlightgray)这个不依赖在线下载能保证基础底图先跑通后面再慢慢补海岸线细节。4.2 数据图层以海冰密集度为例底图画好了接下来就是往图里叠真正的数据。我最常画的是海冰密集度场数据格式通常是一个二维数组配合二维的经度网格和纬度网格坐标原点在北极附近。叠加数据时一定要记得加transformccrs.PlateCarree()否则数据会被错误地当成投影坐标画到图上结果就是你看到海冰分布全都挤在一个角落。下面是完整示例假设你已经读取了数据lon2d和lat2d是二维经纬度网格sic是海冰密集度数组0到1之间的小数import numpy as np import matplotlib.pyplot as plt import matplotlib.path as mpath import matplotlib.ticker as mticker import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter # --- 扇形参数 --- lat_min 60 lon_min, lon_max -90, 60 # --- 投影设置 --- proj ccrs.NorthPolarStereo(central_longitude-30) pc ccrs.PlateCarree() # --- 构建扇形边界 --- arc_lons np.linspace(lon_min, lon_max, 200) arc_lats np.full_like(arc_lons, lat_min) lon_seq [lon_min, lon_min] list(arc_lons) [lon_max, lon_max] lat_seq [90, lat_min] list(arc_lats) [lat_min, 90] verts proj.transform_points(pc, np.array(lon_seq), np.array(lat_seq)) sector_path mpath.Path(verts[:, :2]) # --- 画布 --- fig plt.figure(figsize(9, 9)) ax fig.add_subplot(1, 1, 1, projectionproj) ax.set_boundary(sector_path) ax.set_extent([lon_min - 20, lon_max 20, lat_min - 5, 90], crspc) # --- 底图要素 --- ax.add_feature(cfeature.LAND, colorlightgray, zorder2) ax.add_feature(cfeature.OCEAN, colorwhite, zorder1) ax.coastlines(resolution50m, linewidth0.7, colorblack, zorder3) # --- 模拟一份海冰数据实际使用时替换成自己的数组 --- lons_grid np.linspace(lon_min, lon_max, 120) lats_grid np.linspace(lat_min, 88, 100) lon2d, lat2d np.meshgrid(lons_grid, lats_grid) sic np.clip(np.cos(np.deg2rad(lat2d - 55)) * 0.8 0.1, 0, 1) # --- 画数据 --- mesh ax.pcolormesh(lon2d, lat2d, sic, transformpc, cmapSpectral_r, shadingauto, zorder4) cbar fig.colorbar(mesh, axax, shrink0.7, pad0.05) cbar.set_label(Sea Ice Concentration) # --- 网格线 --- gl ax.gridlines(crspc, linewidth0.5, colorgray, alpha0.6, linestyle--, draw_labelsTrue, xlocsrange(-180, 181, 30), ylocsrange(0, 91, 10)) gl.top_labels False gl.right_labels False gl.xformatter LongitudeFormatter() gl.yformatter LatitudeFormatter() plt.show()这段代码里我做了几件比较重要的事一是把扇形边界通过set_boundary固定住让图层无论是底图还是数据都会被裁到扇形之内二是用pcolormesh叠加数据并用Spectral_r配色展示海冰密集度暖色代表密集度高冷色代表密集度低三是用gridlines的xlocs和ylocs手动控制了经纬度网格线的密度避免线条太密糊成一团。这里需要提醒一下pcolormesh的shadingauto参数是matplotlib比较新的特性能自动判断是网格还是散点旧版本可能需要改成shadingnearest或gouraud。如果你的matplotlib版本比较老遇到“shading auto”不兼容的报错直接改掉就行。4.3 网格线、标注与出图设置网格线的设置看起来是小事但在极区图上特别容易翻车。如果你用默认的gridlines()去画北极扇形图经常会出现经纬度标注挤在扇面周围的奇怪位置上甚至有些标签飘到边界外面。比较好的做法是像我上面那样给gridlines指定crsccrs.PlateCarree()并手动指定经纬度刻度位置然后用LongitudeFormatter和LatitudeFormatter来格式化标签。还要注意top_labels和right_labels在极区图里通常设为False不然图上面和右面的标注会和左边的撞在一起。我自己的习惯是只显示左边纬度和下边经度不但清爽也符合大多数论文图表的阅读习惯。最后是出图设置。扇形图建议用正方形画布也就是figsize(8, 8)或者(9, 9)因为扇形本身在NorthPolarStereo投影下就是一个以中心为顶点的圆饼局部正方形画布能最大程度利用空间。保存时用dpi300格式用png或者pdf都行。如果发给导师或者写报告用我一般保存pdf——矢量格式缩放不失真投到论文里也够清晰fig.savefig(arctic_sector.png, dpi300, bbox_inchestight) fig.savefig(arctic_sector.pdf, bbox_inchestight)5. 高频报错与排错经验能帮你省两小时5.1 问题速查表下面这几类问题是我在画北极扇形图时遇到过也在网上帮别人排查过的。这里整理成表格遇到报错可以对号入座。问题现象可能原因解决思路add_feature时报错或海岸线不出来自然地球矢量数据下载失败先检查网络或者改用OCEAN/LAND基础色块画出的图是一片空白没有数据数据范围不在投影可见范围内检查set_extent范围确认数据和extent是否有交集set_boundary后经纬网格线消失裁剪边界格式不对或transform缺失确认path是投影坐标且不是DataFrame数组set_extent报错“ValueError: CRS mismatch”范围参数没有指定crs写成ax.set_extent([...], crsccrs.PlateCarree())图形出现一条从顶部穿过的斜线多边形路径没有闭合或者采样点太少检查lon_seq和lat_seq首尾是否回到北极点图上横坐标刻度太密集标签挤在一起gridlines自动生成过多刻度手动指定xlocs、ylocs或者用MaxNLocator控制5.2 经度跨边界-180与180之间那道裂缝这个坑非常容易遇到。如果你要画的区域跨过了日期变更线比如从东经170度到西经160度直接用一次linspace得到的是一个横跨整个太平洋的大弧形画出来完全不是你要的区域。cartopy在处理这种跨0度或跨180度的数据时尤其需要小心。解决办法有两个。推荐的做法是把经度范围拆成两段分别构造后再合并。比如要画东经170到西经160也就是在坐标上是从170到200或从-190到-160那你可以统一用“-190到-160”的连续区间来表达lon_min-190, lon_max-160。NorthPolarStereo投影本身不关心经度值是正还是负只要数据坐标与边界线段能不解纠缠就行。如果你的数据本身经度是0到360的格式更好办直接统一到“0到360”范围内计算。我的经验是在做扇形区域图之前先把所有经度统一成同一种习惯要么都是-180到180要么都是0到360。千万别混用不然边界路径和数据图层各画各的出了图就再也对不上了。5.3 海岸线、数据错位与画图横坐标过密的通用解法错位问题多半是transform没写对。比如海岸线是cartopy内置的它知道自己的坐标系统但你自己的海冰数据一定要在pcolormesh里指定transform。如果数据出处不明、网格本身就是投影坐标那就先把数据插值到经纬度网格再画不要偷懒直接塞进pcolormesh。至于热搜里经常有人问“python画图横坐标太密集”这个在cartopy里的表现就是经纬度网格线标签叠在一起。我自己的做法是给gridlines传入xlocs和ylocs手动指定要显示哪些刻度。比如在北极扇形图中我通常让经度每隔30度、纬度每隔10度显示一个标签清晰又不拥挤。如果是普通的matplotlib图则可以使用matplotlib.ticker的MaxNLocator或MultipleLocator控制刻度数量。核心思路是一样的别让绘图库自动决定刻度密度自己动手定好间隔一劳永逸。如果再遇到那些莫名其妙的报错我建议把代码逐步注释掉从底图开始一点点加要素。哪一步加了之后图形开始异常问题就锁定在哪一步。这个方法虽然笨但非常有效。我自己做了这么多次极区图之后最大的体会是扇形图看着复杂本质上并不是某种神奇操作而是把“投影”、“边界路径”、“数据transform”这三个链路理清楚剩下的只是往图里一层一层叠要素。只要你对“经纬度坐标出发、由投影转换到画布坐标”这条主线有把握不管换成北半球还是南半球还是任意一块“切角区域”都能照葫芦画瓢地做出来。如果你正准备画一块扇形区域但又不知道怎么处理边界可以直接把我上面的扇形路径构造函数摘出去改成你要的经度范围和底边纬度就行。这样省下来的时间拿去调配色和排版它不香吗
Cartopy极区扇形图绘制指南：投影、边界裁剪与实战示例
画北极区域的图说难不难说简单也真不简单。最近有个朋友问我用python的cartopy怎么画“北极部分区域”的扇形图他想把特定经度范围内的极区海冰数据展示出来而不是每次都出全场南半球那种大圆饼。这个问题特别典型我干脆把整个思路、代码、以及我踩过的坑都整理出来希望能帮到想画极区扇形图的人。先说结论用cartopy做扇形区域图核心就三件事——选对投影、控制范围、用set_boundary做精确裁剪。实测下来这一套在北极海冰分析、航道研究、气象数据可视化里都非常实用而且不需要多复杂的数学基础只要理解几个坐标系概念半小时就能跑出自己的第一张扇形图。1. 扇形区域图到底难在哪先想清楚再动手1.1 用普通地图投影画北极会发生什么我自己第一次画北极的时候直接用了最常见的PlateCarree等距圆柱投影就是那种把地球直接摊开的平面经纬度地图。结果可想而知高纬度的格陵兰岛和加拿大北极群岛被横向拉伸得又扁又宽北极点附近的数据全部挤成一团根本看不出空间分布规律。这种现象叫“极区畸变”本质是球面坐标强行映射到平面时经线在高纬度区域的间距被严重压缩视觉上给人一种“北冰洋比实际大好几倍”的错觉。所以在画极区图之前必须搞清楚一个基本逻辑不是选“看起来好看”的投影而是选“在该区域变形最小”的投影。北极图最常用的方案是极射赤面投影NorthPolarStereo它的特点是投影中心在北极点变形主要发生在靠近赤道的区域而极区内部的角度和面积关系保持得比较好。这也是为什么气象、海洋领域做极区分析默认都是这类投影。1.2 “扇形”的本质画得好不好全看投影与裁剪理解了投影问题后就要说“扇形”这个词。很多人口中的扇形区域其实包含了两种需求一种是只画一个经纬度范围的矩形区域比如东经30度到150度、北纬60度到90度另一种是真正意义上“带角度扇面”的图也就是从北极点向外辐射的两条经线加上一条纬度圆弧围出来的饼状区域。这两种需求在cartopy里的实现方式完全不一样。第一种非常容易用set_extent指定范围就行第二种就需要自己构造一条闭合的扇形边界路径再通过set_boundary做裁剪。更麻烦的是如果你直接在北半球极射投影下用经纬度矩形去做set_extent出来的图边缘是直线切割的它看起来像一块“被切过的圆饼”而不是标准的“扇面”。具体要哪种得在动手前想清楚是数据范围从视觉上本来就是扇形还是你希望出图后裁成扇形。2. 投影选型与坐标系理解别在这步犯错2.1 三种常见投影方案对比cartopy支持的投影非常多但画北极部分区域常用的其实就那几款我自己实际用下来做了一个对比。投影名称cartopy类名适合场景注意事项北半球极射赤面投影NorthPolarStereo北极区域海冰、气象要素图纬度越低变形越大一般只画50度以北兰伯特等积方位投影LambertAzimuthalEqualArea极区面积对比、统计类分析面积关系精准但形状在高纬边缘有拉伸等距圆柱投影PlateCarree全图概览低纬度区域极区严重变形不适合北极局部图我个人的建议是如果你只是展示北冰洋附近的海冰范围或温度分布NorthPolarStereo是最稳妥的选择如果你要做面积统计或者强调“哪个区域海冰少了多少”那LambertAzimuthalEqualArea在等积性上更专业。至于PlateCarree除非你只是快速预览数据否则尽量别用在北极上。2.2 理解cartopy里的crs与transform很多刚接触cartopy的人会被“crs”和“transform”这两个概念搞晕尤其当代码里一会儿出现ccrs.NorthPolarStereo一会儿出现ccrs.PlateCarree还有人跟你说“数据要指定transformccrs.PlateCarree()”这时候很容易蒙。我打个比方可以把投影坐标系crs理解成“画布”的透视方式而transform是“画笔”在落笔前用的坐标系统。画布是什么投影决定了图底怎么画但你的数据是经纬度网格在画布上落笔前必须告诉cartopy这些经纬度属于哪个坐标系它才能换算到画布坐标里。实战中最常见的数据是经纬度网格而画布经常是NorthPolarStereo所以代码就写成这样ax fig.add_subplot(projectionccrs.NorthPolarStereo()) mesh ax.pcolormesh(lon, lat, data, transformccrs.PlateCarree())记住一个口诀画布用目标投影数据源坐标用transform指定。你只要保证数据的transform和画布投影两者定义正确cartopy会自动完成换算不需要自己手动做坐标变换。3. 从矩形范围到真扇形边界构造与裁剪3.1 先用set_extent搞定“部分区域”如果你只需要画一个经纬度范围内的“北极局部图”那就不用折腾扇形边界了一个set_extent就能解决。以我常用的北极海冰分析为例我想看北欧海和巴伦支海一带范围大概是西经20度向东到东经80度、北纬65度到90度。代码很简单import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature proj ccrs.NorthPolarStereo(central_longitude30) fig, ax plt.subplots(figsize(8, 8), subplot_kw{projection: proj}) ax.set_extent([-45, 80, 60, 90], crsccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linestyle:, linewidth0.5) ax.gridlines(draw_labelsFalse, linestyle--, alpha0.6) plt.show()这里有两个容易踩坑的点。第一set_extent的前四个参数是[min_lon, max_lon, min_lat, max_lat]但经纬度本身要指定为PlateCarree坐标不然会报错或者裁剪位置完全不对。第二central_longitude建议设置在你要展示区域的中心附近这样图的展开范围更自然不会把你关心的区域劈成两半。实测中如果把central_longitude设成0而你想看东经100到160的区域图会被整体旋转一个角度视觉上不够直观。3.2 用set_boundary构造真正的扇形如果是想要真正的扇形区域set_extent就不够了因为它只能做矩形或椭圆的裁剪窗口没办法画一个“弧底两条半径”的扇面。这个时候需要用matplotlib.path下的Path来构造一条闭合边界然后用GeoAxes.set_boundary把它裁出来。我自己最常用的一种方法是在经纬度平面里定义扇形边界上的采样点然后转成投影坐标。这里的关键是要把扇形的顶点放在北极点沿两条经线向外扩展到底部纬度最后沿底边纬线画一段圆弧。代码如下import numpy as np import matplotlib.pyplot as plt import matplotlib.path as mpath import cartopy.crs as ccrs import cartopy.feature as cfeature proj ccrs.NorthPolarStereo(central_longitude45) pc ccrs.PlateCarree() # 扇形参数 lon_min, lon_max 0, 90 lat_min 60 # 在经纬度坐标内构造扇形边界 arc_lons np.linspace(lon_min, lon_max, 200) arc_lats np.full_like(arc_lons, lat_min) lon_seq [lon_min, lon_min] list(arc_lons) [lon_max, lon_max] lat_seq [90, lat_min] list(arc_lats) [lat_min, 90] # 转成NorthPolarStereo投影坐标 verts proj.transform_points(pc, np.array(lon_seq), np.array(lat_seq)) sector_path mpath.Path(verts[:, :2]) fig, ax plt.subplots(figsize(8, 8), subplot_kw{projection: proj}) ax.set_boundary(sector_path) ax.set_extent([-10, 100, 60, 90], crspc) ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linestyle:, linewidth0.5) ax.gridlines(linestyle--, alpha0.6, draw_labelsFalse) plt.show()重点解释一下这段代码的逻辑。lon_seq和lat_seq里我先把起点放在北极点纬度90然后让它沿着lon_min这条经线一路往下到底部纬度lat_min再沿着纬线lat_min从lon_min平滑地走到lon_max最后沿着lon_max经线回到北极点。这样一整圈围出来的就是一个标准的扇面。之所以要采样200个点来做底部圆弧是因为如果只有3个点投影到NorthPolarStereo之后底部就是一条折线视觉上很生硬。采样点越多弧度越圆润。这里选择在经纬度坐标内采样再由transform_points投影到目标坐标系是最不容易出错的方式因为它保证了每两个相邻采样点之间的连线在地理意义上是合理的。3.3 边界裁剪的常见误区用set_boundary的时候我最开始犯过一个错误以为传进去的路径必须在画布投影坐标下定义所以就自己噼里啪啦算了一堆立体几何公式结果不是缝隙就是变形。后来才意识到最简单稳妥的办法就是先在经纬度坐标里构造好路径再交给transform_points转换这个步骤属于常见实践强烈建议直接抄作业。另一个常见误区是set_boundary和set_extent同时使用时的顺序问题。我习惯先set_boundary再set_extent这样既能裁剪扇面又能保证地图的可视范围不会无限延伸到北冰洋以外的区域。如果只set_boundary不set_extent图的范围可能会被自动扩展到很远的地方视觉上扇面会很小。如果只set_extent不set_boundary那就还是一个矩形的裁剪窗口不是扇形。4. 完整实操从底图到数据叠加一步不落4.1 环境准备与基础底图开始跑代码之前先确认环境里装了这几个库python本体之外还需要numpy、matplotlib、cartopy。cartopy的安装偶尔会让你头疼如果使用conda环境会顺畅一些。以我自己的经验到了新电脑上我最常用下面这种方式创建环境conda create -n cartopy_env python3.10 conda activate cartopy_env conda install -c conda-forge cartopy matplotlib numpy装好之后先跑一个最基础的底图脚本确认coastline能正常加载。如果运行时提示“DownloadError”或者“OError: [Errno -2] Name or service not known”多半是海岸线数据下载失败。解决办法有两个一是检查网络二是用cfeature.NaturalEarthFeature的cache目录手动放置下载好的矢量数据。我建议新人在离线或网络不好的环境下直接改用ax.add_feature(cfeature.OCEAN, colorlightblue) ax.add_feature(cfeature.LAND, colorlightgray)这个不依赖在线下载能保证基础底图先跑通后面再慢慢补海岸线细节。4.2 数据图层以海冰密集度为例底图画好了接下来就是往图里叠真正的数据。我最常画的是海冰密集度场数据格式通常是一个二维数组配合二维的经度网格和纬度网格坐标原点在北极附近。叠加数据时一定要记得加transformccrs.PlateCarree()否则数据会被错误地当成投影坐标画到图上结果就是你看到海冰分布全都挤在一个角落。下面是完整示例假设你已经读取了数据lon2d和lat2d是二维经纬度网格sic是海冰密集度数组0到1之间的小数import numpy as np import matplotlib.pyplot as plt import matplotlib.path as mpath import matplotlib.ticker as mticker import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter # --- 扇形参数 --- lat_min 60 lon_min, lon_max -90, 60 # --- 投影设置 --- proj ccrs.NorthPolarStereo(central_longitude-30) pc ccrs.PlateCarree() # --- 构建扇形边界 --- arc_lons np.linspace(lon_min, lon_max, 200) arc_lats np.full_like(arc_lons, lat_min) lon_seq [lon_min, lon_min] list(arc_lons) [lon_max, lon_max] lat_seq [90, lat_min] list(arc_lats) [lat_min, 90] verts proj.transform_points(pc, np.array(lon_seq), np.array(lat_seq)) sector_path mpath.Path(verts[:, :2]) # --- 画布 --- fig plt.figure(figsize(9, 9)) ax fig.add_subplot(1, 1, 1, projectionproj) ax.set_boundary(sector_path) ax.set_extent([lon_min - 20, lon_max 20, lat_min - 5, 90], crspc) # --- 底图要素 --- ax.add_feature(cfeature.LAND, colorlightgray, zorder2) ax.add_feature(cfeature.OCEAN, colorwhite, zorder1) ax.coastlines(resolution50m, linewidth0.7, colorblack, zorder3) # --- 模拟一份海冰数据实际使用时替换成自己的数组 --- lons_grid np.linspace(lon_min, lon_max, 120) lats_grid np.linspace(lat_min, 88, 100) lon2d, lat2d np.meshgrid(lons_grid, lats_grid) sic np.clip(np.cos(np.deg2rad(lat2d - 55)) * 0.8 0.1, 0, 1) # --- 画数据 --- mesh ax.pcolormesh(lon2d, lat2d, sic, transformpc, cmapSpectral_r, shadingauto, zorder4) cbar fig.colorbar(mesh, axax, shrink0.7, pad0.05) cbar.set_label(Sea Ice Concentration) # --- 网格线 --- gl ax.gridlines(crspc, linewidth0.5, colorgray, alpha0.6, linestyle--, draw_labelsTrue, xlocsrange(-180, 181, 30), ylocsrange(0, 91, 10)) gl.top_labels False gl.right_labels False gl.xformatter LongitudeFormatter() gl.yformatter LatitudeFormatter() plt.show()这段代码里我做了几件比较重要的事一是把扇形边界通过set_boundary固定住让图层无论是底图还是数据都会被裁到扇形之内二是用pcolormesh叠加数据并用Spectral_r配色展示海冰密集度暖色代表密集度高冷色代表密集度低三是用gridlines的xlocs和ylocs手动控制了经纬度网格线的密度避免线条太密糊成一团。这里需要提醒一下pcolormesh的shadingauto参数是matplotlib比较新的特性能自动判断是网格还是散点旧版本可能需要改成shadingnearest或gouraud。如果你的matplotlib版本比较老遇到“shading auto”不兼容的报错直接改掉就行。4.3 网格线、标注与出图设置网格线的设置看起来是小事但在极区图上特别容易翻车。如果你用默认的gridlines()去画北极扇形图经常会出现经纬度标注挤在扇面周围的奇怪位置上甚至有些标签飘到边界外面。比较好的做法是像我上面那样给gridlines指定crsccrs.PlateCarree()并手动指定经纬度刻度位置然后用LongitudeFormatter和LatitudeFormatter来格式化标签。还要注意top_labels和right_labels在极区图里通常设为False不然图上面和右面的标注会和左边的撞在一起。我自己的习惯是只显示左边纬度和下边经度不但清爽也符合大多数论文图表的阅读习惯。最后是出图设置。扇形图建议用正方形画布也就是figsize(8, 8)或者(9, 9)因为扇形本身在NorthPolarStereo投影下就是一个以中心为顶点的圆饼局部正方形画布能最大程度利用空间。保存时用dpi300格式用png或者pdf都行。如果发给导师或者写报告用我一般保存pdf——矢量格式缩放不失真投到论文里也够清晰fig.savefig(arctic_sector.png, dpi300, bbox_inchestight) fig.savefig(arctic_sector.pdf, bbox_inchestight)5. 高频报错与排错经验能帮你省两小时5.1 问题速查表下面这几类问题是我在画北极扇形图时遇到过也在网上帮别人排查过的。这里整理成表格遇到报错可以对号入座。问题现象可能原因解决思路add_feature时报错或海岸线不出来自然地球矢量数据下载失败先检查网络或者改用OCEAN/LAND基础色块画出的图是一片空白没有数据数据范围不在投影可见范围内检查set_extent范围确认数据和extent是否有交集set_boundary后经纬网格线消失裁剪边界格式不对或transform缺失确认path是投影坐标且不是DataFrame数组set_extent报错“ValueError: CRS mismatch”范围参数没有指定crs写成ax.set_extent([...], crsccrs.PlateCarree())图形出现一条从顶部穿过的斜线多边形路径没有闭合或者采样点太少检查lon_seq和lat_seq首尾是否回到北极点图上横坐标刻度太密集标签挤在一起gridlines自动生成过多刻度手动指定xlocs、ylocs或者用MaxNLocator控制5.2 经度跨边界-180与180之间那道裂缝这个坑非常容易遇到。如果你要画的区域跨过了日期变更线比如从东经170度到西经160度直接用一次linspace得到的是一个横跨整个太平洋的大弧形画出来完全不是你要的区域。cartopy在处理这种跨0度或跨180度的数据时尤其需要小心。解决办法有两个。推荐的做法是把经度范围拆成两段分别构造后再合并。比如要画东经170到西经160也就是在坐标上是从170到200或从-190到-160那你可以统一用“-190到-160”的连续区间来表达lon_min-190, lon_max-160。NorthPolarStereo投影本身不关心经度值是正还是负只要数据坐标与边界线段能不解纠缠就行。如果你的数据本身经度是0到360的格式更好办直接统一到“0到360”范围内计算。我的经验是在做扇形区域图之前先把所有经度统一成同一种习惯要么都是-180到180要么都是0到360。千万别混用不然边界路径和数据图层各画各的出了图就再也对不上了。5.3 海岸线、数据错位与画图横坐标过密的通用解法错位问题多半是transform没写对。比如海岸线是cartopy内置的它知道自己的坐标系统但你自己的海冰数据一定要在pcolormesh里指定transform。如果数据出处不明、网格本身就是投影坐标那就先把数据插值到经纬度网格再画不要偷懒直接塞进pcolormesh。至于热搜里经常有人问“python画图横坐标太密集”这个在cartopy里的表现就是经纬度网格线标签叠在一起。我自己的做法是给gridlines传入xlocs和ylocs手动指定要显示哪些刻度。比如在北极扇形图中我通常让经度每隔30度、纬度每隔10度显示一个标签清晰又不拥挤。如果是普通的matplotlib图则可以使用matplotlib.ticker的MaxNLocator或MultipleLocator控制刻度数量。核心思路是一样的别让绘图库自动决定刻度密度自己动手定好间隔一劳永逸。如果再遇到那些莫名其妙的报错我建议把代码逐步注释掉从底图开始一点点加要素。哪一步加了之后图形开始异常问题就锁定在哪一步。这个方法虽然笨但非常有效。我自己做了这么多次极区图之后最大的体会是扇形图看着复杂本质上并不是某种神奇操作而是把“投影”、“边界路径”、“数据transform”这三个链路理清楚剩下的只是往图里一层一层叠要素。只要你对“经纬度坐标出发、由投影转换到画布坐标”这条主线有把握不管换成北半球还是南半球还是任意一块“切角区域”都能照葫芦画瓢地做出来。如果你正准备画一块扇形区域但又不知道怎么处理边界可以直接把我上面的扇形路径构造函数摘出去改成你要的经度范围和底边纬度就行。这样省下来的时间拿去调配色和排版它不香吗
