2026/10/17 10:15:56

Basemap气象绘图实战:底图叠加与白化maskoceans完整代码解析

Basemap气象绘图实战:底图叠加与白化maskoceans完整代码解析 做气象绘图的谁没被basemap折磨过这个藏在Matplotlib工具箱里的老家伙虽然是十年前的“旧物”但至今还有大量科研脚本、毕业论文代码在用尤其在做大气科学可视化时底图叠加和“白化”这两件事几乎绕不开它。这个系列第二篇我就把实战中最常用的一套流程掰开讲清楚怎么用basemap叠加地图底图怎么给气象要素场做白化效果并给出一份可以直接跑的完整代码。文章适合正在做天气分析、气候统计、模式后处理的学生和科研人员也适合那些接手了老代码、被迫继续维护basemap脚本的朋友。1. 为什么是basemap为什么需要白化1.1 basemap在气象绘图里的真实地位basemap本质上是一套把地理投影、海岸线、行政边界和经纬网整合进Matplotlib的工具包。早期做气象图大家几乎离不开它先用Basemap初始化一个带投影的地图对象再调用drawcoastlines、drawcountries、drawparallels这些方法把地理要素铺到底层最后把模式数据或者站点观测画上去。它的优点是API非常直接写起来像“搭积木”对于非地理信息专业出身的气象研究者来说门槛比手写投影变换低得多。这里就产生了一个现实问题basemap官方早已停止更新很多新装环境里已经很难直接pip安装。但存量代码实在太多高校课题组里随便翻一个老项目可能就是from mpl_toolkits.basemap import Basemap开头。这些脚本2015年能跑2025年只要环境对也能跑。我自己的体会是与其劝人彻底忘掉basemap不如把它讲透毕竟底图投影、坐标变换、掩膜白化这些思路迁移到cartopy之后是完全相通的学会basemap等于先把绘图“骨架”建了一遍。1.2 “白化”到底是在白化什么气象要素场经常覆盖整个经纬度网格但海洋上的观测资料往往是缺测的或者模式输出的海温、降水在海上没有实际业务意义。如果不做任何处理等值线填色会一路画到海面上图上看起来就像大陆的降水“长”到了海洋里既不符合科学事实也影响读图判断。所谓白化就是把不需要显示的区域用掩膜遮掉让填色只落在目标区域。在basemap里这个操作的对象不是绘图结果而是数据本身先根据经纬度判断某个格点是在陆地还是海洋再把海洋格点标记为无效值等值线遇到无效值就不会绘制视觉上就形成了“只保留陆地”的效果。有人把这种操作叫mask有人叫maskout也有人沿用了气象绘图里的习惯叫“白化”其实都是同一件事。理解了这一点后面所有基于掩膜的操作就都清楚了。1.3 方案选型为什么还要用basemap既然basemap已经停止维护如果完全从零开始一个新项目我通常建议直接上cartopy。但现实世界里大量任务是复现历史图、维护老代码、沿用课题组既有模板这种情况下硬搬到cartopy反而会引入新的问题。另外basemap自带maskoceans、maskland这类非常顺手的白化方法老版本代码里只要一行调用迁移时还得自己实现一遍同等功能投入不小。所以我的建议是如果只是临时画几张图且环境干净优先cartopy如果手头已有成熟basemap脚本或者要在旧环境里批量出图就安心把basemap用明白。技术选型从来不是“必须用最新”而是“在约束条件下把活干好”。这篇文章后面所有代码都基于basemap但我也会在末尾交代一些迁移到cartopy时的对应思路。2. 环境准备与数据组织2.1 安装basemap的正确姿势与版本约束安装basemap是很多人遇到的第一个拦路虎。直接pip install basemap在较新的Python版本里大概率失败原因在于basemap依赖GEOS地理库并且旧版包与新版numpy、Python的ABI不兼容。实际中我试过最稳的做法是用conda创建独立环境指定Python 3.8以内然后通过conda-forge安装conda create -n py38_basemap python3.8 conda activate py38_basemap conda install -c conda-forge basemap如果项目还需要读取NetCDF数据建议一起装上xarray和netcdf4conda install -c conda-forge xarray netcdf4 matplotlib这里有一个很多人忽略的细节basemap提供了mpl_toolkits.basemap这个模块但底层依赖GEOS的Python封装。通过conda安装时conda会帮我们把geos的库文件一起装好而pip安装往往缺这一步导致import时报libgeos_c.so找不到。所以不要在pip这边硬磕换conda会省很多时间。2.2 数据准备经纬度网格、维序和缺测处理做白化之前先把数据组织好。气象网格数据最常见的格式是NetCDF里面通常有longitude、latitude和变量三个维度变量的维度顺序可能是(time, level, lat, lon)。为了绘图方便我会先用xarray读入再调整成二维网格需要的形状import xarray as xr ds xr.open_dataset(precip_example.nc) var ds[precip].isel(time0) # 取某个时次 lons var[lon].values lats var[lat].values data var.values # 此时是 2D: (lat, lon)这里要特别注意维度顺序。basemap的maskoceans需要二维的经纬度网格如果数据是(lat, lon)二维数组构造网格时要保证和数组一致lon2d, lat2d np.meshgrid(lons, lats)另外如果数据里有NaN或者其他缺测标记最好在画图前先处理干净。maskoceans主要处理海陆掩膜不会替你把所有缺测修好。常见做法是先判断合法范围再把无效值设为NaN这样等值线图也会自动跳过这些区域。还有一个容易踩的点是经度范围。有的数据用0到360表示东经有的用-180到180。如果数据和底图投影区域不一致白化效果会完全错位。建议在进入绘图流程前统一经度范围必要时用np.where(lon 180, lon - 360, lon)把经度转到-180到180。2.3 投影选择与中心参数怎么定basemap支持很多投影气象上最常用的三个cyl等距圆柱投影适合全球或低纬大范围merc墨卡托投影适合热带和赤道附近区域lcc兰勃特等角圆锥投影适合中纬度天气图尤其是冷暖气团、锋面、降水这类过程分析。区域图我基本都用lcc。初始化的关键参数包括中心经度lon_0、中心纬度lat_0、地图宽度width和高度height以及标准纬度lat_ts。m Basemap(projectionlcc, lon_0105, lat_035, width2400000, height2400000, resolutionh, lat_ts35)很多人不明白width和height的单位。这两个参数是米不是经纬度。粗略换算可以记一个数纬度方向1度大约111公里经度方向在纬度为φ时大约111*cos(φ)公里。如果我想画一个东西向跨20度、南北向跨15度的区域中心纬度35度估算宽度大约是20*111*cos(35°) ≈ 1800公里高度大约是15*111 ≈ 1665公里。实际画图时宁可留点余量把width和height稍微放大10%到20%避免数据被地图边界裁掉。resolution控制海岸线精度可选c粗糙、l低、i中等、h高、f精细。区域图建议至少用h否则海岸线锯齿很明显放大了更难看。lat_ts一般与lat_0一致也可以根据自己的标准纬线调整。3. 先画底图再画数据核心实现流程3.1 底图初始化与地图要素叠加底图初始化的本质是把经纬度坐标转换到一个投影平面坐标系。Basemap对象创建后可以用它返回的转换函数把经纬度数组变成画图坐标x, y m(lon2d, lat2d)这一步做好后后面所有绘图操作都在x、y上进行而不是直接用经纬度。很多人画完图发现数据位置偏了、拉伸了十有八九是忘了做这一步。地图要素的叠加有自己的固定套路通常包括海岸线、国界、经纬网和边界m.drawcoastlines(linewidth0.8, colorblack) m.drawcountries(linewidth0.8, colorblack) m.drawparallels(np.arange(20, 50, 5), labels[1, 0, 0, 0], fontsize10) m.drawmeridians(np.arange(80, 130, 5), labels[0, 0, 0, 1], fontsize10) m.drawmapboundary(fill_colorlightsteelblue)drawparallels和drawmeridians里的labels参数控制哪一侧显示刻度标签四个值分别对应左、右、上、下。区域图我只开左侧纬度和下侧经度这样图上不会到处飘数字。有一个顺序问题必须提醒底图要素之间相对顺序会影响观感。海岸线和国界建议在填色图之后、标注之前叠加这样数据不会盖住边界线边界线又不会压住最上层信息。3.2 数据投影与海洋白化核心代码这是整篇文章最核心的一段。以降水场为例如果用contourf直接画原始数据填色会铺满整个矩形区域海洋上也全是颜色。要在basemap里做白化标准做法是用maskoceans方法data_masked m.maskoceans(lon2d, lat2d, data)这个方法内部会根据经纬度判断每个格点是否在海洋上是海洋就置为掩膜值。mask过后的数组传给contourf等值线只在陆地格点上绘制。完整示例代码如下import numpy as np import xarray as xr import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap # 1. 读取网格数据 ds xr.open_dataset(precip_example.nc) var ds[precip].isel(time0) lons var[lon].values lats var[lat].values data var.values # 2. 构造二维经纬度网格 lon2d, lat2d np.meshgrid(lons, lats) # 3. 处理无效值 data np.where(np.isfinite(data), data, np.nan) # 4. 初始化底图区域中心经度和纬度请替换成你的实际参数 fig plt.figure(figsize(8, 8)) m Basemap(projectionlcc, lon_0105, lat_035, width2400000, height2400000, resolutionh, lat_ts35) # 5. 经纬度转投影坐标 x, y m(lon2d, lat2d) # 6. 海洋白化 data_masked m.maskoceans(lon2d, lat2d, data) # 7. 绘制等值线填色 levels np.arange(0, 60, 5) cs m.contourf(x, y, data_masked, levelslevels, cmapBlues, extendboth) # 8. 叠加地图要素 m.drawcoastlines(linewidth0.8) m.drawcountries(linewidth0.8) m.drawparallels(np.arange(20, 50, 5), labels[1, 0, 0, 0], fontsize10) m.drawmeridians(np.arange(80, 130, 5), labels[0, 0, 0, 1], fontsize10) m.drawmapboundary(fill_colorlightsteelblue) # 9. 色标 cbar m.colorbar(cs, locationright, size5%, pad5%) cbar.set_label(Precipitation (mm)) # 10. 保存 plt.savefig(precip_with_mask.png, dpi150, bbox_inchestight)很多人会问maskoceans为什么输入的是经纬度而不是投影坐标因为海洋判断依赖真实地理信息只有经纬度坐标才能定位到海陆位置。实测下来如果传二维网格mask效果会比较稳定只传一维坐标可能会导致shape匹配问题建议都先np.meshgrid再调用。3.3 等值线、色标与出图参数怎么调contourf的几个参数直接影响成图质量。levels决定了等值线的分级气象图最好用整十整五的间隔方便业务判读。cmap的选择也要考虑场景降水用Blues或者BrBG气温用RdYlBu_r风场高度场用jet也常见但科研投稿越来越忌讳彩虹色我一般优先用Matplotlib里那些色觉友好型色板。extendboth很关键。当数据最大值超过levels最大值时色标顶部会多一个箭头表示“还有更大值”数据最小值低于levels最小值时底部也有箭头。如果不加这个参数超出范围的数据会被静默丢弃容易让人误读成缺测。色标有两种写法。基础写法是plt.colorbar(cs, shrink0.8, pad0.05)比较通用basemap环境里也可以用m.colorbar(cs, locationright, size5%, pad5%)好处是自动把色标贴在地图边上。我个人更常用后者因为在区域图里色标占用的相对位置更可控。保存图片时dpi150基本能满足论文初稿投稿可以到300。bbox_inchestight会把多余白边裁掉不过要注意如果还想再叠其他小图建议先不加这个参数否则坐标系边界会变化。4. 白化效果的进阶玩法与细节打磨4.1 反向白化与区域级联裁剪maskoceans是“把海洋mask掉”反过来如果要画海温、海洋变量需要“把陆地mask掉”basemap同样提供了对应方法maskland。用法完全一致sst_masked m.maskland(lon2d, lat2d, sst_data)这种反向白化在海温图、海浪图里非常常见。需要注意maskland的陆判精度同样受海岸线分辨率影响如果区域沿岸地形复杂建议把resolution提高到h以上再用。有时候需要更精细的区域裁剪比如只保留某个流域、某条海岸线以内。basemap没有直接提供的“多边形裁剪”函数但可以通过Matplotlib的路径对象实现先把目标区域边界做成一个多边形再判断每个网格点是否落在多边形内不在的点置为NaN。核心代码如下from matplotlib.path import Path # poly_lon, poly_lat 是目标区域边界经纬度 polygon Path(list(zip(poly_lon, poly_lat))) grid_points np.column_stack([lon2d.ravel(), lat2d.ravel()]) in_region polygon.contains_points(grid_points).reshape(lon2d.shape) data_region np.where(in_region, data, np.nan)这种做法的好处是不依赖海陆判断完全由你自己定义“有效区域”。做流域面雨量、城市群污染带分析时非常实用。边界点越多裁剪结果越精细但计算量也会上升网格特别大的时候建议先把边界抽稀到几百个点以内。4.2 中文字体、标注和数据叠加basemap本身不管中文字体中文乱码是Matplotlib的经典问题。用之前先设置字体import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [SimHei] # Windows 常见 plt.rcParams[axes.unicode_minus] False # 修正负号显示Linux环境下没有SimHei可以改成Noto Sans CJK SC或者WenQuanYi Micro Hei要确保系统已安装对应字体。字体设置要放在绘图脚本最开始否则个别元素已经渲染成方框了再改也来不及。地图上叠加站点观测是很常见的需求。站点经纬度是散点不能直接扔进plot要先通过底图对象转换slons [108, 110.5, 113] slats [30, 31.2, 33.5] sx, sy m(slons, slats) m.plot(sx, sy, o, colorred, markersize6)如果还要标站名用ax.text(sx[i], sy[i], 站名)时注意偏移一点避免文字盖住圆点。标题和色标标签里也可以放中文但保存时务必确认输出文件里的中文正常不行就换plt.savefig后再用图片工具检查一下。4.3 图层顺序为什么会影响成图观感同样的数据图层顺序不同观感差异很大。我见过不少新手把fillcontinents放在填色图之后画结果陆地色块把数据整个盖住只剩海洋上一圈颜色完全看不出降水分布。常规的图层顺序是先画海洋底色或陆地底色再画数据填色再画海岸线和国界最后画站点和标注。按照这个顺序陆地在数据层下面海岸线在数据层上面读图时既能看到数据分布又能看清地形边界。具体到代码m.fillcontinents(colorlightgray, lake_colorlightsteelblue) cs m.contourf(x, y, data_masked, levelslevels, cmapBlues) m.drawcoastlines(linewidth0.8)如果不想让陆地底色抢数据视觉可以把陆地颜色调得很浅比如#f5f5f5。海陆对比不要太强否则填色的色彩会显得灰蒙蒙。5. 踩坑实录常见问题与排查技巧5.1 安装和导入阶段的报错速查basemap的坑集中在两个阶段装不上和画不对。先整理一份高频报错速查表。报错或现象常见原因解决办法ModuleNotFoundError: No module named mpl_toolkits.basemapbasemap未安装或装到了别的环境确认当前conda环境按上文conda-forge方式安装ImportError: libgeos_c.so.1: cannot open shared object fileGEOS动态库缺失常发生在pip安装后改用conda安装basemap或手动补装geos安装时报Failed to build basemapPython版本过新旧版basemap源码不兼容创建Python 3.8环境再用conda-forgeAttributeError: Basemap object has no attribute colorbarbasemap版本过老或用户混淆了方法名改用plt.colorbar(cs, shrink0.8, pad0.05)如果你是在某个共享服务器上维护老项目建议不要轻易动系统级Python环境为项目单独建一个conda环境把basemap、xarray、matplotlib版本都固定好。这样后续别人复现脚本一克隆环境就能跑不会互相踩。5.2 白化与绘图阶段的典型坑装好环境只是开始画图阶段才是考验细节的地方。我把自己踩过以及帮别人排查过的坑整理成几个典型情况。第一个是白化边界呈锯齿状。这个太常见了尤其用resolutionl或i时岸线本身就不平滑掩膜也跟着粗糙。解决方法是提高分辨率到h同时把网格数据插值到更细的格点。如果还觉得边界不够好看可以在填色图之上再加一层drawcoastlines用线图遮住一部分锯齿视觉上会干净很多。第二个是maskoceans之后整块图全空白。出现这个先别怀疑方法本身多半是经纬度维序不匹配。比如数据是(lon, lat)顺序但你先meshgrid(lons, lats)得到(lat, lon)形状两者错位之后每个格点对应的经纬度全乱了mask自然把所有点都吞掉。排查技巧很简单打印data.shape和lon2d.shape、lat2d.shape确保三者完全一致。第三个是等值线填色范围“炸了”颜色分布和预想不一致。常见原因是levels设置没覆盖数据实际范围。比如数据最小值-5最大值100你却从0开始设levels小于0的部分会全部落在第一级色标底部还带箭头。建议在画图前先看np.nanmin(data)和np.nanmax(data)再决定levels。第四个是数据位置整体偏移或者被拉伸。这个问题几乎都是忘记做投影转换直接把经纬度数组传给了contourf或者传了二维经纬度但没用m()转换。底图坐标系是投影平面绘制任何数据都要先经过m(lon, lat)变换这个步骤不能省。第五个是经度跨0度或跨180度的数据画出来特别怪。比如数据经度范围是0到360底图经度范围是-180到180海陆判断就会错位看起来降水不在正确的位置上。处理办法在数据准备小节已经提过统一经度表示方式把大于180的值减360。5.3 如何用最快的方式定位问题绘图脚本出错时与其逐行看报错不如分段排查。我的经验是三层递进先跑一个最简底图不画数据只画海岸线和经纬网确认投影参数没问题再加数据但不加mask确认数据位置、色标、levels都正常最后再加maskoceans确认白化效果。这样每一步的输出都容易判断出了问题也知道是底图、数据还是mask环节。在这个排查过程中print是你的好朋友。打印数据的shape、经纬度范围、mask前后有效点数量变化往往一眼就能看出问题出在哪。比如mask前有效点1万个mask后只剩2000个那说明陆地格点占比本身就低如果mask后剩100个那就要怀疑经纬度顺序是不是错了。6. 写在最后老脚本与新脚本之间怎么选我个人在实际操作里的体会是basemap的“老”不代表“没有价值”它的掩膜思路、投影选择、图层组织方式放到今天依然成立。我手上有大量老脚本后来为了维护方便我把常用的basemap绘图逻辑封装成了一套函数库maskoceans那步做成了装饰器一行调用就能叠加白化。再往后接触cartopy时我发现只要理解了“经纬度转投影坐标→画数据→按海陆掩膜→叠加地理要素”这条链路迁移起来非常顺畅并没有想象中那么痛苦。最后再分享一个小技巧如果以后换环境、换新库别急着推翻老代码先把basemap里相关函数的参数对照表列出来比如drawcoastlines对应cartopy的feature.COASTLINEmaskoceans则可以用形态学填充的思路自己实现。这样你维护老项目的同时也能逐步把核心绘图逻辑抽离成与库无关的模块。说到底画图工具会变但你对数据的理解、对科学表达的要求才是真正值钱的那部分。