做氣候數(shù)據(jù)分析的同行應(yīng)該都繞不開這樣一批數(shù)據(jù)中國(guó)全國(guó)地面站點(diǎn)1981-2010年月平均氣象數(shù)據(jù)包含月平均降水、月平均相對(duì)濕度、月最大降水、月最高/最低溫度以及月最高溫平均值、月最低溫平均值最后整理成一份全國(guó)高低溫站點(diǎn)基礎(chǔ)數(shù)據(jù)表。這類數(shù)據(jù)是氣候統(tǒng)計(jì)、建筑節(jié)能分區(qū)、農(nóng)業(yè)區(qū)劃和災(zāi)害評(píng)估里最常用的底數(shù)但也是我見過被問得最多、最容易踩坑的數(shù)據(jù)之一。這篇文章不講空泛的大道理只講我實(shí)際整理、清洗、統(tǒng)計(jì)這份數(shù)據(jù)時(shí)遇到的口徑問題、單位坑以及一套可復(fù)用的Python處理流程。無論你是要畫中國(guó)區(qū)域的氣候平均圖還是要給某個(gè)行業(yè)項(xiàng)目提取站點(diǎn)氣候指標(biāo)都可以照著這個(gè)思路往下走至少能少走幾個(gè)月的彎路。1. 這批數(shù)據(jù)到底是什么幾個(gè)要素名稱不能想當(dāng)然1.1 1981-2010為何被當(dāng)成“標(biāo)準(zhǔn)氣候期”先說時(shí)間范圍。1981-2010這30年是氣象行業(yè)常用的一個(gè)“標(biāo)準(zhǔn)氣候期”。世界氣象組織建議用連續(xù)30年的觀測(cè)記錄來計(jì)算氣候平均值目的是把年際波動(dòng)盡量平滑掉讓“某地夏季多熱”“冬季多冷”這種描述有一個(gè)穩(wěn)定的參照系。每過10年這個(gè)基準(zhǔn)會(huì)整體向后滑動(dòng)一次所以早期研究里常見1961-1990最近這些年很多新標(biāo)準(zhǔn)已經(jīng)切到1991-2020但1981-2010仍然大量出現(xiàn)在建筑規(guī)范、能源設(shè)計(jì)參數(shù)和歷史氣候?qū)Ρ戎?。我拿到這類數(shù)據(jù)時(shí)第一件事永遠(yuǎn)是確認(rèn)時(shí)間粒度標(biāo)題寫著1981-2010但到底給的是“1981年1月到2010年12月的逐月實(shí)測(cè)記錄”還是“1981-2010年每個(gè)月多年平均后的氣候態(tài)數(shù)據(jù)”這兩種數(shù)據(jù)長(zhǎng)得像用法完全不一樣。前者是時(shí)間序列可以做年際變化、趨勢(shì)分析后者是統(tǒng)計(jì)產(chǎn)品只能當(dāng)基準(zhǔn)值用。不少項(xiàng)目沖突就是這么來的兩邊都在說“月平均降水”一個(gè)說的是某個(gè)具體月份的實(shí)測(cè)降水另一個(gè)說的是30年同月平均降水對(duì)數(shù)必然對(duì)不上。1.2 月平/極值/平均高低溫同一個(gè)“溫度”有四種讀法標(biāo)題里“月最高溫度最低溫度”和“月最高溫平均值最低溫平均值”這兩串詞是最容易歧義的地方。按我接觸月值數(shù)據(jù)集的習(xí)慣需要拆成四個(gè)概念月平均氣溫某個(gè)站點(diǎn)當(dāng)月所有日平均氣溫再求平均代表這個(gè)月整體的冷熱水平。月平均最高氣溫當(dāng)月每天日最高氣溫的平均值標(biāo)題里說的“月最高溫平均值”大概率指它反映白天的一般高溫水平。月平均最低氣溫當(dāng)月每天日最低氣溫的平均值對(duì)應(yīng)“月最低溫平均值”反映夜間的一般低溫水平。月極端最高/最低氣溫當(dāng)月出現(xiàn)過的最高、最低瞬時(shí)氣溫對(duì)應(yīng)標(biāo)題里“月最高溫度最低溫度”。前兩個(gè)是“平均態(tài)”后兩個(gè)是“極值”完全不是一回事。做站點(diǎn)排序時(shí)要特別小心如果要把“全國(guó)最熱站點(diǎn)”排出來到底比的是月平均最高氣溫還是月極端最高氣溫結(jié)果會(huì)差很多。吐魯番的極端高溫可以接近50但它的月平均最高氣溫通常沒有這么夸張反過來一個(gè)夏天氣溫起伏不大的海濱城市月平均最高氣溫可能不低但極端最高溫度反而排不到前面去。這類差別如果不寫清楚最終報(bào)告是站不住腳的。還有“月最大降水”這個(gè)概念。它通常指當(dāng)月內(nèi)某一天出現(xiàn)的最大日降水量而不是當(dāng)月降水總量。比如某站7月降水總量300毫米其中最大的一場(chǎng)暴雨下了120毫米那“月最大降水”是120毫米不是300毫米。標(biāo)題里的“月平均降水”更接近多年平均后的月降水量或日平均降水量使用時(shí)務(wù)必看原始字段說明。2. 數(shù)據(jù)預(yù)處理把原始月值文件收拾干凈2.1 讀取站點(diǎn)信息和月值數(shù)據(jù)絕大多數(shù)地面站點(diǎn)月值數(shù)據(jù)的原始組織方式都很“古早”。站點(diǎn)基礎(chǔ)信息一般是單獨(dú)一張表包含區(qū)站號(hào)、緯度、經(jīng)度、海拔高度月值數(shù)據(jù)可能是按要素分開的多列寬表也可能是長(zhǎng)表一行代表一個(gè)站點(diǎn)某年某月后面跟著氣溫、降水、濕度等若干列。我的習(xí)慣是先把所有內(nèi)容讀成長(zhǎng)表統(tǒng)一成四列核心主鍵站點(diǎn)代號(hào)、年、月、要素名。站點(diǎn)代號(hào)一定要用字符串讀取不能用整數(shù)否則“53068”這種編號(hào)開頭是0的站號(hào)會(huì)丟失前導(dǎo)零。文件編碼也得注意早期數(shù)據(jù)很多用GBK直接用pandas默認(rèn)編碼讀會(huì)報(bào)亂碼我一般這樣起步import pandas as pd # 站點(diǎn)元數(shù)據(jù)站號(hào)、緯度、經(jīng)度、海拔 station pd.read_csv( station_info.csv, dtype{station_id: str}, encodingutf-8 ) # 月值數(shù)據(jù)先用文本方式讀便于排查分隔符問題 df pd.read_csv( monthly_raw.csv, dtype{station_id: str}, encodinggbk, na_values[32744, 99999, 999999] )這里把異常缺測(cè)碼直接放進(jìn)na_values是為了讓數(shù)據(jù)在后續(xù)聚合統(tǒng)計(jì)中自然跳過而不是被當(dāng)成有效值參與平均。不同數(shù)據(jù)源缺測(cè)碼不一樣有的用32744表示“無觀測(cè)”有的用32766表示“降水微量”還有的用一系列很怪的數(shù)字所以強(qiáng)烈建議拿到數(shù)據(jù)后先統(tǒng)計(jì)每個(gè)字段的最小值、最大值和頻數(shù)分布再?zèng)Q定缺測(cè)碼清單。2.2 缺失值、異常值與單位換算清洗時(shí)最大也最隱蔽的坑是單位。月值數(shù)據(jù)里的氣溫常用0.1攝氏度存儲(chǔ)降水量常用0.1毫米存儲(chǔ)不知道的會(huì)直接拿原始值去分析和畫圖最后得到的氣溫全部是幾十幾百攝氏度降水則整體放大10倍。這個(gè)鍋不怪?jǐn)?shù)據(jù)怪?jǐn)?shù)據(jù)字典沒讀。我在處理任何字段前都會(huì)先查一遍字段表哪怕數(shù)據(jù)說明里寫著“氣溫單位0.1℃”也要用真實(shí)樣本驗(yàn)證一下# 假設(shè)原始溫度列叫 temp_raw降水列叫 prec_raw # 先看原始統(tǒng)計(jì)量判斷是不是0.1倍關(guān)系 print(df[[temp_raw, prec_raw]].describe()) # 確認(rèn)后統(tǒng)一換算成常規(guī)單位 df[temp] df[temp_raw] / 10.0 # 單位攝氏度 df[prec] df[prec_raw] / 10.0 # 單位毫米還要處理經(jīng)緯度格式。有些數(shù)據(jù)給的是度分秒比如“09545”表示95度45分有些給的是十進(jìn)制度比如“95.75”。如果直接按數(shù)值畫圖站點(diǎn)能偏出上千公里。我一般寫一個(gè)簡(jiǎn)單函數(shù)把度分秒轉(zhuǎn)成十進(jìn)制度再畫散點(diǎn)圖看站點(diǎn)分布是否落在中國(guó)版圖范圍內(nèi)。異常值檢測(cè)也不能只看“低于-60度就刪”。中國(guó)冬季站點(diǎn)出現(xiàn)過-52.3攝氏度的極端低溫漠河夏季新疆和吐魯番一帶超過50度也是真實(shí)記錄。一上來就按常識(shí)砍極端值等于把最有價(jià)值的極值刪了。正確做法是先換算單位再按每個(gè)要素做箱線圖和月度分布結(jié)合氣象學(xué)常識(shí)判斷哪些值屬于“可疑”。比如降水量為負(fù)值相對(duì)濕度大于100或小于0這些可以直接判異常氣溫極高值先和鄰近站點(diǎn)對(duì)比如果只有一個(gè)孤立高值且周圍沒有類似現(xiàn)象再考慮是否測(cè)站故障。3. 統(tǒng)計(jì)口徑拆解平均和極值的算法別混用3.1 月平均氣溫與月平均最高/最低氣溫很多剛接觸數(shù)據(jù)的人以為“月平均氣溫 月平均最高氣溫加月平均最低氣溫除以2”這是近似不是定義。真正的月平均氣溫是把當(dāng)月每一天的日平均氣溫都加一起再除以當(dāng)月天數(shù)。日平均氣溫的算法在不同觀測(cè)標(biāo)準(zhǔn)下還不完全一樣有的用一天24小時(shí)整點(diǎn)氣溫平均有的用“最高加最低再除以2”估算但作為月值產(chǎn)品數(shù)據(jù)生產(chǎn)者已經(jīng)把這些細(xì)節(jié)處理好了。拿到月值數(shù)據(jù)集后直接使用“月平均氣溫”字段就好不要自己拿最高最低去拼?！霸缕骄罡邭鉁亍焙汀霸缕骄畹蜌鉁亍眲t非常直觀把每天的日最高氣溫求月平均把每天的日最低氣溫求月平均。它們之間的差值就是學(xué)術(shù)界和工程界常說的“氣溫日較差”。建筑熱工設(shè)計(jì)里日較差影響圍護(hù)結(jié)構(gòu)蓄熱策略農(nóng)業(yè)上日較差大有利于糖分積累。所以這兩個(gè)字段不是冗余它們至少有獨(dú)立的物理意義。我在寫氣候統(tǒng)計(jì)報(bào)告時(shí)會(huì)把這三列同時(shí)輸出月平均氣溫、月平均最高氣溫、月平均最低氣溫。缺少任何一個(gè)后續(xù)計(jì)算都容易被迫用近似公式。3.2 月最大降水與月極端溫度的屬性“月最大降水”和“月降水量”的區(qū)別前面已經(jīng)提過但實(shí)際檢查時(shí)還是經(jīng)常出問題。一個(gè)站月降水量可能是200毫米月最大降水可能是80毫米這表示整月下了幾天雨其中最強(qiáng)的一天貢獻(xiàn)了80毫米。如果只想算月度水資源總量用月降水量如果想分析短時(shí)強(qiáng)降水和暴雨過程要關(guān)注月最大降水。兩者在極端氣候評(píng)估里都有價(jià)值但不能互相替代。“月極端最高氣溫”和“月極端最低氣溫”同樣如此。它們記錄的是某個(gè)站點(diǎn)在這個(gè)月里抓住的最極端值可能只存在了幾分鐘也可能持續(xù)了幾個(gè)小時(shí)。像“全國(guó)6月最高氣溫排行榜”這種統(tǒng)計(jì)本質(zhì)上排的是月極端最高氣溫而“夏季體感熱不熱”應(yīng)該看月平均最高氣溫。把這兩類指標(biāo)混著用會(huì)得到很誤導(dǎo)人的結(jié)論。3.3 多年平均怎么算才合理如果目標(biāo)是生成1981-2010年的氣候態(tài)正確算法是“同月求平均”而不是“全部混在一起求平均”。也就是說要先把每年1月的月平均氣溫單獨(dú)取出30個(gè)值求平均得到1月氣候平均氣溫每年2月同樣處理得到2月氣候平均氣溫以此類推。降水和濕度也是一樣。這里有三個(gè)容易出錯(cuò)的地方。第一如果某站某年某月缺測(cè)直接把這一個(gè)月剔除再算剩余年份的平均值而不是當(dāng)作0值參與平均。第二要考慮樣本量。一個(gè)站點(diǎn)如果30年里只有3年有數(shù)據(jù)算出來的“氣候平均”代表性很差應(yīng)該加篩選條件。第三極端溫度的多年平均有特殊含義比如月極端最高氣溫的30年平均意思是“每年這個(gè)月里最高氣溫的典型水平”它會(huì)比單一極端值更穩(wěn)健但仍不等于月平均最高氣溫。我在下面的實(shí)操里會(huì)給出具體代碼這個(gè)步驟用pandas很容易實(shí)現(xiàn)難的是意識(shí)到口徑差異。4. 實(shí)操從1981-2010月值數(shù)據(jù)到氣候統(tǒng)計(jì)產(chǎn)品4.1 讀取轉(zhuǎn)換與逐站統(tǒng)計(jì)的代碼骨架假設(shè)我已經(jīng)把數(shù)據(jù)清洗成一張長(zhǎng)表字段包括station_id、year、month、tavg、tmax_avg、tmin_avg、tmax_extreme、tmin_extreme、prec、prec_max、rh。下面這個(gè)流程可以計(jì)算每個(gè)站點(diǎn)在1981-2010年間的逐月氣候平均值# 篩選研究時(shí)段 period df[(df[year] 1981) (df[year] 2010)].copy() # 先統(tǒng)計(jì)每個(gè)站點(diǎn)每個(gè)月份有多少個(gè)有效年份 def monthly_climate(sub): # 這里對(duì)同一個(gè)月做平均且只保留至少25年有效數(shù)據(jù)的月份 months [] for m in range(1, 13): month_data sub[sub[month] m] if len(month_data.dropna(subset[tavg])) 25: continue months.append({ month: m, tavg: month_data[tavg].mean(), tmax_avg: month_data[tmax_avg].mean(), tmin_avg: month_data[tmin_avg].mean(), tmax_extreme: month_data[tmax_extreme].mean(), tmin_extreme: month_data[tmin_extreme].mean(), prec: month_data[prec].mean(), prec_max: month_data[prec_max].mean(), rh: month_data[rh].mean(), valid_years: len(month_data.dropna(subset[tavg])), }) return pd.DataFrame(months) result period.groupby(station_id, group_keysFalse).apply(monthly_climate).reset_index()這段代碼的核心思想是“先分站點(diǎn)再分月份最后多年平均”。valid_years字段很重要后面做站點(diǎn)篩選時(shí)可以直接用。25年這個(gè)閾值是我常用的經(jīng)驗(yàn)值既考慮早期數(shù)據(jù)缺測(cè)較多的現(xiàn)實(shí)又保證平均值有足夠的代表性。如果你做的是嚴(yán)格的氣候標(biāo)準(zhǔn)值建議直接要求30年完整缺一年都不行如果只是做區(qū)域背景分析25年也可以接受。4.2 有效站點(diǎn)篩選與高低溫站點(diǎn)基礎(chǔ)數(shù)據(jù)構(gòu)建“高低溫站點(diǎn)基礎(chǔ)數(shù)據(jù)”本質(zhì)上要回答兩個(gè)問題每個(gè)站點(diǎn)的正常冷暖水平是多少以及歷史極值是多少。前者看tmax_avg和tmin_avg后者看tmax_extreme和tmin_extreme。我會(huì)再做一張年尺度匯總表把12個(gè)月平均成年值annual result.groupby(station_id).agg( tavg_annual(tavg, mean), tmax_avg_annual(tmax_avg, mean), tmin_avg_annual(tmin_avg, mean), tmax_extreme_annual(tmax_extreme, max), tmin_extreme_annual(tmin_extreme, min), prec_annual(prec, sum), ).reset_index()注意年降水不能直接對(duì)“月平均降水”求平均再乘12必須用原始月降水量在30年同月平均后再累加。我這里只是演示字段關(guān)系實(shí)際項(xiàng)目里要回到period表先算各月降水量的30年平均再把12個(gè)月加起來。為了防止站點(diǎn)在30年間遷站后經(jīng)緯度變化還要和station表做一次關(guān)聯(lián)每個(gè)站點(diǎn)只保留一個(gè)最終坐標(biāo)。做完匯總后我通常輸出三張表站點(diǎn)基礎(chǔ)信息、逐月氣候平均值、年尺度極端值排名。這就是很多項(xiàng)目里說的“站點(diǎn)基礎(chǔ)數(shù)據(jù)”產(chǎn)品。用CSV保存會(huì)被人反復(fù)索要用Parquet保存更高效我會(huì)兩個(gè)都輸出annual.merge(station, onstation_id, howleft).to_csv(station_climate_1981_2010.csv, indexFalse) annual.merge(station, onstation_id, howleft).to_parquet(station_climate_1981_2010.parquet)4.3 把結(jié)果落到地圖上統(tǒng)計(jì)完表格后很多人下一步就是畫站點(diǎn)分布圖。我不建議一上來就做空間插值先畫散點(diǎn)圖看數(shù)據(jù)是否合理是最快的質(zhì)檢方式import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(8, 7)) sc ax.scatter( annual[lon], annual[lat], cannual[tavg_annual], cmapcoolwarm, s10, linewidths0 ) plt.colorbar(sc, labelAnnual mean temperature (C)) ax.set_xlabel(Longitude) ax.set_ylabel(Latitude) plt.savefig(station_tavg_map.png, dpi300)如果畫出來的站點(diǎn)點(diǎn)位明顯有聚集、有空缺且空缺不是自然區(qū)域比如青藏高原西部站點(diǎn)稀疏是正常的那就要回頭查數(shù)據(jù)。散點(diǎn)圖能發(fā)現(xiàn)經(jīng)緯度格式錯(cuò)誤、站號(hào)合并錯(cuò)誤、坐標(biāo)缺失等問題。更嚴(yán)格的做法是加載省級(jí)行政區(qū)邊界檢查站點(diǎn)是否落在陸地上但日常調(diào)試用散點(diǎn)圖已經(jīng)能攔住大部分坑。5. 實(shí)際項(xiàng)目中常見的問題與排查實(shí)錄5.1 最常翻車的三類問題第一類單位或缺測(cè)碼沒有處理就進(jìn)入統(tǒng)計(jì)。癥狀是平均氣溫高得離譜或者降水總量全部放大10倍。第二類把“月極端最高氣溫”當(dāng)成“月平均最高氣溫”使用導(dǎo)致部分內(nèi)陸站點(diǎn)“平均高溫”排名異?？壳?。第三類站點(diǎn)合并時(shí)忽略遷移記錄使同一站號(hào)在30年間經(jīng)緯度跳了好幾百公里畫出來的軌跡像一條斷線。這三類問題都不是代碼復(fù)雜而是對(duì)數(shù)據(jù)背景不夠敏感。我一般會(huì)在項(xiàng)目開始時(shí)建一個(gè)“數(shù)據(jù)字典”文檔把每個(gè)字段的名稱、單位、有效范圍、缺失碼、來源說明寫清楚后續(xù)所有腳本都依賴這個(gè)字典生成而不是每次手動(dòng)核對(duì)。5.2 一個(gè)單位換算錯(cuò)誤引發(fā)的排查過程有一年我做一個(gè)區(qū)域降水分析讀入某省14個(gè)站點(diǎn)的月值數(shù)據(jù)直接畫出年降水量圖。結(jié)果圖上一片大紅色好幾個(gè)站年降水量超過3萬毫米明顯不可能是現(xiàn)實(shí)記錄。我當(dāng)時(shí)第一反應(yīng)是數(shù)據(jù)本身異常后來用describe查看prec_raw的統(tǒng)計(jì)量發(fā)現(xiàn)最大值是32767像是某個(gè)上下限邊界。進(jìn)一步查詢字段說明才發(fā)現(xiàn)降水量原始單位是0.1毫米我漏了除以10。這個(gè)錯(cuò)誤很好排查但也說明一個(gè)道理任何行業(yè)數(shù)據(jù)單位轉(zhuǎn)換都不能靠“默認(rèn)”。我會(huì)在統(tǒng)計(jì)前后各加一個(gè)“單位校驗(yàn)”步驟比如對(duì)比站點(diǎn)的多年平均年降水量是否落在200-2500毫米這個(gè)合理區(qū)間。如果跑出5000毫米以上或50毫米以下先檢查換算因子再檢查缺測(cè)處理最后才懷疑極端氣候事件。5.3 數(shù)據(jù)質(zhì)量快速檢查清單下面這張表是我每次處理完數(shù)據(jù)后必過的檢查項(xiàng)效率很高檢查項(xiàng)判斷標(biāo)準(zhǔn)常見問題站號(hào)完整性站號(hào)不重復(fù)且均為字符串前導(dǎo)零丟失導(dǎo)致部分站點(diǎn)合并經(jīng)緯度范圍經(jīng)度73-135緯度18-54度分秒格式未轉(zhuǎn)十進(jìn)制度氣溫單位月平均氣溫在-40到35之間0.1℃未除以10降水量單位年降水量大致在50-3000毫米0.1mm未除以10相對(duì)濕度0-100%原數(shù)據(jù)可能是整數(shù)百分比月份完整性每站每年最多12條記錄重復(fù)月份、跨年篩選錯(cuò)誤站點(diǎn)有效性有效年數(shù)不少于25年缺測(cè)年誤當(dāng)0值參與平均平均與極端關(guān)系月平均最高氣溫大于月平均氣溫字段選擇錯(cuò)誤這張清單看著簡(jiǎn)單實(shí)際上每一條都是從項(xiàng)目加班里熬出來的。把它們寫進(jìn)自動(dòng)化測(cè)試后我的數(shù)據(jù)返工率明顯下降。6. 最后分享一點(diǎn)個(gè)人維護(hù)這類數(shù)據(jù)的習(xí)慣我自己的體會(huì)是氣候數(shù)據(jù)項(xiàng)目最值錢的從來不是那張最終結(jié)果表而是你處理數(shù)據(jù)的過程和口徑說明。現(xiàn)在我每做一個(gè)這樣的項(xiàng)目都會(huì)在輸出目錄下放三個(gè)文件原始數(shù)據(jù)不做任何修改的備份、已經(jīng)清洗過的長(zhǎng)表、一份幾百字的數(shù)據(jù)字典。數(shù)據(jù)字典里會(huì)寫清楚“月最大降水”用的是不是日最大值“月平均最高氣溫”有沒有做過0.1倍換算哪些站點(diǎn)因?yàn)橛行陻?shù)不足被剔除。這樣做的好處是幾個(gè)月后別人拿著結(jié)果問我“為什么某站和某站數(shù)據(jù)對(duì)不上”我可以直接翻出字典一分鐘定位到是統(tǒng)計(jì)口徑還是缺測(cè)處理的問題。如果你也準(zhǔn)備整理這批1981-2010地面站點(diǎn)數(shù)據(jù)建議從一開始就保持這個(gè)習(xí)慣。下次再遇到類似需求你就能直接復(fù)用上一輪的腳本和字典把精力放在真正需要判斷的業(yè)務(wù)問題上。