尧图网页设计 
1. 扇形北极图是怎样一种需求先弄清楚你要画的是什么做气象、海洋或者极地研究的朋友应该都有过类似的经历手里的数据覆盖了整个北半球甚至全球但真正关心的区域就在北纬60度以北那一圈。这时候你在屏幕上拖来拖去想把地图范围限定在北极附近结果用常规的矩形框截出来要么是一块被拉伸到几乎认不出的变形图要么就是一大片空白海域看着特别浪费。“扇形区域图”这个名字听起来有点玄实际上它就是专门为极地/near-polar区域设计的一种可视化方案。把北极点放在图上某个中心位置纬度从北极点向外一圈一圈扩展经度从一侧边界到另一侧边界像扇子一样张开。效果就说把北冰洋周边一圈陆地北美北部、格陵兰、北欧、西伯利亚完整包裹进视野没有任何几何失真带来的难看压迫感。我最早接触这个需求是帮一个做海冰预报的课题组画“北极海冰密集度分布图”。他们要求图的长宽比大约1:1.2北极点不能像默认设置那样偏到左边或者被裁掉整个场景要像雷达屏幕一样把极区“罩”进来。用Matplotlib自带的Basemap试了一版虽然能出图但速度慢、样式老旧后来切到Cartopy之后才真正顺畅起来。这也是为什么今天的文章核心是基于Cartopy来展开。简单说这篇文章适合谁看你手上有北半球极区相关的数据海冰、气温、风场、位势高度都可以想用Python画出一张能直接放进论文或者汇报PPT里的北极扇形区域图同时还想搞清楚背后的投影、范围和裁剪逻辑而不是每次网上搜一段代码碰运气。2. 投影坐标系是核心为什么北极图必须用特殊投影2.1 Cartopy里和北极投影有关的那几个类Cartopy是整个可视化链条里最关键的一环它不负责画数据本身画数据依然是Matplotlib的事Cartopy做的事情是把地理坐标比如经纬度转换到画布上的平面坐标。这一步要是错了后面的海岸线、等值线、填色图就全部对不上位置。和北极扇形图直接相关的投影类主要有这几个投影类投影中心典型效果ccrs.Orthographic可以指定中心点坐标从太空中某个点看地球边缘明显弯曲能看到整个半球但靠近边缘压缩严重ccrs.NorthPolarStereo()北极点以北极点为中心的极射赤面投影极区形态横向展开适合画扇形区域ccrs.Stereographic(central_latitude90)等价于NorthPolarStereo本质是同一个东西ccrs.LambertConformal(central_latitude70)可以指定中高纬度变形小但全北极一圈覆盖起来边缘会有些奇怪我实测下来画“北极部分区域”的扇形图NorthPolarStereo是最恰当的默认选择。它的纬线结构是从北极点向外一圈圈放大经线则从一点向外辐射天然就是扇形的骨架。至于为什么不用Orthographic因为Orthographic的中心视角决定了一旦数据范围超过中心点约70度图的外缘就开始“滑落”到球体背面画面会出现一团模糊的边缘而且靠近边缘的投影网格极度压缩不适合放入定量信息。2.2 你不需要从零学习投影数学但这三个坐标系概念必须分清很多人第一次在Cartopy里画极区图报错“Data has no positive latitude”或者图形扭曲得离谱问题几乎出在坐标系混淆上。Cartopy里至少有三个层面的坐标系统我按使用频率从高到低排列数据坐标系你的数据本身所在的空间绝大多数是PlateCarree就是经纬度数据别管它叫“地理坐标”还是什么在Cartopy里统一用ccrs.PlateCarree()表示。地图投影坐标系就是你选定的NorthPolarStereo也叫目标投影。轴对象在创建时通过projection参数指定。显示坐标系Matplotlib的显示像素坐标一般不需要直接操作。记住一个硬规则绘图时凡是带经纬度的数据都要在后面通过transform参数声明它的原始坐标系。比如ax.pcolormesh(lon, lat, data, transformccrs.PlateCarree())这个transform参数其实就是告诉Cartopy“喂我给你的这些lon/lat数据是用经纬度表示的默认的PlateCarree你给我按这个解释。” 如果不写或者写错Cartopy就会默认你的数据已经处于目标投影下结果就是整片数据被挪到不知名角落。2.3 为什么“扇形”的本质是裁剪而不是投影很多第一次接触的人会以为只要选了NorthPolarStereo画出来的图天然就是扇形。实际上完全不是。NorthPolarStereo只是改变坐标计算方式如果你不限制范围Cartopy默认会把整个北极点周边可见的圆盘范围都显示出来也就是北半球全部出来的是一张“圆饼”不是“扇形”。“扇形”还有一层特殊含义它有一个起始经度和一个终止经度纬度通常是北纬某个下限比如60度。换句话说扇形区域图等于“极射赤面投影整盘显示”“按角度范围裁剪出一块”。这个裁剪动作才是文章标题里“扇形”的关键。再具体一点假如我想画的是“北纬60度以北、东经90度到西经90度覆盖的那一瓣”那就意味着在地图投影的球面上我们要把范围框成一个从北极点伸出去的楔形并只显示这部分。手动裁剪本身很麻烦幸好Cartopy检查了自定义边界boundary通过重写set_boundary可以精准切成扇形。3. 实操第一版让扇形边界真正落在图上3.1 全过程代码先跑通再理解我直接给出一段能跑的代码。这段代码是我平时做极区可视化时的常用底座你大概率可以直接改一改参数就用上不用从零写import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import numpy as np # 定义扇形区域的范围 lon_min, lon_max 60, 120 # 感兴趣的经度窗口比如东西伯利亚海 lat_min 60 # 起始纬度北纬60度以北 # 创建极射赤面投影坐标 proj ccrs.NorthPolarStereo(central_longitude0) fig plt.figure(figsize(6, 8), dpi150) ax plt.axes(projectionproj) # 最关键的一步把地图轴限制在“北极点向外、lat_min为半径”的圆形区域内 # 原理通过构造一个在投影坐标下的圆多边形然后set_boundary裁剪 import shapely.geometry as sgeom # 在经纬度坐标下生成一个覆盖[lon_min, lon_max] × [lat_min, 90]的方框 # 但我方框的南北极附近做法比较特殊先在整个极射投影下划定圆形范围 # 更稳妥先设定ax的显示范围再裁剪成扇形等一下上面这段只是框架。真正的扇形裁剪有一个比较经典的三步套路第一步设定经纬度范围后先用ax.set_extent把地图大致框到北半球高纬区域。但set_extent配合NorthPolarStereo时有个坑——它希望你填写的范围是该投影下的边界框如果在高纬度地区填一个普通的经纬度矩形出来的图会有大块空白或者被拉扁。第二步构造一个matplotlib.path.Path对象路径在数据坐标经纬度下是扇形轮廓即纬度从lat_min到90度、经度从lon_min到lon_max边沿在lat_min处有一条弧线。我通常的做法是先用np.linspace把扇形轮廓的顶点全部生成出来然后转换到投影坐标再用ax.set_boundary应用这个轮廓。第三步把生成的轮廓路径和Axes的set_boundary结合起来。这一步是真正的裁剪魔法。以下是能直接运行的完整代码import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import numpy as np import matplotlib.path as mpath # 1. 定义扇形范围 lon_min, lon_max 60, 120 lat_min 60 # 2. 投影 proj ccrs.NorthPolarStereo(central_longitude0) fig plt.figure(figsize(6, 8), dpi150) ax plt.axes(projectionproj) # 3. 生成扇形轮廓顶点经纬度坐标 # 从lon_min到lon_max沿lat_min画弧然后沿经线回到北极点 theta np.linspace(lon_min, lon_max, 100) verts [] # 先沿纬线圈lat_min从lon_min到lon_max for lon in theta: verts.append((lon, lat_min)) # 再沿经线lon_max从lat_min到90 for lat in np.linspace(lat_min, 90, 30): verts.append((lon_max, lat)) # 再沿经线lon_min从90回到lat_min for lon in theta[::-1]: verts.append((lon, 90)) # 这样会形成闭合轮廓吗注意我们其实需要从90度到lat_min沿lon_min补充一段其实这个顶点顺序有点乱。更稳妥的做法是用shapely生成一个多边形对象然后转成Path但CarTop边界又要求必须是matplotlib.path.Path对象。下面这个办法最直接、也最不容易出错def make_sector_path(lon_min, lon_max, lat_min, n100): 生成从lat_min到北极点、经度范围为[lon_min, lon_max]的扇形路径。 返回的Path在经纬度坐标系下。 # 沿着纬线圈从lon_min到lon_max lons1 np.linspace(lon_min, lon_max, n) lats1 np.full_like(lons1, lat_min) # 沿着lon_max经线从lat_min到90 lats2 np.linspace(lat_min, 90, n) lons2 np.full_like(lats2, lon_max) # 沿着纬线圈接近北极点从lon_max回到lon_min实际收缩成一点 lons3 np.linspace(lon_max, lon_min, n) lats3 np.full_like(lons3, 90.0) # 沿着lon_min经线从90回到lat_min lats4 np.linspace(90, lat_min, n) lons4 np.full_like(lats4, lon_min) lons np.concatenate([lons1, lons2, lons3, lons4]) lats np.concatenate([lats1, lats2, lats3, lats4]) # 用Path构造 return mpath.Path(np.column_stack([lons, lats]))注意上面这个Path是在经纬度坐标系近似PlateCarree下构建的所以创建Axes之后还需要把它转换到投影坐标系path make_sector_path(lon_min, lon_max, lat_min) # 在投影坐标下重新生成这个Path proj_path proj.project_geometry(sgeom.LineString(path.vertices))等一下project_geometry需要一个shapely geometry对象而LineString需要至少两个点并且要求所有点的经度不能产生跨180度的问题。如果你的扇形跨了日期变更线比如从160度到-160度直接构造LineString会得到一条横穿全球的线这时候需要把多边形拆成两段分别处理。这是后话我先给一个不跨日期线的一般用法import shapely.geometry as sgeom def sector_polygon(lon_min, lon_max, lat_min): 返回经纬度坐标下的扇形shapely Polygon n 100 lons1 np.linspace(lon_min, lon_max, n) lats1 np.full_like(lons1, lat_min) lons2 np.linspace(lon_max, lon_min, n) lats2 np.full_like(lons2, 90.0) coords list(zip(lons1, lats1)) list(zip(lons2, lats2)) return sgeom.Polygon(coords) polygon sector_polygon(lon_min, lon_max, lat_min) proj_polygon proj.project_geometry(polygon) # 投影到目标坐标系 boundary_path mpath.Path(np.asarray(proj_polygon.exterior.coords))这还不够因为proj.project_geometry返回的坐标已经是投影坐标米而ax.set_boundary需要在数据坐标即投影坐标下应用这样刚好对得上。最后应用边界ax.set_boundary(boundary_path, transformproj)等等这种用法有个版本差异问题旧版Cartopyset_boundary不强制transform参数新版如果传了transform而数据坐标不是投影数据会导致奇怪偏移。实际我测试下来最稳妥的写法是直接ax.set_boundary(boundary_path)set_boundary默认就采用当前Axes的投影坐标系来解读路径坐标所以不需要再传transform。这一点很多人会踩坑特意标注一下。加上地图要素ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linestyle:, linewidth0.6) ax.gridlines(draw_labelsTrue, linewidth0.5, alpha0.5)最终出来的效果北极点附近一个完整的扇形边界是笔直的经线侧边和一段弧线底边海岸线、国界和网格线都被裁剪在这个扇形里。这个就是标题要的“扇形区域图”。3.2 细节打磨为什么我的图右上角会出现“破洞”或者“飞线”用上面的方法跑通之后很多人会遇到一个老问题扇形以外的海岸线、网格线还是会穿过边界“漏”出来一部分像是什么东西从扇形边缘飞出去。这是因为set_boundary只限制了Axes的绘制范围但Gridline、Coastline这些Feature的绘制方式在某些版本里并不会严格被Boundary裁剪。最简单的解决方案关掉draw_labels或者手动控制网格标签的位置。我在实际项目中更喜欢的做法是关闭自动label叠加手动标注gl ax.gridlines(draw_labelsFalse, linewidth0.5, alpha0.5, linestyle--)然后手工加上想要的经纬度标注特别是扇形边缘的经线、底边的纬线。具体做法在后面的章节里再说。3.3 数据填充把格点数据放到扇形图里画边界只是第一步实际项目里真正核心的是往扇形里填你的物理场。比如我有一个T2M2米气温的全球格点数据lon和lat都是1度间隔的二维数组现在只想显示北纬60到90度、东经60到120度的区域。第一步把数据坐标转换到投影坐标之前先裁剪数据本身的范围以减小绘制压力避免NorthPolarStereo投影时把南北极附近已经重复的格点重复投影。第二步绘制时一定要带transformccrs.PlateCarree()我见过太多人因为漏了这个参数导致填色错了几万公里lon2d, lat2d np.meshgrid(lons, lats) # 数据裁剪 mask (lat2d lat_min) (lon2d lon_min) (lon2d lon_max) data_masked np.ma.masked_where(~mask, data) cf ax.contourf(lon2d, lat2d, data_masked, levels20, cmapRdBu_r, transformccrs.PlateCarree()) fig.colorbar(cf, orientationhorizontal, pad0.05, aspect40)这里注意虽然是北纬60度以上区域但如果数据分辨率很高比如0.1度一次性投影全部格点会非常慢。一个常用的技巧是只数据裁剪后用block方式降低分辨率绘制或者直接scatter少量站点。等值线填色其实对分辨率敏感度不高可以先做1度格点再contourf。对于海冰、SST等有明显边缘的数据contourf比pcolormesh更平滑并且可以直接加extendboth控制溢出色标cf ax.contourf(lon2d, lat2d, data_masked, levelsnp.linspace(0, 100, 21), cmapBlues, extendboth, transformccrs.PlateCarree())4. 地图要素与网格标签让扇形图“像”一张学术地图4.1 海岸线和边界线亮度和宽度的控制画扇形图时地图要素需要比常规地图更谨慎因为极区地形复杂海岸线支离破碎默认设置往往要么太粗影响数据读图要么太淡导致根本看不清。我常用的一套基线配置ax.add_feature(cfeature.LAND, color#e8e0d0, zorder1) ax.add_feature(cfeature.OCEAN, color#b0d4e8, zorder1) ax.add_feature(cfeature.COASTLINE, edgecolorblack, linewidth0.8, zorder2) ax.add_feature(cfeature.BORDERS, linestyle:, edgecolorgray, linewidth0.5, zorder2)几个重点提醒LAND/OCEAN的颜色别用纯白纯蓝因为后面一旦叠加半透明数据比如透明度0.7的填色底图颜色会透出来形成奇怪的混合色。使用浅的暖色和冷色能保证数据层的辨识度。如果只是做灰度打印直接用color0.7作为陆地色更稳妥。zorder务必设置尤其当数据填色zorder默认是1如果底图要素也是1可能出现海岸线被数据盖住的尴尬局面。4.2 经纬网与标签手动标注边缘经度自动gridlines(draw_labelsTrue)在极坐标投影下经常出现标签叠成一团或者出现在扇形外部的情况。我的经验是全部手动控制。先关闭自动标签只画网格gl ax.gridlines(crsccrs.PlateCarree(), draw_labelsFalse, linewidth0.4, colorgray, alpha0.6, linestyle--)然后标注扇形边界位置。比如扇形左边经度是60E右边经度是120E底边纬度是60N。我可以在投影坐标系下的对应位置通过ax.text手工写# 在底边中点标注纬度 ax.text((lon_minlon_max)/2, lat_min, 60°N, transformccrs.PlateCarree(), hacenter, vabottom)也可以把这些标注放在扇形左下和右下角ax.text(lon_min, lat_min2, 60°E, transformccrs.PlateCarree(), hacenter, vacenter) ax.text(lon_max, lat_min2, 120°E, transformccrs.PlateCarree(), hacenter, vacenter)这里有个小坑ax.text的transformccrs.PlateCarree()表示我指定的坐标是经纬度它会自动换算到投影坐标。你需要确保文本没被边界裁剪掉。默认情况下在Axes边界外的text是不会显示的所以文本位置要稍微向扇形内部偏移一点比如lat_min2。4.3 北极点标记与极点周边“死区”另一个常见需求是在图上标出北极点。北极点在NorthPolarStereo投影下就是圆心位置直接用ax.plot(0, 90, o, markersize6, colorred, transformccrs.PlateCarree()) ax.text(0, 89.5, NP, transformccrs.PlateCarree(), hacenter)但要注意如果你的数据网格里有准确的90N格点很多全球数据在90N处只有一个点投影后是一个极小的区域放大后会看到一个孤立的色块或者被渲染成一条线。这是正常现象。处理办法是在绘制之前把数据中的90N那一行极其靠近北极点的数据mask掉或者用contourf绘制使它成为扇形最顶端的一个平滑色块而不是一个点。5. 极端边界情况跨180度经线和半球范围5.1 经度窗口跨日期变更线怎么办很多时候你关心的区域并不老老实实待在0到360度的范围里。比如研究楚科奇海和白令海需要显示从东经160度到西经150度这个窗口恰好跨了180度经线。直接用前面的np.linspace(lon_min, lon_max)会得到一段横穿太平洋全球的长线导致多边形异常。解决办法有两种第一种将经度窗口统一写为0到360度。比如东经160度记为160西经150度记为210于是窗口是[160, 210]这样就不跨180度了。相应地进行经纬度范围判断时也需要把数据经度做统一变换lon360 np.where(lon 0, lon 360, lon) mask (lat lat_min) (lon360 160) (lon360 210)第二种拆分成两个窗口比如160到180、-180到-150分别画但这样处理起来繁琐且容易在接口处出现拼接缝隙。我推荐第一种操作简单性能损耗最低。5.2 宽扇形从90E到120W如果经度范围跨度太大比如150度以上扇形明显变宽整个图会变得非常扁。这时需要调整画布figsize的比例。比如经度跨度为100度时用figsize(8, 5)经度跨度为40度时用figsize(4, 6)。比例失控会导致经线看起来倾斜严重不美观。一个经验公式扇形“展开角”越大画布宽度比例应相应增加但别超过21否则视觉上整个北极被拉得不成形。5.3 当投影中心不在本初子午线默认NorthPolarStereo(central_longitude0)把0度经线放在正下方。如果你希望扇形区域看起来“正”。比如你关注的区域是180度经线附近白令海那最好设置central_longitude180这样180度经线会位于扇形正下方而不是卷曲到侧面proj ccrs.NorthPolarStereo(central_longitude180)实测效果变化很大尤其是跨大经度范围时。推荐做法是先整体看数据分布的质心经度把central_longitude设为那个值附近扇形图会显得更“正”。6. 进阶把扇形图嵌入更大的极区场景6.1 与整个北半球图拼接有些论文需要在一张图里同时展示“全北极范围”和“局部放大扇形”。我用两种方式实现过子图方式用fig.add_subplot(1, 2, 1)画全北极fig.add_subplot(1, 2, 2)画局部扇形两个Axes使用同一个投影ccrs.NorthPolarStereo但第二个Axes的边界范围不同。插图方式在扇形图内部添加一个小坐标轴用于指路类似“地图中的地图”。第二种方式在论文中更常见因为可以直观地告诉读者扇形区域在北极的位置。实现代码ax_sub fig.add_axes([0.72, 0.75, 0.2, 0.2], projectionccrs.NorthPolarStereo()) ax_sub.add_feature(cfeature.COASTLINE, linewidth0.4) ax_sub.set_extent([-180, 180, 65, 90], crsccrs.PlateCarree()) # 画一个表示扇形边界的矩形 ax_sub.add_patch(plt.Rectangle((lon_min, lat_min), lon_max-lon_min, 90-lat_min, facecolorred, alpha0.3, transformccrs.PlateCarree()))注意这个add_patch里的transformccrs.PlateCarree()是必需的否则矩形会被放在投影坐标里位置完全错乱。6.2 高分辨率背景与在线瓦片如果觉得Cartopy自带的GSHHS海岸线不够精细可以叠加Cartopy的img_tiles功能加载在线底图。但从实际体验看在一些网络不稳的环境中在线瓦片加载很慢而且极区瓦片覆盖并不完整。我更推荐本地自然地球Natural Earth矢量数据或分辨率更高的GADM国界数据。Cartopy自带的是NaturalEarthFeatureax.add_feature(cfeature.NaturalEarthFeature(physical, land, 10m, edgecolorblack, facecolor#e8e0d0))注意下载10m分辨率数据需要联网如果离线环境默认用50m即可。6.3 加入站点散点和轨迹线网格数据之后站点和航线散点是最常见的叠加元素。站点绘制时同样必须指定transformccrs.PlateCarree()。比如ax.scatter(site_lon, site_lat, s30, cred, edgecolorsblack, linewidths0.5, transformccrs.PlateCarree(), zorder5)航线直接把经纬度序列连线ax.plot(route_lons, route_lats, colorblue, linewidth1.5, transformccrs.PlateCarree(), zorder4)这里有一个极区特有的大坑如果你的轨迹线有跨180度的经度跳变比如从179到-179在投影坐标下会画出一条横穿整张图的直线。处理办法与5.1相同先统一成0~360或者拆分线段。7. 常见报错与排查链路7.1ValueError: Cannot handle non-finite values或Invalid transform sphere这类报错通常有两种原因一是数据中出现了NaN二是经纬度网格中有不符合PlateCarree的坐标如纬度超过90或者经度超过360。排查路径先检查数据np.isfinite(data).all()如果False进行差值或剔除。检查经纬度范围纬度必须在-90到90之间经度在-180到360之间。检查是否有坐标数组与数据shape不一致遮罩矩阵写错导致坐标轴对齐出错。7.2 扇形边界裁剪后内侧出现白色缺口经常看到扇形底边有一段白色空隙海岸线但是海陆底色没有填进去。这多半是LAND要素的zorder低于扇形边界背景层。先尝试把ax.add_feature(cfeature.LAND, ...)放在set_boundary之后执行如果不行就在set_boundary之后用ax.add_patch补一层扇形内部背景色。7.3 输出PDF时扇形边缘出现细白线在PDF输出中set_boundary裁剪的路径边缘经常会出现一条极细的空白线这是渲染器的抗锯齿问题。给Axes加一个宽度足够的patch边框可以压制ax.patch.set_linewidth(1) ax.patch.set_edgecolor(black) # 或者和背景色一致7.4 Cartopy和Shapely版本冲突project_geometry在不同版本上的行为略有差异。旧版Cartopy0.18左右在投影某些跨越180度的多边形时可能会生成异常形状。遇到这种情况先升级Cartopy到最新版本0.21再把数据经度统一到0~360。我已经多次靠这个组合解决“多边形突然扭曲成一个圆”的问题。8. 我的选型判断与一些私人经验最后聊几句选型层面的体会。用Cartopy画北极扇形区域图核心工作量其实不在“画”本身而在“把投影和裁剪彻底想明白”。很多人绕不开Basemap因为Basemap的API更直观设置llcrnrlat, urcrnrlat就能直接限定范围但它早在2020年就停止维护了在新版Matplotlib上频繁报np.float相关的兼容性错误。Cartopy虽然初期接触时总觉得“transform这东西真麻烦”可一旦形成心智模型后面画任何投影图都是一通百通。我个人在实际项目中最终沉淀下的固定配置组合是这样的投影优先ccrs.NorthPolarStereo(central_longitude需要对齐的经度)扇形路径用shapely.geometry.Polygon构造数据坐标用0~360的经度体系避免跨180度问题数据绘制时所有经纬度散点、填色、等值线全部带transformccrs.PlateCarree()网格标签用手动ax.text标注不用gridlines(draw_labelsTrue)省去极区标签乱飞的烦恼;输出时优先fig.savefig(xxx.pdf, bbox_inchestight)矢量格式能避免放大后锯齿。这两个“手动”替代“自动”的习惯帮我省下过无数用来调样式的时间。反正极区图的标签就那么几根手写几分钟自动调整反而常常半小时还没看顺眼。如果你只是临时画一张给自己看的数据快照那上面章节里的代码足够了但如果你要做得更细比如论文出版图、甚至要出动画纬度带选择、经度中心对齐、标注位置、色彩方案这些细节是值得一调再调的。极区图最怕一带而过它比中低纬度图更容易因为投影变形引发误读。做图的人和看图的人得在同一个坐标系里对话你的图才真正站得住。 用Cartopy绘制北极扇形区域图：投影、裁剪与可视化实战
1. 扇形北极图是怎样一种需求先弄清楚你要画的是什么做气象、海洋或者极地研究的朋友应该都有过类似的经历手里的数据覆盖了整个北半球甚至全球但真正关心的区域就在北纬60度以北那一圈。这时候你在屏幕上拖来拖去想把地图范围限定在北极附近结果用常规的矩形框截出来要么是一块被拉伸到几乎认不出的变形图要么就是一大片空白海域看着特别浪费。“扇形区域图”这个名字听起来有点玄实际上它就是专门为极地/near-polar区域设计的一种可视化方案。把北极点放在图上某个中心位置纬度从北极点向外一圈一圈扩展经度从一侧边界到另一侧边界像扇子一样张开。效果就说把北冰洋周边一圈陆地北美北部、格陵兰、北欧、西伯利亚完整包裹进视野没有任何几何失真带来的难看压迫感。我最早接触这个需求是帮一个做海冰预报的课题组画“北极海冰密集度分布图”。他们要求图的长宽比大约1:1.2北极点不能像默认设置那样偏到左边或者被裁掉整个场景要像雷达屏幕一样把极区“罩”进来。用Matplotlib自带的Basemap试了一版虽然能出图但速度慢、样式老旧后来切到Cartopy之后才真正顺畅起来。这也是为什么今天的文章核心是基于Cartopy来展开。简单说这篇文章适合谁看你手上有北半球极区相关的数据海冰、气温、风场、位势高度都可以想用Python画出一张能直接放进论文或者汇报PPT里的北极扇形区域图同时还想搞清楚背后的投影、范围和裁剪逻辑而不是每次网上搜一段代码碰运气。2. 投影坐标系是核心为什么北极图必须用特殊投影2.1 Cartopy里和北极投影有关的那几个类Cartopy是整个可视化链条里最关键的一环它不负责画数据本身画数据依然是Matplotlib的事Cartopy做的事情是把地理坐标比如经纬度转换到画布上的平面坐标。这一步要是错了后面的海岸线、等值线、填色图就全部对不上位置。和北极扇形图直接相关的投影类主要有这几个投影类投影中心典型效果ccrs.Orthographic可以指定中心点坐标从太空中某个点看地球边缘明显弯曲能看到整个半球但靠近边缘压缩严重ccrs.NorthPolarStereo()北极点以北极点为中心的极射赤面投影极区形态横向展开适合画扇形区域ccrs.Stereographic(central_latitude90)等价于NorthPolarStereo本质是同一个东西ccrs.LambertConformal(central_latitude70)可以指定中高纬度变形小但全北极一圈覆盖起来边缘会有些奇怪我实测下来画“北极部分区域”的扇形图NorthPolarStereo是最恰当的默认选择。它的纬线结构是从北极点向外一圈圈放大经线则从一点向外辐射天然就是扇形的骨架。至于为什么不用Orthographic因为Orthographic的中心视角决定了一旦数据范围超过中心点约70度图的外缘就开始“滑落”到球体背面画面会出现一团模糊的边缘而且靠近边缘的投影网格极度压缩不适合放入定量信息。2.2 你不需要从零学习投影数学但这三个坐标系概念必须分清很多人第一次在Cartopy里画极区图报错“Data has no positive latitude”或者图形扭曲得离谱问题几乎出在坐标系混淆上。Cartopy里至少有三个层面的坐标系统我按使用频率从高到低排列数据坐标系你的数据本身所在的空间绝大多数是PlateCarree就是经纬度数据别管它叫“地理坐标”还是什么在Cartopy里统一用ccrs.PlateCarree()表示。地图投影坐标系就是你选定的NorthPolarStereo也叫目标投影。轴对象在创建时通过projection参数指定。显示坐标系Matplotlib的显示像素坐标一般不需要直接操作。记住一个硬规则绘图时凡是带经纬度的数据都要在后面通过transform参数声明它的原始坐标系。比如ax.pcolormesh(lon, lat, data, transformccrs.PlateCarree())这个transform参数其实就是告诉Cartopy“喂我给你的这些lon/lat数据是用经纬度表示的默认的PlateCarree你给我按这个解释。” 如果不写或者写错Cartopy就会默认你的数据已经处于目标投影下结果就是整片数据被挪到不知名角落。2.3 为什么“扇形”的本质是裁剪而不是投影很多第一次接触的人会以为只要选了NorthPolarStereo画出来的图天然就是扇形。实际上完全不是。NorthPolarStereo只是改变坐标计算方式如果你不限制范围Cartopy默认会把整个北极点周边可见的圆盘范围都显示出来也就是北半球全部出来的是一张“圆饼”不是“扇形”。“扇形”还有一层特殊含义它有一个起始经度和一个终止经度纬度通常是北纬某个下限比如60度。换句话说扇形区域图等于“极射赤面投影整盘显示”“按角度范围裁剪出一块”。这个裁剪动作才是文章标题里“扇形”的关键。再具体一点假如我想画的是“北纬60度以北、东经90度到西经90度覆盖的那一瓣”那就意味着在地图投影的球面上我们要把范围框成一个从北极点伸出去的楔形并只显示这部分。手动裁剪本身很麻烦幸好Cartopy检查了自定义边界boundary通过重写set_boundary可以精准切成扇形。3. 实操第一版让扇形边界真正落在图上3.1 全过程代码先跑通再理解我直接给出一段能跑的代码。这段代码是我平时做极区可视化时的常用底座你大概率可以直接改一改参数就用上不用从零写import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import numpy as np # 定义扇形区域的范围 lon_min, lon_max 60, 120 # 感兴趣的经度窗口比如东西伯利亚海 lat_min 60 # 起始纬度北纬60度以北 # 创建极射赤面投影坐标 proj ccrs.NorthPolarStereo(central_longitude0) fig plt.figure(figsize(6, 8), dpi150) ax plt.axes(projectionproj) # 最关键的一步把地图轴限制在“北极点向外、lat_min为半径”的圆形区域内 # 原理通过构造一个在投影坐标下的圆多边形然后set_boundary裁剪 import shapely.geometry as sgeom # 在经纬度坐标下生成一个覆盖[lon_min, lon_max] × [lat_min, 90]的方框 # 但我方框的南北极附近做法比较特殊先在整个极射投影下划定圆形范围 # 更稳妥先设定ax的显示范围再裁剪成扇形等一下上面这段只是框架。真正的扇形裁剪有一个比较经典的三步套路第一步设定经纬度范围后先用ax.set_extent把地图大致框到北半球高纬区域。但set_extent配合NorthPolarStereo时有个坑——它希望你填写的范围是该投影下的边界框如果在高纬度地区填一个普通的经纬度矩形出来的图会有大块空白或者被拉扁。第二步构造一个matplotlib.path.Path对象路径在数据坐标经纬度下是扇形轮廓即纬度从lat_min到90度、经度从lon_min到lon_max边沿在lat_min处有一条弧线。我通常的做法是先用np.linspace把扇形轮廓的顶点全部生成出来然后转换到投影坐标再用ax.set_boundary应用这个轮廓。第三步把生成的轮廓路径和Axes的set_boundary结合起来。这一步是真正的裁剪魔法。以下是能直接运行的完整代码import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import numpy as np import matplotlib.path as mpath # 1. 定义扇形范围 lon_min, lon_max 60, 120 lat_min 60 # 2. 投影 proj ccrs.NorthPolarStereo(central_longitude0) fig plt.figure(figsize(6, 8), dpi150) ax plt.axes(projectionproj) # 3. 生成扇形轮廓顶点经纬度坐标 # 从lon_min到lon_max沿lat_min画弧然后沿经线回到北极点 theta np.linspace(lon_min, lon_max, 100) verts [] # 先沿纬线圈lat_min从lon_min到lon_max for lon in theta: verts.append((lon, lat_min)) # 再沿经线lon_max从lat_min到90 for lat in np.linspace(lat_min, 90, 30): verts.append((lon_max, lat)) # 再沿经线lon_min从90回到lat_min for lon in theta[::-1]: verts.append((lon, 90)) # 这样会形成闭合轮廓吗注意我们其实需要从90度到lat_min沿lon_min补充一段其实这个顶点顺序有点乱。更稳妥的做法是用shapely生成一个多边形对象然后转成Path但CarTop边界又要求必须是matplotlib.path.Path对象。下面这个办法最直接、也最不容易出错def make_sector_path(lon_min, lon_max, lat_min, n100): 生成从lat_min到北极点、经度范围为[lon_min, lon_max]的扇形路径。 返回的Path在经纬度坐标系下。 # 沿着纬线圈从lon_min到lon_max lons1 np.linspace(lon_min, lon_max, n) lats1 np.full_like(lons1, lat_min) # 沿着lon_max经线从lat_min到90 lats2 np.linspace(lat_min, 90, n) lons2 np.full_like(lats2, lon_max) # 沿着纬线圈接近北极点从lon_max回到lon_min实际收缩成一点 lons3 np.linspace(lon_max, lon_min, n) lats3 np.full_like(lons3, 90.0) # 沿着lon_min经线从90回到lat_min lats4 np.linspace(90, lat_min, n) lons4 np.full_like(lats4, lon_min) lons np.concatenate([lons1, lons2, lons3, lons4]) lats np.concatenate([lats1, lats2, lats3, lats4]) # 用Path构造 return mpath.Path(np.column_stack([lons, lats]))注意上面这个Path是在经纬度坐标系近似PlateCarree下构建的所以创建Axes之后还需要把它转换到投影坐标系path make_sector_path(lon_min, lon_max, lat_min) # 在投影坐标下重新生成这个Path proj_path proj.project_geometry(sgeom.LineString(path.vertices))等一下project_geometry需要一个shapely geometry对象而LineString需要至少两个点并且要求所有点的经度不能产生跨180度的问题。如果你的扇形跨了日期变更线比如从160度到-160度直接构造LineString会得到一条横穿全球的线这时候需要把多边形拆成两段分别处理。这是后话我先给一个不跨日期线的一般用法import shapely.geometry as sgeom def sector_polygon(lon_min, lon_max, lat_min): 返回经纬度坐标下的扇形shapely Polygon n 100 lons1 np.linspace(lon_min, lon_max, n) lats1 np.full_like(lons1, lat_min) lons2 np.linspace(lon_max, lon_min, n) lats2 np.full_like(lons2, 90.0) coords list(zip(lons1, lats1)) list(zip(lons2, lats2)) return sgeom.Polygon(coords) polygon sector_polygon(lon_min, lon_max, lat_min) proj_polygon proj.project_geometry(polygon) # 投影到目标坐标系 boundary_path mpath.Path(np.asarray(proj_polygon.exterior.coords))这还不够因为proj.project_geometry返回的坐标已经是投影坐标米而ax.set_boundary需要在数据坐标即投影坐标下应用这样刚好对得上。最后应用边界ax.set_boundary(boundary_path, transformproj)等等这种用法有个版本差异问题旧版Cartopyset_boundary不强制transform参数新版如果传了transform而数据坐标不是投影数据会导致奇怪偏移。实际我测试下来最稳妥的写法是直接ax.set_boundary(boundary_path)set_boundary默认就采用当前Axes的投影坐标系来解读路径坐标所以不需要再传transform。这一点很多人会踩坑特意标注一下。加上地图要素ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linestyle:, linewidth0.6) ax.gridlines(draw_labelsTrue, linewidth0.5, alpha0.5)最终出来的效果北极点附近一个完整的扇形边界是笔直的经线侧边和一段弧线底边海岸线、国界和网格线都被裁剪在这个扇形里。这个就是标题要的“扇形区域图”。3.2 细节打磨为什么我的图右上角会出现“破洞”或者“飞线”用上面的方法跑通之后很多人会遇到一个老问题扇形以外的海岸线、网格线还是会穿过边界“漏”出来一部分像是什么东西从扇形边缘飞出去。这是因为set_boundary只限制了Axes的绘制范围但Gridline、Coastline这些Feature的绘制方式在某些版本里并不会严格被Boundary裁剪。最简单的解决方案关掉draw_labels或者手动控制网格标签的位置。我在实际项目中更喜欢的做法是关闭自动label叠加手动标注gl ax.gridlines(draw_labelsFalse, linewidth0.5, alpha0.5, linestyle--)然后手工加上想要的经纬度标注特别是扇形边缘的经线、底边的纬线。具体做法在后面的章节里再说。3.3 数据填充把格点数据放到扇形图里画边界只是第一步实际项目里真正核心的是往扇形里填你的物理场。比如我有一个T2M2米气温的全球格点数据lon和lat都是1度间隔的二维数组现在只想显示北纬60到90度、东经60到120度的区域。第一步把数据坐标转换到投影坐标之前先裁剪数据本身的范围以减小绘制压力避免NorthPolarStereo投影时把南北极附近已经重复的格点重复投影。第二步绘制时一定要带transformccrs.PlateCarree()我见过太多人因为漏了这个参数导致填色错了几万公里lon2d, lat2d np.meshgrid(lons, lats) # 数据裁剪 mask (lat2d lat_min) (lon2d lon_min) (lon2d lon_max) data_masked np.ma.masked_where(~mask, data) cf ax.contourf(lon2d, lat2d, data_masked, levels20, cmapRdBu_r, transformccrs.PlateCarree()) fig.colorbar(cf, orientationhorizontal, pad0.05, aspect40)这里注意虽然是北纬60度以上区域但如果数据分辨率很高比如0.1度一次性投影全部格点会非常慢。一个常用的技巧是只数据裁剪后用block方式降低分辨率绘制或者直接scatter少量站点。等值线填色其实对分辨率敏感度不高可以先做1度格点再contourf。对于海冰、SST等有明显边缘的数据contourf比pcolormesh更平滑并且可以直接加extendboth控制溢出色标cf ax.contourf(lon2d, lat2d, data_masked, levelsnp.linspace(0, 100, 21), cmapBlues, extendboth, transformccrs.PlateCarree())4. 地图要素与网格标签让扇形图“像”一张学术地图4.1 海岸线和边界线亮度和宽度的控制画扇形图时地图要素需要比常规地图更谨慎因为极区地形复杂海岸线支离破碎默认设置往往要么太粗影响数据读图要么太淡导致根本看不清。我常用的一套基线配置ax.add_feature(cfeature.LAND, color#e8e0d0, zorder1) ax.add_feature(cfeature.OCEAN, color#b0d4e8, zorder1) ax.add_feature(cfeature.COASTLINE, edgecolorblack, linewidth0.8, zorder2) ax.add_feature(cfeature.BORDERS, linestyle:, edgecolorgray, linewidth0.5, zorder2)几个重点提醒LAND/OCEAN的颜色别用纯白纯蓝因为后面一旦叠加半透明数据比如透明度0.7的填色底图颜色会透出来形成奇怪的混合色。使用浅的暖色和冷色能保证数据层的辨识度。如果只是做灰度打印直接用color0.7作为陆地色更稳妥。zorder务必设置尤其当数据填色zorder默认是1如果底图要素也是1可能出现海岸线被数据盖住的尴尬局面。4.2 经纬网与标签手动标注边缘经度自动gridlines(draw_labelsTrue)在极坐标投影下经常出现标签叠成一团或者出现在扇形外部的情况。我的经验是全部手动控制。先关闭自动标签只画网格gl ax.gridlines(crsccrs.PlateCarree(), draw_labelsFalse, linewidth0.4, colorgray, alpha0.6, linestyle--)然后标注扇形边界位置。比如扇形左边经度是60E右边经度是120E底边纬度是60N。我可以在投影坐标系下的对应位置通过ax.text手工写# 在底边中点标注纬度 ax.text((lon_minlon_max)/2, lat_min, 60°N, transformccrs.PlateCarree(), hacenter, vabottom)也可以把这些标注放在扇形左下和右下角ax.text(lon_min, lat_min2, 60°E, transformccrs.PlateCarree(), hacenter, vacenter) ax.text(lon_max, lat_min2, 120°E, transformccrs.PlateCarree(), hacenter, vacenter)这里有个小坑ax.text的transformccrs.PlateCarree()表示我指定的坐标是经纬度它会自动换算到投影坐标。你需要确保文本没被边界裁剪掉。默认情况下在Axes边界外的text是不会显示的所以文本位置要稍微向扇形内部偏移一点比如lat_min2。4.3 北极点标记与极点周边“死区”另一个常见需求是在图上标出北极点。北极点在NorthPolarStereo投影下就是圆心位置直接用ax.plot(0, 90, o, markersize6, colorred, transformccrs.PlateCarree()) ax.text(0, 89.5, NP, transformccrs.PlateCarree(), hacenter)但要注意如果你的数据网格里有准确的90N格点很多全球数据在90N处只有一个点投影后是一个极小的区域放大后会看到一个孤立的色块或者被渲染成一条线。这是正常现象。处理办法是在绘制之前把数据中的90N那一行极其靠近北极点的数据mask掉或者用contourf绘制使它成为扇形最顶端的一个平滑色块而不是一个点。5. 极端边界情况跨180度经线和半球范围5.1 经度窗口跨日期变更线怎么办很多时候你关心的区域并不老老实实待在0到360度的范围里。比如研究楚科奇海和白令海需要显示从东经160度到西经150度这个窗口恰好跨了180度经线。直接用前面的np.linspace(lon_min, lon_max)会得到一段横穿太平洋全球的长线导致多边形异常。解决办法有两种第一种将经度窗口统一写为0到360度。比如东经160度记为160西经150度记为210于是窗口是[160, 210]这样就不跨180度了。相应地进行经纬度范围判断时也需要把数据经度做统一变换lon360 np.where(lon 0, lon 360, lon) mask (lat lat_min) (lon360 160) (lon360 210)第二种拆分成两个窗口比如160到180、-180到-150分别画但这样处理起来繁琐且容易在接口处出现拼接缝隙。我推荐第一种操作简单性能损耗最低。5.2 宽扇形从90E到120W如果经度范围跨度太大比如150度以上扇形明显变宽整个图会变得非常扁。这时需要调整画布figsize的比例。比如经度跨度为100度时用figsize(8, 5)经度跨度为40度时用figsize(4, 6)。比例失控会导致经线看起来倾斜严重不美观。一个经验公式扇形“展开角”越大画布宽度比例应相应增加但别超过21否则视觉上整个北极被拉得不成形。5.3 当投影中心不在本初子午线默认NorthPolarStereo(central_longitude0)把0度经线放在正下方。如果你希望扇形区域看起来“正”。比如你关注的区域是180度经线附近白令海那最好设置central_longitude180这样180度经线会位于扇形正下方而不是卷曲到侧面proj ccrs.NorthPolarStereo(central_longitude180)实测效果变化很大尤其是跨大经度范围时。推荐做法是先整体看数据分布的质心经度把central_longitude设为那个值附近扇形图会显得更“正”。6. 进阶把扇形图嵌入更大的极区场景6.1 与整个北半球图拼接有些论文需要在一张图里同时展示“全北极范围”和“局部放大扇形”。我用两种方式实现过子图方式用fig.add_subplot(1, 2, 1)画全北极fig.add_subplot(1, 2, 2)画局部扇形两个Axes使用同一个投影ccrs.NorthPolarStereo但第二个Axes的边界范围不同。插图方式在扇形图内部添加一个小坐标轴用于指路类似“地图中的地图”。第二种方式在论文中更常见因为可以直观地告诉读者扇形区域在北极的位置。实现代码ax_sub fig.add_axes([0.72, 0.75, 0.2, 0.2], projectionccrs.NorthPolarStereo()) ax_sub.add_feature(cfeature.COASTLINE, linewidth0.4) ax_sub.set_extent([-180, 180, 65, 90], crsccrs.PlateCarree()) # 画一个表示扇形边界的矩形 ax_sub.add_patch(plt.Rectangle((lon_min, lat_min), lon_max-lon_min, 90-lat_min, facecolorred, alpha0.3, transformccrs.PlateCarree()))注意这个add_patch里的transformccrs.PlateCarree()是必需的否则矩形会被放在投影坐标里位置完全错乱。6.2 高分辨率背景与在线瓦片如果觉得Cartopy自带的GSHHS海岸线不够精细可以叠加Cartopy的img_tiles功能加载在线底图。但从实际体验看在一些网络不稳的环境中在线瓦片加载很慢而且极区瓦片覆盖并不完整。我更推荐本地自然地球Natural Earth矢量数据或分辨率更高的GADM国界数据。Cartopy自带的是NaturalEarthFeatureax.add_feature(cfeature.NaturalEarthFeature(physical, land, 10m, edgecolorblack, facecolor#e8e0d0))注意下载10m分辨率数据需要联网如果离线环境默认用50m即可。6.3 加入站点散点和轨迹线网格数据之后站点和航线散点是最常见的叠加元素。站点绘制时同样必须指定transformccrs.PlateCarree()。比如ax.scatter(site_lon, site_lat, s30, cred, edgecolorsblack, linewidths0.5, transformccrs.PlateCarree(), zorder5)航线直接把经纬度序列连线ax.plot(route_lons, route_lats, colorblue, linewidth1.5, transformccrs.PlateCarree(), zorder4)这里有一个极区特有的大坑如果你的轨迹线有跨180度的经度跳变比如从179到-179在投影坐标下会画出一条横穿整张图的直线。处理办法与5.1相同先统一成0~360或者拆分线段。7. 常见报错与排查链路7.1ValueError: Cannot handle non-finite values或Invalid transform sphere这类报错通常有两种原因一是数据中出现了NaN二是经纬度网格中有不符合PlateCarree的坐标如纬度超过90或者经度超过360。排查路径先检查数据np.isfinite(data).all()如果False进行差值或剔除。检查经纬度范围纬度必须在-90到90之间经度在-180到360之间。检查是否有坐标数组与数据shape不一致遮罩矩阵写错导致坐标轴对齐出错。7.2 扇形边界裁剪后内侧出现白色缺口经常看到扇形底边有一段白色空隙海岸线但是海陆底色没有填进去。这多半是LAND要素的zorder低于扇形边界背景层。先尝试把ax.add_feature(cfeature.LAND, ...)放在set_boundary之后执行如果不行就在set_boundary之后用ax.add_patch补一层扇形内部背景色。7.3 输出PDF时扇形边缘出现细白线在PDF输出中set_boundary裁剪的路径边缘经常会出现一条极细的空白线这是渲染器的抗锯齿问题。给Axes加一个宽度足够的patch边框可以压制ax.patch.set_linewidth(1) ax.patch.set_edgecolor(black) # 或者和背景色一致7.4 Cartopy和Shapely版本冲突project_geometry在不同版本上的行为略有差异。旧版Cartopy0.18左右在投影某些跨越180度的多边形时可能会生成异常形状。遇到这种情况先升级Cartopy到最新版本0.21再把数据经度统一到0~360。我已经多次靠这个组合解决“多边形突然扭曲成一个圆”的问题。8. 我的选型判断与一些私人经验最后聊几句选型层面的体会。用Cartopy画北极扇形区域图核心工作量其实不在“画”本身而在“把投影和裁剪彻底想明白”。很多人绕不开Basemap因为Basemap的API更直观设置llcrnrlat, urcrnrlat就能直接限定范围但它早在2020年就停止维护了在新版Matplotlib上频繁报np.float相关的兼容性错误。Cartopy虽然初期接触时总觉得“transform这东西真麻烦”可一旦形成心智模型后面画任何投影图都是一通百通。我个人在实际项目中最终沉淀下的固定配置组合是这样的投影优先ccrs.NorthPolarStereo(central_longitude需要对齐的经度)扇形路径用shapely.geometry.Polygon构造数据坐标用0~360的经度体系避免跨180度问题数据绘制时所有经纬度散点、填色、等值线全部带transformccrs.PlateCarree()网格标签用手动ax.text标注不用gridlines(draw_labelsTrue)省去极区标签乱飞的烦恼;输出时优先fig.savefig(xxx.pdf, bbox_inchestight)矢量格式能避免放大后锯齿。这两个“手动”替代“自动”的习惯帮我省下过无数用来调样式的时间。反正极区图的标签就那么几根手写几分钟自动调整反而常常半小时还没看顺眼。如果你只是临时画一张给自己看的数据快照那上面章节里的代码足够了但如果你要做得更细比如论文出版图、甚至要出动画纬度带选择、经度中心对齐、标注位置、色彩方案这些细节是值得一调再调的。极区图最怕一带而过它比中低纬度图更容易因为投影变形引发误读。做图的人和看图的人得在同一个坐标系里对话你的图才真正站得住。