工作需要用到我國(guó)環(huán)境監(jiān)測(cè)站點(diǎn)的污染物濃度數(shù)據(jù)(感謝大佬的分享), 數(shù)據(jù)很全, 不過(guò)csv格式在分析的時(shí)候尤其大量數(shù)據(jù)分析的時(shí)候并不友好, 所以一般要二次處理一下
上學(xué)的時(shí)候搞過(guò)一次, 當(dāng)時(shí)為了查詢數(shù)據(jù)方便, 塞到了sqlite3的單文件數(shù)據(jù)庫(kù)里, 不過(guò)制作起來(lái)很慢(可能是我沒(méi)用并行支持好的數(shù)據(jù)庫(kù)), 近期更新了一下數(shù)據(jù)且又重新調(diào)整了一下數(shù)據(jù)處理邏輯, 這里記錄一下
主要的思路是以netcdf保存(此前做出來(lái)db之后還要再轉(zhuǎn)nc, 麻煩的很), 然后數(shù)據(jù)規(guī)模的話, 物種一共15個(gè), 時(shí)間是2014年到現(xiàn)在逐小時(shí), 站點(diǎn)大概在1000多
采用每月一個(gè)nc文件存儲(chǔ), 每個(gè)物種單獨(dú)變量, 每個(gè)變量的維度是(hour, sites)或者是(hour, cities), 因?yàn)閿?shù)據(jù)除了站點(diǎn)的還有城市的值
然后希望處理流程盡量通用一些, 以后一些其他的觀測(cè)數(shù)據(jù)也基于類似的邏輯和形式去做
那么處理流程大概是三步走
1. 制作站點(diǎn)列表
現(xiàn)有的站點(diǎn)信息文件有點(diǎn)亂, 需要制作一份完整的站點(diǎn)列表信息, 并且最好涵蓋地理位置信息和行政區(qū)信息, 以便處理的時(shí)候篩選用, 如下 (risk后續(xù)討論):
code,name,city,prov,nation,lat,lon,risk
1013A,市監(jiān)測(cè)中心(停運(yùn)150702),天津市,天津市,中國(guó),39.097,117.151,1
1014A,南口路(停運(yùn)161122),天津市,天津市,中國(guó),39.173,117.193,0
在制作過(guò)程中大概會(huì)遇到一些問(wèn)題:
- 數(shù)據(jù)中的站點(diǎn)編號(hào)對(duì)應(yīng)的名字會(huì)變, 這里以新名字為準(zhǔn)
- 數(shù)據(jù)中的站點(diǎn)編號(hào)對(duì)應(yīng)的經(jīng)緯度會(huì)變, 這里以新的經(jīng)緯為準(zhǔn), 但是不知道是什么原因, 是站點(diǎn)搬遷了還是老的定位不準(zhǔn)咋的, 所以基于經(jīng)緯度變化幅度(0~0.1度), 給risk(風(fēng)險(xiǎn)等級(jí))設(shè)定2-11
- 數(shù)據(jù)中站點(diǎn)的經(jīng)緯度信息缺失, 根據(jù)名稱自己在經(jīng)緯度地圖上找的點(diǎn), risk等級(jí)設(shè)定12
- 數(shù)據(jù)中站點(diǎn)編號(hào)對(duì)應(yīng)的我國(guó)行政區(qū)信息不統(tǒng)一, 這個(gè)比較麻煩, 有時(shí)候是叫法問(wèn)題, 有時(shí)候是到不同級(jí)別的問(wèn)題, 有時(shí)候甚至是行政區(qū)變遷的問(wèn)題, 為了保持統(tǒng)一, 我直接根據(jù)最新的shapefile地圖配合經(jīng)緯度信息標(biāo)注站點(diǎn)的行政區(qū)信息, 具體到地級(jí)市/自治州/直轄市
對(duì)于risk問(wèn)題, 穩(wěn)定的站點(diǎn)是0, 僅改名字也是0, 經(jīng)緯度信息有/無(wú)變化的設(shè)為1, 2-12如上, 目的主要是在處理的時(shí)候可以做個(gè)篩選, risk越高, 觀測(cè)結(jié)果和地理信息對(duì)應(yīng)偏差的可能性越高, 可以考慮在對(duì)比模式驗(yàn)證的時(shí)候舍去試試
2. 制作nc數(shù)據(jù)文件
由于要逐月去生成, 且為了同步初始化數(shù)據(jù)和后續(xù)更新數(shù)據(jù)的邏輯, 采用的方式是先init好這個(gè)月的文件, 然后就都是根據(jù)csv文件去update了
注意的點(diǎn)要是:
- 每個(gè)月的文件要以這個(gè)月的sites為準(zhǔn), 不要用上一步制作的全部sites, 不過(guò)站點(diǎn)信息可以從上一步的結(jié)果里獲取, 這一步用wsl跑個(gè)bash腳本來(lái)實(shí)現(xiàn)
- sites和cities的結(jié)果放一個(gè)文件里
- 數(shù)據(jù)較多, 所以基于MS-MPI和mpi4py做了個(gè)并行支持, 單機(jī)4核一起跑也能快很多! 不用py自帶的多進(jìn)程是因?yàn)閙pi更熟悉...而且顯然更好使
最終得到的結(jié)果如下, 包含了站點(diǎn)信息和觀測(cè)數(shù)值
netcdf CNMEE.aqo.site-city.201808 {
dimensions:
site = 1605 ;
city = 367 ;
ymdh = 744 ;
variables:
string site(site) ;
string city(city) ;
string ymdh(ymdh) ;
string site_name(site) ;
string site_city(site) ;
string site_prov(site) ;
string site_nation(site) ;
double site_lat(site) ;
site_lat:_FillValue = NaN ;
double site_lon(site) ;
site_lon:_FillValue = NaN ;
int64 site_risk(site) ;
float AQI(ymdh, site) ;
float PM2.5(ymdh, site) ;
PM2.5:units = "ug/m3" ;
float PM2.5_24h(ymdh, site) ;
PM2.5_24h:units = "ug/m3" ;
float PM10(ymdh, site) ;
PM10:units = "ug/m3" ;
... ...
float CO(ymdh, site) ;
CO:units = "mg/m3" ;
float CO_24h(ymdh, site) ;
CO_24h:units = "mg/m3" ;
float AQI_city(ymdh, city) ;
float PM2.5_city(ymdh, city) ;
PM2.5_city:units = "ug/m3" ;
float PM2.5_24h_city(ymdh, city) ;
PM2.5_24h_city:units = "ug/m3" ;
... ...
// global attributes:
:source = "https://quotsoft.net/air/" ;
}
3. 后處理查詢腳本支持
本質(zhì)上數(shù)據(jù)整理的目的還是為了用, 查詢或者后處理等等, 目前先搞了個(gè)簡(jiǎn)單的查詢并且屏幕打印的功能
通過(guò)config.py來(lái)控制參數(shù), 如下
# coding=utf-8
name = 'test'
mode = 'slice'
sites = ['1223A', '1224A']
datehours = ['2022070100', '2022070101', '2021080100']
variables = ['PM2.5', 'O3']
后續(xù)還得持續(xù)優(yōu)化, 爭(zhēng)取把什么數(shù)據(jù)提取到nc/db/csv, 然后對(duì)應(yīng)wrf-cmaq網(wǎng)格化, 計(jì)算MDA8等等常規(guī)處理流程都搞起來(lái)
代碼放在了這里, 僅供參考
