How barometer.today draws its isobars

How barometer.today draws its isobars

barometer.today 如何绘制等压线

barometer.today shows air pressure for any city: the reading now, the change since yesterday, and the week ahead. It also has a world isobar map, animated hour by hour from a day back to four days ahead. This post explains what isobars are and the algorithms that we use to draw them. The code is Python with numpy, scipy, contourpy and shapely, and anyone drawing contour maps from gridded data will meet the same problems.

barometer.today 可以显示任何城市的气压:当前的读数、自昨天以来的变化以及未来一周的趋势。它还提供了一张全球等压线地图,以小时为单位进行动画演示,时间跨度从过去一天到未来四天。本文将解释什么是等压线,以及我们绘制它们所使用的算法。代码基于 Python,使用了 numpy、scipy、contourpy 和 shapely 库。任何从网格数据绘制等值线图的人都会遇到同样的问题。

What an isobar is

什么是等压线

Air pressure varies from place to place, and a weather map shows this with isobars - lines joining points of equal pressure. Where lines crowd together pressure changes fast over a short distance and the wind is usually strong. Where they spread out the air is calmer. The rings close around centers of high pressure marked H and low pressure, marked L. A deep low with tight rings is a storm. A broad high with loose rings is settled weather.

气压因地而异,天气图通过等压线(连接气压相等点的线)来展示这一点。当等压线密集时,气压在短距离内变化剧烈,通常风力强劲;当等压线稀疏时,空气则较为平静。环线围绕着高压中心(标记为 H)和低压中心(标记为 L)闭合。环线紧密的深低压通常意味着风暴,而环线稀疏的广阔高压则代表天气稳定。

The pressure on such a map is sea-level pressure. A barometer in Denver reads about 835 hPa because it sits 1,600 m up, which says nothing about the weather. So the reading is converted to what it would be at sea level, and the map becomes comparable everywhere. With one exception, covered below.

地图上的气压是指海平面气压。丹佛的气压计读数约为 835 hPa,因为它位于海拔 1,600 米处,这并不能反映天气情况。因此,读数会被转换为海平面气压,从而使地图在各地具有可比性。但有一个例外,下文将详细说明。

The data

数据来源

The map is drawn from the ECMWF global forecast. It arrives as a grid of sea-level pressure values every 0.25 degrees, about 28 km, for the whole planet. The model’s own time step is three hours; Open-Meteo interpolates that to hourly values. We fetch it from Open-Meteo, which we run on our own server for all SimpleMeteo.com sites. uvi.today and pollen.today use the same instance. The model updates four times a day. After each update a job produces one frame per hour, about 120 in total.

该地图基于欧洲中期天气预报中心(ECMWF)的全球预报数据绘制。数据以全球网格形式提供,分辨率为 0.25 度(约 28 公里)。模型本身的时间步长为三小时,Open-Meteo 将其插值为每小时数值。我们从 Open-Meteo 获取数据,并将其运行在我们自己的服务器上,供所有 SimpleMeteo.com 站点使用。uvi.today 和 pollen.today 也使用相同的实例。该模型每天更新四次。每次更新后,后台任务会生成每小时一帧的图像,总计约 120 帧。

Smoothing

平滑处理

A raw model grid is jagged cell to cell, so lines traced straight from it wobble. We use a Gaussian blur with a sigma of one cell, 28 km. Wider blurs give calmer lines but flatten the extremes: with a 110 km blur, a hurricane with a 936 hPa core comes out as a mild 993 hPa low. It could produce a map that looks fine and is wrong.

原始模型网格在单元格之间存在锯齿,因此直接绘制的线条会产生抖动。我们使用 sigma 为一个单元格(28 公里)的高斯模糊进行处理。更宽的模糊虽然能使线条更平滑,但会削弱极值:使用 110 公里的模糊,一个核心气压为 936 hPa 的飓风可能会被平滑成一个 993 hPa 的温和低压。这会导致地图看起来很美观,但实际上是错误的。

The blur has to know about mountains. Over a high plateau sea-level pressure is a calculation about air that does not exist, and it usually comes out several hPa too high. A plain Gaussian would spread those inflated values into the real field around the plateau. Along the Himalayan front this moves lowland values by more than an isobar interval.

模糊处理必须考虑山脉的影响。在高海拔高原上,海平面气压是对不存在的空气进行的计算,通常会比实际高出几个 hPa。普通的高斯模糊会将这些虚高的数值扩散到高原周围的真实区域。在喜马拉雅山脉前沿,这会导致低地数值的偏移超过一个等压线间隔。

The fix is a normalized convolution. Blur the field with the plateau cells set to zero, blur a mask of ones and zeros the same way, and divide one by the other. Each output cell is then the weighted mean of only the valid cells nearby.

解决方法是归一化卷积。将高原单元格设为零对场进行模糊,同时以相同方式对由 1 和 0 组成的掩码进行模糊,然后将两者相除。这样,每个输出单元格就仅是附近有效单元格的加权平均值。

Where the mask is

掩码的位置

The mask is ground above 2,500 m. It comes from a digital elevation model, smoothed to 55 km so that a single peak does not count, then only massifs of at least 265,000 km² are kept. Four regions qualify: the Andes, the Tibetan plateau with the Himalaya and the Pamir, the Greenland ice sheet and Antarctica.

掩码覆盖了海拔 2,500 米以上的地面。它源自数字高程模型,经过 55 公里的平滑处理以排除孤立山峰,并仅保留面积至少为 265,000 平方公里的山地。符合条件的区域有四个:安第斯山脉、青藏高原(含喜马拉雅山脉和帕米尔高原)、格陵兰冰盖和南极洲。

Over the mask the isobars are still drawn, but dashed. The map stays continuous and the reader is not left with a hole in Asia, while the dashes say the values here are extrapolated. Each line segment is tested against the mask individually, so an isobar crossing the Andes is dashed for exactly the part over high ground.

在掩码区域内,等压线依然会被绘制,但以虚线表示。这样地图保持了连续性,读者不会在亚洲看到一个“空洞”,同时虚线也提示此处数值为外推值。每个线段都会单独与掩码进行比对,因此穿过安第斯山脉的等压线只有在高海拔区域的部分才会显示为虚线。

Tracing the lines

绘制线条

Contours are traced with marching squares, the standard algorithm that matplotlib uses. Each grid cell has a value at its four corners. For a given level, say 1012 hPa, the algorithm marks each corner as above or below, which gives 16 possible patterns, and each pattern says where the line enters and leaves the cell. The crossing point along an edge is found by linear interpolation between the corner values. Done for every cell, the pieces join into continuous lines. Done for every level, you have the map.

等值线使用“移动方块”(marching squares)算法绘制,这是 matplotlib 使用的标准算法。每个网格单元的四个角都有一个数值。对于给定的气压水平(例如 1012 hPa),算法会标记每个角是高于还是低于该值,从而产生 16 种可能的模式,每种模式决定了线条进入和离开单元格的位置。边缘上的交叉点通过角点数值的线性插值求得。对每个单元格执行此操作,碎片就会连接成连续的线条;对每个气压水平执行此操作,就得到了完整的地图。

A world map adds a twist: longitude wraps around. The algorithm sees a rectangle with left and right edges and would stop every line at 180 degrees. So the grid is padded with two degrees of columns copied from the opposite edge, contoured, and the lines clipped back to the world. They then cross the Pacific date line without a seam.

世界地图增加了一个复杂之处:经度是循环的。算法将地图视为一个具有左右边缘的矩形,会在 180 度经线处中断所有线条。因此,我们在网格两侧填充了从对面边缘复制的两度列数据,进行等值线计算,然后再将线条裁剪回世界范围。这样,线条就能无缝穿过太平洋日期变更线。

Each line is then simplified with the Douglas-Peucker algorithm at 0.05 degrees. It removes points that lie within that distance of a straight line between their neighbours, so a nearly straight 900 km isobar collapses to three points while a tight ring around a storm keeps its shape. A frame ends up as a few hundred kilobytes of GeoJSON, under 100 KB compressed.

随后,每条线都使用 Douglas-Peucker 算法以 0.05 度的精度进行简化。该算法会移除那些位于相邻点连线一定距离内的点,因此一条近乎笔直的 900 公里等压线会被简化为三个点,而围绕风暴的紧密环线则能保持其形状。最终,每一帧数据约为几百 KB 的 GeoJSON,压缩后不到 100 KB。

Finding the H and L

寻找高压(H)和低压(L)中心

Centers are found on the smoothed field with min and max filters: a cell is a candidate if it is the extreme within a window of about 850 km. The window is built per band of latitude, so it is a real distance and not a cell count. A fixed cell count would cover only half the distance at 60 degrees north that it covers at the equator.

中心点是在平滑后的场上通过最小值和最大值滤波器找到的:如果一个单元格在约 850 公里的窗口内是极值,它就是一个候选点。窗口是按纬度带构建的,因此它代表的是实际距离而非单元格数量。如果使用固定的单元格数量,在北纬 60 度覆盖的距离将仅为赤道处的一半。

A candidate must also stand out. The field within the window must range by at least 3 hPa, or 5 hPa inside 20 degrees of the equator, where the tropical atmosphere has a natural twice-daily tide of 2 to 4 hPa that would otherwise put a letter at every ripple.

候选点还必须足够显著。窗口内的气压范围必须至少达到 3 hPa,在赤道 20 度以内则需达到 5 hPa,因为热带大气存在自然的每日两次、幅度为 2 到 4 hPa 的气压潮汐,否则每一个波动都会被标记上字母。

Two details matter in practice. The masked cells are filled with minus infinity for the max filter and plus infinity for the min filter, never NaN. Scipy’s filters treat a NaN differently depending on where it falls in the window, and the result is that a real low on the Antarctic coast can silently vanish. And no letter is placed on ground above…

在实践中,有两个细节很重要。掩码单元格在最大值滤波器中被填充为负无穷大,在最小值滤波器中被填充为正无穷大,绝不使用 NaN。Scipy 的滤波器根据 NaN 在窗口中的位置对其处理方式不同,这可能导致南极海岸的一个真实低压中心悄无声息地消失。此外,字母不会放置在海拔高于……的地面上。