簡介一份面向氣象數(shù)據(jù)分析初學(xué)者與Python學(xué)習(xí)者的輕量實(shí)操示例包聚焦如何用Python完成氣象數(shù)據(jù)的讀取、清洗、分析與可視化展示覆蓋從CSV站點(diǎn)數(shù)據(jù)導(dǎo)入到溫度等要素變化圖繪制的完整鏈路。包內(nèi)共15個(gè)文件包括10個(gè)氣象觀測CSV樣本、4個(gè)Python腳本及1個(gè)說明文檔壓縮包整體僅14KB結(jié)構(gòu)小巧、便于快速上手。已有9028人學(xué)習(xí)下載適合想通過短小代碼了解pandas、matplotlib、seaborn等庫在氣象場景中實(shí)際用法的讀者。腳本按不同城市分組組織數(shù)據(jù)示例中既包含用pandas讀取站點(diǎn)記錄并繪制時(shí)間變化曲線也涉及scipy時(shí)間序列分析與周期性觀察可通過對照多個(gè)城市樣本理解地理差異對氣象要素的影響。項(xiàng)目目錄清晰代碼與數(shù)據(jù)一一對應(yīng)運(yùn)行后即得可視化結(jié)果邊看邊改即可遷移到自己的數(shù)據(jù)集是一份入門氣象數(shù)據(jù)分析與可視化流程的實(shí)用參考。1. 氣象數(shù)據(jù)分析與可視化先從NetCDF到一張能用的天氣圖接到一份氣象數(shù)據(jù)分析與可視化的活兒甲方給的是全國逐日氣溫NetCDF文件要求出三張圖今年與去年的距平對比、夏季熱浪過程的溫度演變、冬季冷空氣南下的空間分布。很多人拿到數(shù)據(jù)第一步就打開matplotlib畫圖結(jié)果是坐標(biāo)軸亂掉、缺測值把色條拉爆、時(shí)間序列里混進(jìn)多年數(shù)據(jù)。這篇文章講的就是如何用Python把這套氣象數(shù)據(jù)分析與可視化流程走通從讀取NetCDF原始數(shù)據(jù)到清洗缺測、校準(zhǔn)時(shí)間坐標(biāo)再到計(jì)算距平和區(qū)域平均最后用Cartopy畫帶底圖的空間場用matplotlib動(dòng)畫呈現(xiàn)時(shí)間演變讓圖表能放進(jìn)報(bào)告或大屏。適合剛把python入門過的分析崗新人也適合氣象、農(nóng)業(yè)、環(huán)境行業(yè)想自建分析流程的從業(yè)者。2. 環(huán)境與數(shù)據(jù)準(zhǔn)備裝好Python庫之后先讀懂?dāng)?shù)據(jù)再動(dòng)手做氣象分析最先要過的關(guān)不是算法而是數(shù)據(jù)讀取。氣象行業(yè)的原始數(shù)據(jù)大多是NetCDF格式一個(gè)文件里同時(shí)存在時(shí)間、緯度、經(jīng)度、氣壓層等多個(gè)維度和平時(shí)處理的表格結(jié)構(gòu)完全不是一回事。這一章先把環(huán)境怎么裝講清楚再規(guī)范讀取流程最后把數(shù)據(jù)清洗的幾個(gè)邊界條件交代完。2.1 依賴環(huán)境為什么是xarray加netCDF4加cartopy我常用的環(huán)境組合是Anaconda基礎(chǔ)環(huán)境加三個(gè)核心庫xarray負(fù)責(zé)多維數(shù)組和NetCDF文件的讀寫底層驅(qū)動(dòng)是netCDF4庫繪圖部分用matplotlib加cartopy。之所以不推薦直接用netCDF4庫讀數(shù)據(jù)是因?yàn)樗祷氐氖堑讓訑?shù)組對象手動(dòng)管理維度坐標(biāo)非常麻煩xarray把這層包了一層讀進(jìn)來就是帶維度名的DataArray后續(xù)按時(shí)間、按區(qū)域切數(shù)據(jù)都在明面上。pip install xarray netCDF4 matplotlib cartopy在Windows上建議順手把cftime也裝上后面處理帶hours since 1900-01-01這種參考時(shí)間格式的坐標(biāo)會(huì)用到pip install cftime裝完后在python交互式環(huán)境里驗(yàn)證一下各庫版本能否配合。實(shí)際工作中遇到過xarray升級到新版后和高版本cartopy在投影計(jì)算上的兼容性問題表現(xiàn)是運(yùn)行不報(bào)錯(cuò)但出圖時(shí)地圖邊界是空的這種黑匣子問題最耗時(shí)間。建議用一個(gè)固定版本的虛擬環(huán)境跑完整套流程別頻繁升級。cftime這個(gè)庫單獨(dú)說明一下CF約定時(shí)間坐標(biāo)是氣象數(shù)據(jù)的標(biāo)準(zhǔn)寫法但pandas原生不認(rèn)識(shí)它需要cftime做轉(zhuǎn)換層后面在避坑章節(jié)會(huì)展開講。2.2 NetCDF讀取三步維度、變量、屬性逐項(xiàng)確認(rèn)拿到任何一份氣象數(shù)據(jù)文件不要急著轉(zhuǎn)DataFrame也不要急著畫圖。先做三個(gè)確認(rèn)動(dòng)作看維度、看變量、看全局屬性。這三個(gè)動(dòng)作能省掉后面90%的排錯(cuò)時(shí)間尤其是當(dāng)你接手的是別人傳過來的數(shù)據(jù)集時(shí)。import xarray as xr ds xr.open_dataset(temp_daily_2023.nc) print(ds.dims) # 維度信息time/lat/lon分別多長 print(ds.data_vars) # 變量列表具體有哪些物理量 print(ds.attrs) # 全局屬性數(shù)據(jù)來源、單位、參考時(shí)間這段代碼的邏輯是先把數(shù)據(jù)集對象打開然后分別打印三部分元數(shù)據(jù)。dims告訴你每個(gè)維度的長度比如time是365、lat是181、lon是360說明這是一份全球1度分辨率的逐日數(shù)據(jù)data_vars列出變量常見的有temp、precip、pressure也可能叫t2m、tp這類短名需要查屬性確認(rèn)attrs里最關(guān)鍵的是單位說明氣溫是開爾文還是攝氏度降水是毫米還是千克每平方米差一個(gè)數(shù)量級圖就全錯(cuò)了。參數(shù)說明上xr.open_dataset的第一個(gè)參數(shù)是文件路徑還可以加engine參數(shù)顯式指定讀取引擎比如enginenetcdf4。如果文件特別大比如幾十GB的再分析資料建議加上chunks參數(shù)做懶加載這一步對機(jī)器的內(nèi)存壓力會(huì)小很多。ds xr.open_dataset(temp_daily_2023.nc, chunks{time: 30})加了chunks之后文件不會(huì)一次性讀進(jìn)內(nèi)存而是按time維度每30天一個(gè)塊真正計(jì)算時(shí)才去讀。處理ERA5這類再分析數(shù)據(jù)時(shí)這個(gè)參數(shù)幾乎是必加的否則16GB內(nèi)存的機(jī)器根本扛不住全球逐小時(shí)數(shù)據(jù)。2.3 數(shù)據(jù)清洗邊界缺測值、時(shí)間戳和經(jīng)緯度單位讀取流程走完后下一步是清洗。氣象數(shù)據(jù)的臟和業(yè)務(wù)數(shù)據(jù)的臟不太一樣常見的有三種情況。第一種是缺測值占位缺測值通常用9999、1.0e36這類極端值填充如果沒處理就參與統(tǒng)計(jì)均值和方差會(huì)被拉到一個(gè)離譜的數(shù)值上。第二種是時(shí)間戳格式NetCDF里的time維度經(jīng)常是hours since 1900-01-01 00:00:00這種參考時(shí)間格式需要轉(zhuǎn)成datetime類型才能做按月聚合。第三種是經(jīng)緯度坐標(biāo)單位有些數(shù)據(jù)集的lat/lon單位是弧度而不是度不轉(zhuǎn)換直接切片會(huì)切出不存在的位置。# 處理缺測值占位 ds ds.where(ds[temp] ! -9999) # 統(tǒng)一時(shí)間坐標(biāo)為datetime類型 ds xr.decode_cf(ds) # 確認(rèn)經(jīng)緯度范圍 print(float(ds.lat.min()), float(ds.lat.max()))處理缺測值時(shí)where會(huì)把等于9999的位置置為NaN后續(xù)計(jì)算默認(rèn)會(huì)跳過NaN這比手動(dòng)填0要安全得多。decode_cf是把CF約定的時(shí)間編碼轉(zhuǎn)成可讀的datetime轉(zhuǎn)完后可以直接用ds[temp].sel(time2023-01-15)這種寫法取數(shù)。經(jīng)緯度范圍打印一下如果min和max是0.5到1.5這種量級說明單位是弧度需要乘以180再除以圓周率換算回度。清洗這一步看起來瑣碎但它決定了后面所有圖會(huì)不會(huì)翻車。常見的情況是數(shù)據(jù)源說明里寫著溫度單位是K實(shí)際文件里存的是攝氏度畫出來的全國7月平均氣溫只有5度色標(biāo)全藍(lán)錯(cuò)得離譜。所以每次拿到新數(shù)據(jù)集我習(xí)慣先打印一個(gè)具體點(diǎn)的數(shù)值和常識(shí)對一下數(shù)量級再往下走。3. 核心分析氣溫、降水、氣壓的時(shí)間序列與空間場計(jì)算環(huán)境就緒、數(shù)據(jù)干凈之后才進(jìn)入真正意義上的氣象分析。這一章圍繞三個(gè)最常見訴求展開時(shí)間序列上算月平均和距平空間場上做區(qū)域平均和季節(jié)平均場以及把極端事件量化成簡單指數(shù)。每一步都給出最小可復(fù)現(xiàn)代碼并說明參數(shù)調(diào)整的邊界。3.1 時(shí)間序列用resample算月平均用groupby算氣候態(tài)距平溫度序列分析最基礎(chǔ)的兩個(gè)操作是降采樣和距平計(jì)算。降采樣把逐日數(shù)據(jù)聚合成逐月減少噪音距平表示某個(gè)時(shí)刻相對氣候態(tài)的偏移是判斷偏暖還是偏冷的核心指標(biāo)python數(shù)據(jù)分析與可視化里這塊也是最常用的。import xarray as xr ds xr.open_dataset(temp_daily_2023.nc) temp ds[temp].where(ds[temp] ! -9999) # 逐日數(shù)據(jù)降采樣為逐月平均 monthly_mean temp.resample(time1ME).mean() # 計(jì)算逐月氣候態(tài)假設(shè)文件內(nèi)含多年數(shù)據(jù)按月份分組平均 climatology temp.groupby(time.month).mean(dimtime) # 距平 當(dāng)月值 - 氣候態(tài) anomaly monthly_mean - climatology理解這段代碼的關(guān)鍵在于resample和groupby兩個(gè)操作的分工。resample是按時(shí)間頻率重采樣1ME表示按月取末值作為分組錨點(diǎn)常用寫法還有1MS表示按月初。mean是聚合函數(shù)表示組內(nèi)平均也可以用max、min換成月最高溫、月最低溫。groupby更靈活它按時(shí)間坐標(biāo)的month分量分組把多年數(shù)據(jù)中所有1月聚到一起求平均得到的就是1月氣候態(tài)dimtime明確表示縮減時(shí)間維度只保留month作為新的分組標(biāo)簽。參數(shù)上要注意的是time.month返回的是1到12的整數(shù)所以climatology的維度是month長度12。如果你手里的數(shù)據(jù)只有單年算不了氣候態(tài)可以退一步用多年歷史數(shù)據(jù)單獨(dú)算一個(gè)climatology.nc文件再和當(dāng)年數(shù)據(jù)做對齊運(yùn)算。這個(gè)對齊需要保證兩個(gè)文件網(wǎng)格分辨率完全一致否則xarray會(huì)按坐標(biāo)自動(dòng)插值這在氣象分析里常常是隱患最好在數(shù)據(jù)準(zhǔn)備階段就把分辨率統(tǒng)一。距平的工程化價(jià)值在于它消除了季節(jié)周期。對比1月和7月的原始溫度沒有意義但對比它們的距平有物理含義。做可視化時(shí)距平圖用的色標(biāo)通常是對稱的藍(lán)紅兩色這種圖的視覺沖擊力比原始溫度圖強(qiáng)很多也是氣象報(bào)告里最常出現(xiàn)的圖表類型。3.2 空間場計(jì)算區(qū)域平均、季節(jié)平均場和多層差值當(dāng)分析目標(biāo)是整個(gè)華北平原今夏平均氣溫比常年高了多少就不能只看單點(diǎn)需要做空間聚合?？臻g聚合的常見做法是先按經(jīng)緯度切片選出區(qū)域再對lat和lon兩個(gè)維度求平均。# 選區(qū)域華北平原大致在35—40N, 110—120E region temp.sel(latslice(35, 40), lonslice(110, 120)) # 區(qū)域時(shí)間序列對空間維求平均 region_series region.mean(dim[lat, lon]) # 多年同一季節(jié)的空間平均場 seasonal_mean temp.sel(timetemp[time.season] JJA).mean(dimtime)sel方法用slice做切片時(shí)要求坐標(biāo)軸是單調(diào)遞增或遞減的如果lat坐標(biāo)在數(shù)據(jù)里是南緯在下、北緯在上的遞增排列直接切片即可如果是幾十年前數(shù)據(jù)的北緯在上排列需要先sortby處理。slice(35,40)表示取35到40度之間的所有格點(diǎn)注意這個(gè)區(qū)間是閉區(qū)間邊界值都會(huì)被包含。mean(dim[lat,lon])的效果是把二維空間壓成一個(gè)點(diǎn)得到一條純粹的時(shí)間序列適合做后續(xù)的趨勢回歸。季節(jié)平均場這里的寫法利用了season這個(gè)派生坐標(biāo)xarray會(huì)自動(dòng)從標(biāo)準(zhǔn)時(shí)間坐標(biāo)生成季節(jié)標(biāo)簽DJF是冬季MAM是春季JJA是夏季SON是秋季。如果需要限定夏季三個(gè)月也可以寫得更明確summer_months temp.sel(timetemp[time.month].isin([6, 7, 8]))這段代碼等價(jià)于先篩選月份再求平均邏輯上更直白適合新手理解?？臻g場分析還有一個(gè)常用操作是沿氣壓層差值比如分析500hPa和850hPa的高度差來推斷冷暖平流做法是用interp對level維度線性插值interp_level temp.interp(level[500, 850])interp是xarray內(nèi)置的線性插值方法目標(biāo)level序列放在列表里插值后對象的維度會(huì)變得和原始level維度一致。要注意的是如果原始數(shù)據(jù)里level坐標(biāo)不是等間距的線性插值的誤差會(huì)放大這時(shí)應(yīng)該考慮用scipy的帶權(quán)插值函數(shù)做更高階處理。這一步在無資料地區(qū)的氣象分析里經(jīng)常單獨(dú)拿出來評估。3.3 簡單氣候指數(shù)把極端高溫/低溫事件量化出來分析做完常規(guī)平均還有一個(gè)高頻需求是把極端事件量化。氣候指數(shù)本質(zhì)上是對原始序列做閾值判斷再統(tǒng)計(jì)頻次或累計(jì)量。這里以高溫日數(shù)和熱浪過程為例這兩個(gè)指標(biāo)在夏季高溫評估和農(nóng)業(yè)災(zāi)害風(fēng)險(xiǎn)里非常常用。import numpy as np # 定義高溫日日最高溫超過35℃數(shù)據(jù)為開爾文時(shí)需換算 threshold 308.15 # 35℃ 308.15K hot_flag (temp threshold).astype(int) # 連續(xù)高溫日數(shù)用累計(jì)計(jì)數(shù)識(shí)別連片時(shí)段 run_length hot_flag.groupby( (hot_flag ! hot_flag.shift(time1)).cumsum() ).cumsum()高溫日數(shù)的定義極其簡單布爾比較加類型轉(zhuǎn)換就夠。關(guān)鍵點(diǎn)是單位換算數(shù)據(jù)是開爾文35℃對應(yīng)308.15K換算錯(cuò)了所有指數(shù)都會(huì)錯(cuò)得離譜。連續(xù)高溫日數(shù)的計(jì)算稍微繞一點(diǎn)先構(gòu)造一個(gè)01標(biāo)志序列再用shift判斷前后兩天是否連續(xù)不連續(xù)則分組計(jì)數(shù)加1最后用cumsum輸出連續(xù)長度。這個(gè)邏輯是氣象統(tǒng)計(jì)里識(shí)別過程事件的通用做法同樣適用于連續(xù)降水日數(shù)、連續(xù)無降水日數(shù)。實(shí)際業(yè)務(wù)中做這類時(shí)間序列分析最容易出問題的不是算法而是邊界條件數(shù)據(jù)第一天的前一天是什么狀態(tài)我們并不知道。所以計(jì)算run_length時(shí)第一天的連續(xù)值天然不可信使用時(shí)要考慮截掉序列首尾各一個(gè)點(diǎn)或用多年數(shù)據(jù)把邊界補(bǔ)齊。這個(gè)坑在后面避坑章節(jié)里會(huì)再提一遍。4. 可視化落地靜態(tài)天氣圖到動(dòng)態(tài)氣象面板分析結(jié)果是數(shù)字可視化才讓它變成能對外溝通的圖。這一章按三個(gè)層次來寫先用cartopy畫一張帶底圖的靜態(tài)空間分布圖再用matplotlib動(dòng)畫把時(shí)間演變串起來最后用Panel搭一個(gè)可拖拽、可切換參數(shù)的交互看板。可視化部分的坑往往不在畫圖本身而在投影和坐標(biāo)轉(zhuǎn)換上這兩個(gè)問題會(huì)在第5章統(tǒng)一排查。4.1 用Cartopy畫一個(gè)帶底圖的溫度空間分布圖單張靜態(tài)圖是最基礎(chǔ)的產(chǎn)出一張合格的空間分布圖需要包含四個(gè)要素地圖底圖、等值線填色、色標(biāo)、投影信息。cartopy負(fù)責(zé)前三者的地理框架matplotlib負(fù)責(zé)渲染細(xì)節(jié)。import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature temp_jan ds[temp].sel(time2023-01-15) fig plt.figure(figsize(10, 6)) ax plt.axes(projectionccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linestyle:) im temp_jan.plot.pcolormesh( axax, transformccrs.PlateCarree(), cmapRdYlBu_r, cbar_kwargs{label: 2m temperature (K)} ) ax.set_extent([70, 140, 15, 55]) ax.set_title(2023-01-15 Surface Temperature) plt.savefig(temp_2023_jan15.png, dpi300, bbox_inchestight)這里最關(guān)鍵的是要理解投影坐標(biāo)和數(shù)據(jù)坐標(biāo)是兩套系統(tǒng)。PlateCarree是最簡單的等距圓柱投影經(jīng)緯度直接映射到平面坐標(biāo)對于中國區(qū)域的1度網(wǎng)格足夠用。transformccrs.PlateCarree()告訴matplotlib數(shù)據(jù)本身是經(jīng)緯度坐標(biāo)畫布ax的投影是另一個(gè)投影時(shí)cartopy會(huì)自動(dòng)做轉(zhuǎn)換。pcolormesh比contourf更適合格點(diǎn)數(shù)據(jù)它不做插值直接按網(wǎng)格畫色塊速度也快。cmap用RdYlBu_r這個(gè)色標(biāo)的紅藍(lán)方向是暖色表示高溫冷色表示低溫和氣象行業(yè)的讀圖習(xí)慣一致。set_extent控制地圖顯示范圍這里的經(jīng)緯度范圍正好覆蓋中國全境。savefig的dpi建議300報(bào)文章或者報(bào)告里的圖這個(gè)清晰度才夠。如果數(shù)據(jù)是規(guī)則1度經(jīng)緯度網(wǎng)格pcolormesh可以直接用如果投影比較復(fù)雜比如極地立體投影數(shù)據(jù)坐標(biāo)轉(zhuǎn)換后會(huì)出現(xiàn)空白區(qū)域處理辦法在第5章的避坑里會(huì)講到。4.2 制作逐日溫度演變動(dòng)畫matplotlib.animation的正確用法靜態(tài)圖只能看一個(gè)時(shí)刻逐日動(dòng)畫才能看到過程的演進(jìn)。制作方法是用matplotlib.animation的FuncAnimation逐幀更新數(shù)據(jù)把每一幀存成PNG再合成視頻。實(shí)際項(xiàng)目中我常用兩種產(chǎn)出方式直接保存mp4或者生成HTML5的動(dòng)畫嵌入網(wǎng)頁里后者在演示時(shí)非常方便。from matplotlib.animation import FuncAnimation temp ds[temp].where(ds[temp] ! -9999) fig plt.figure(figsize(10, 6)) ax plt.axes(projectionccrs.PlateCarree()) ax.add_feature(cfeature.COASTLINE, linewidth0.5) def update(frame): ax.clear() ax.add_feature(cfeature.COASTLINE, linewidth0.5) data temp.isel(timeframe) im data.plot.pcolormesh(axax, transformccrs.PlateCarree(), cmapRdYlBu_r, vmin250, vmax310, add_colorbarFalse) ax.set_title(fTime: {str(data.time.values)[:10]}) return im anim FuncAnimation(fig, update, frameslen(temp.time), interval200) anim.save(temp_2023_animation.gif, writerpillow, dpi120)FuncAnimation的核心是update函數(shù)這個(gè)函數(shù)在每一幀執(zhí)行一次先clear掉上一次的繪圖再重新繪制新一幀的數(shù)據(jù)。這里的注意點(diǎn)有兩個(gè)。一是vmin和vmax手動(dòng)固定不要讓它自動(dòng)縮放否則每幀的色標(biāo)范圍都在變動(dòng)畫看起來會(huì)閃爍二是add_colorbarFalse動(dòng)畫里每一幀都加色標(biāo)會(huì)非常慢顏色圖例統(tǒng)一放在最后單獨(dú)輸出或固定。interval200表示每幀間隔200毫秒即每秒5幀。保存gif用pillow writer保存mp4則用ffmpeg writer需要系統(tǒng)裝有ffmpeg。如果你只需要在Jupyter里看效果不用保存文件update函數(shù)返回圖像對象后用HTML(anim.to_html5_video())渲染即可。還有一個(gè)性能優(yōu)化點(diǎn)值得注意——如果數(shù)據(jù)維度很大每幀都重新繪制整個(gè)網(wǎng)格很慢可以重用上一次的quadmesh對象只在set_array里替換數(shù)據(jù)速度能快3倍以上。4.3 用Panel搭一個(gè)可拖拽的氣象數(shù)據(jù)小看板走到這一步單張圖和動(dòng)畫都能產(chǎn)出了。如果分析結(jié)果要交給非技術(shù)同事看或者要放到科室的大屏上輪播就需要一個(gè)交互式的可視化小面板。這里我會(huì)用Panel框架它和Jupyter綁定得好不需要額外寫前端代碼就能搭出可拖拽、可切換參數(shù)的看板。import panel as pn pn.extension() ds xr.open_dataset(temp_daily_2023.nc) selector pn.widgets.Select(name觀測日期, optionslist(ds.time.values)[:10]) pn.depends(selector.param.value) def plot_map(date): data ds[temp].sel(timedate) fig plt.figure(figsize(8, 5)) ax plt.axes(projectionccrs.PlateCarree()) data.plot.pcolormesh(axax, transformccrs.PlateCarree(), cmapRdYlBu_r) plt.close(fig) return fig layout pn.Column(selector, plot_map) layout.servable()用Panel搭看板的思路是用一個(gè)widget組件綁定日期pn.depends裝飾器讓plot_map函數(shù)在參數(shù)變化時(shí)自動(dòng)重算pn.Column組合成縱向布局。servable方法可以讓這個(gè)面板通過panel serve script.py直接起一個(gè)本地服務(wù)瀏覽器打開就能交互。做交互看板時(shí)要注意一個(gè)性能問題每次拖動(dòng)日期下拉框都會(huì)重新運(yùn)行plot_map如果底層數(shù)據(jù)是懶加載的大型NetCDF建議把需要展示的變量先降采樣到1度或2.5度網(wǎng)格或者用緩存機(jī)制存住最近計(jì)算過的結(jié)果。我還習(xí)慣把色標(biāo)范圍固定理由和動(dòng)畫完全一樣用戶拖動(dòng)時(shí)顏色跳動(dòng)會(huì)讓人誤以為數(shù)據(jù)有問題。最后這個(gè)面板里如果要展示多個(gè)氣象變量可以把Select換成多選或RadioButton然后用循環(huán)生成多個(gè)輸出。項(xiàng)目的落地形態(tài)就是一個(gè)小型氣象可視化工具它的邊界在于交互邏輯不宜做得很復(fù)雜真正復(fù)雜的分析計(jì)算寫在Panel外面。5. 避坑氣象數(shù)據(jù)處理的5個(gè)常見翻車點(diǎn)這一章把前面各環(huán)節(jié)里最容易讓初學(xué)者翻車的五個(gè)場景單獨(dú)拎出來按現(xiàn)象、原因、解決的順序梳理清楚。這些問題我基本都在實(shí)際項(xiàng)目中踩過有些甚至踩了兩遍才總結(jié)出規(guī)律。5.1 cftime和pandas時(shí)間索引打架現(xiàn)象讀取NetCDF后執(zhí)行resample或sel(time2023-01-15)報(bào)錯(cuò)提示time坐標(biāo)不是datetime類型或者index類型不匹配。原因CF約定的時(shí)間坐標(biāo)是hours since 1900-01-01 00:00:00xarray解碼后生成的是cftime.DatetimeNoLeap類型不是標(biāo)準(zhǔn)的pandas.DatetimeIndex。當(dāng)數(shù)據(jù)里包含非標(biāo)準(zhǔn)日歷如無閏年日歷時(shí)xarray不會(huì)自動(dòng)轉(zhuǎn)成pandas類型直接用字符串切片就會(huì)失敗。解決讀取時(shí)顯式使用xr.decode_cf(ds)如果日歷仍是非標(biāo)準(zhǔn)類型就把time坐標(biāo)手動(dòng)轉(zhuǎn)換后寫回。轉(zhuǎn)換用ds[time] ds.indexes[time].to_datetimeindex()。如果數(shù)據(jù)量太大建議在清洗階段一次性轉(zhuǎn)換并保存成一個(gè)新的nc文件后續(xù)所有分析都用這個(gè)清洗后的文件避免每個(gè)腳本都要重復(fù)處理。5.2 經(jīng)緯度維度名寫反地圖直接變形現(xiàn)象畫出來的圖經(jīng)緯度范圍完全不對中國區(qū)域的地圖變成橫向拉伸的長條形或者色塊和海岸線對不上。原因不同數(shù)據(jù)源的維度命名不一致有的是lat/lon有的是latitude/longitude還有的是y/x而且數(shù)組的行列順序可能是經(jīng)緯度反的。用plot.pcolormesh時(shí)x和y參數(shù)如果傳反了圖自然就歪。解決先打印ds.dims確認(rèn)維度名和排列順序然后在取數(shù)時(shí)用sel方法按名字切片而非按位置取數(shù)比如sel(latitudeslice(15,55), longitudeslice(70,140))。如果數(shù)據(jù)本身是y/x這種網(wǎng)格坐標(biāo)先重命名維度再處理ds ds.rename({y: lat, x: lon})。養(yǎng)成按維度名操作的習(xí)慣后基本不會(huì)再出現(xiàn)維度寫反的翻車。5.3 缺測值混進(jìn)統(tǒng)計(jì)計(jì)算現(xiàn)象區(qū)域平均結(jié)果出現(xiàn)NaN或者距平圖上一大片區(qū)域變成空白而不是正常的色塊。原因缺測值在數(shù)據(jù)里表現(xiàn)為特殊的_FillValue例如9999如果不做掩膜處理NaN會(huì)通過均值傳播到所有下游計(jì)算。解決在進(jìn)入任何統(tǒng)計(jì)計(jì)算之前先統(tǒng)一處理缺測值。用ds[temp] ds[temp].where(ds[temp] ! -9999)把占位值替換成NaN之后所有聚合函數(shù)默認(rèn)skipnaTrue自動(dòng)跳過NaN。需要特別注意的地方是進(jìn)行空間平均時(shí)如果某個(gè)格點(diǎn)在整段時(shí)間都缺測平均結(jié)果仍會(huì)是NaN這時(shí)可以用region.mean(dim[lat,lon], skipnaTrue)顯式跳過。更穩(wěn)妥的做法是統(tǒng)計(jì)每個(gè)格點(diǎn)的有效觀測數(shù)把低于閾值的格點(diǎn)直接標(biāo)為NaN。5.4 Cartopy投影邊界撕裂現(xiàn)象使用非PlateCarree投影如Albers等積投影或Lambert投影繪圖時(shí)中國東北或南海區(qū)域的海岸線出現(xiàn)橫線撕裂或色塊邊緣參差不齊。原因數(shù)據(jù)坐標(biāo)是經(jīng)緯度網(wǎng)格投影到目標(biāo)投影后原本一條經(jīng)線上連續(xù)的格點(diǎn)被切到了畫布兩側(cè)或某條邊界上cartopy無法自動(dòng)判斷格點(diǎn)的連續(xù)性。解決先用data.plot.pcolormesh(axax, transformccrs.PlateCarree())明確告訴繪圖函數(shù)數(shù)據(jù)原始坐標(biāo)是經(jīng)緯度再在ax上設(shè)置目標(biāo)投影。如果撕裂依然存在把數(shù)據(jù)先插值到目標(biāo)投影的規(guī)則網(wǎng)格上再畫圖。插值用scipy.interpolate.griddata目標(biāo)網(wǎng)格用投影坐標(biāo)下的np.meshgrid生成。對于中國區(qū)域我還習(xí)慣把a(bǔ)x的邊界設(shè)置成圓角矩形既好看也能掩蓋某些邊界處的毛刺。5.5 文件句柄未關(guān)閉內(nèi)存越用越大現(xiàn)象運(yùn)行多個(gè)分析腳本后內(nèi)存占用持續(xù)上升不退出python進(jìn)程就不釋放甚至出現(xiàn)OSError: [Errno 24] Too many open files。原因xr.open_dataset默認(rèn)是懶加載文件句柄不會(huì)立即釋放。在循環(huán)里反復(fù)打開文件卻不關(guān)閉句柄和懶加載數(shù)據(jù)會(huì)累積在內(nèi)存里。尤其在讀取多個(gè)NetCDF文件合成一個(gè)時(shí)間序列時(shí)最容易觸發(fā)。解決養(yǎng)成用上下文管理器的習(xí)慣或者顯式調(diào)用close。with xr.open_dataset(temp_daily_2023.nc) as ds: result ds[temp].mean(dimtime) # 文件在此處自動(dòng)關(guān)閉 # 多個(gè)文件合并時(shí)合并完成后立即釋放句柄 ds_list [xr.open_dataset(f) for f in file_list] combined xr.concat(ds_list, dimtime) for ds in ds_list: ds.close()這兩個(gè)例子覆蓋了日常使用和批處理兩種場景。前者的with塊在代碼退出后自動(dòng)關(guān)閉文件后者的循環(huán)close確保合并操作的臨時(shí)句柄被釋放。如果在實(shí)際運(yùn)行中內(nèi)存還是控制不住就要考慮用open_dataset(..., chunks{})配合dask做分布式處理而不是盲目加內(nèi)存。6. 進(jìn)階把分析流程封裝成可重復(fù)執(zhí)行的小工具走到這一步讀、洗、算、畫四條鏈路你都跑通了。接下來值得做的一件事是把這套流程封裝成一個(gè)命令行可執(zhí)行的Python腳本讓換數(shù)據(jù)源、換年份、換區(qū)域時(shí)不用再改一堆參數(shù)。我一般會(huì)按函數(shù)粒度拆分load_data負(fù)責(zé)讀文件和清洗cal_temp_index負(fù)責(zé)時(shí)間序列和指數(shù)計(jì)算plot_map負(fù)責(zé)畫圖main函數(shù)負(fù)責(zé)串起整個(gè)流程。函數(shù)之間只通過DataArray對象傳參不共享全局變量。這樣做的最大好處是臨時(shí)接到一個(gè)把去年數(shù)據(jù)重新分析一遍的需求時(shí)只需要改main里的文件路徑和一個(gè)年份參數(shù)十分鐘內(nèi)能出全套圖。def main(data_path: str, year: int, region: tuple) - None: ds load_data(data_path) monthly cal_monthly_stats(ds, year) anomaly cal_anomaly(monthly) plot_map(anomaly, region, outputfanomaly_{year}.png)這個(gè)封裝是不是過度設(shè)計(jì)要看場景。如果你只是自己寫一次性的分析腳本函數(shù)拆分反而增加負(fù)擔(dān)但如果要讓同事復(fù)用或者要定時(shí)產(chǎn)出一組圖這個(gè)封裝就是值得的。參數(shù)校驗(yàn)和異常處理建議在main入口做比如文件路徑不存在時(shí)給出明確提示而不是拋一個(gè)看不懂的堆棧。另外一個(gè)我常用的習(xí)慣是畫完圖后把計(jì)算的中間結(jié)果同步導(dǎo)出成CSV方便非Python用戶直接拿Excel做二次分析。所謂可視化大屏最后的落地也往往是這些CSV加PNG再加少量交互組件組合出來的底層的計(jì)算邏輯并不復(fù)雜?；乜催@套流程最讓我記憶深刻的教訓(xùn)是氣象數(shù)據(jù)分析中80%的錯(cuò)誤不是分析本身的問題而是數(shù)據(jù)格式、單位、缺測和投影這些周邊問題。先把數(shù)據(jù)邊界摸清楚分析的信心才有基礎(chǔ)。希望這篇筆記能幫你少走一遍這些彎路。本文還有配套的精品資源點(diǎn)擊獲取