观测数据处理-国控环境监测站点

2024-04-10 03:30:43 阅读数 727

工作需要用到我国环境监测站点的污染物浓度数据(感谢大佬的分享), 数据很全, 不过csv格式在分析的时候尤其大量数据分析的时候并不友好, 所以一般要二次处理一下

上学的时候搞过一次, 当时为了查询数据方便, 塞到了sqlite3的单文件数据库里, 不过制作起来很慢(可能是我没用并行支持好的数据库), 近期更新了一下数据且又重新调整了一下数据处理逻辑, 这里记录一下

主要的思路是以netcdf保存(此前做出来db之后还要再转nc, 麻烦的很), 然后数据规模的话, 物种一共15个, 时间是2014年到现在逐小时, 站点大概在1000多

采用每月一个nc文件存储, 每个物种单独变量, 每个变量的维度是(hour, sites)或者是(hour, cities), 因为数据除了站点的还有城市的值

然后希望处理流程尽量通用一些, 以后一些其他的观测数据也基于类似的逻辑和形式去做

那么处理流程大概是三步走

1. 制作站点列表

现有的站点信息文件有点乱, 需要制作一份完整的站点列表信息, 并且最好涵盖地理位置信息和行政区信息, 以便处理的时候筛选用, 如下 (risk后续讨论):

code,name,city,prov,nation,lat,lon,risk
1013A,市监测中心(停运150702),天津市,天津市,中国,39.097,117.151,1
1014A,南口路(停运161122),天津市,天津市,中国,39.173,117.193,0

在制作过程中大概会遇到一些问题:

数据中的站点编号对应的名字会变, 这里以新名字为准
数据中的站点编号对应的经纬度会变, 这里以新的经纬为准, 但是不知道是什么原因, 是站点搬迁了还是老的定位不准咋的, 所以基于经纬度变化幅度(0~0.1度), 给risk(风险等级)设定2-11
数据中站点的经纬度信息缺失, 根据名称自己在经纬度地图上找的点, risk等级设定12
数据中站点编号对应的我国行政区信息不统一, 这个比较麻烦, 有时候是叫法问题, 有时候是到不同级别的问题, 有时候甚至是行政区变迁的问题, 为了保持统一, 我直接根据最新的shapefile地图配合经纬度信息标注站点的行政区信息, 具体到地级市/自治州/直辖市

对于risk问题, 稳定的站点是0, 仅改名字也是0, 经纬度信息有/无变化的设为1, 2-12如上, 目的主要是在处理的时候可以做个筛选, risk越高, 观测结果和地理信息对应偏差的可能性越高, 可以考虑在对比模式验证的时候舍去试试

2. 制作nc数据文件

由于要逐月去生成, 且为了同步初始化数据和后续更新数据的逻辑, 采用的方式是先init好这个月的文件, 然后就都是根据csv文件去update了

注意的点要是:

每个月的文件要以这个月的sites为准, 不要用上一步制作的全部sites, 不过站点信息可以从上一步的结果里获取, 这一步用wsl跑个bash脚本来实现
sites和cities的结果放一个文件里
数据较多, 所以基于MS-MPI和mpi4py做了个并行支持, 单机4核一起跑也能快很多! 不用py自带的多进程是因为mpi更熟悉...而且显然更好使

最终得到的结果如下, 包含了站点信息和观测数值

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. 后处理查询脚本支持

本质上数据整理的目的还是为了用, 查询或者后处理等等, 目前先搞了个简单的查询并且屏幕打印的功能

通过config.py来控制参数, 如下

# coding=utf-8

name = 'test'
mode = 'slice'
sites = ['1223A', '1224A']
datehours = ['2022070100', '2022070101', '2021080100']
variables = ['PM2.5', 'O3']

后续还得持续优化, 争取把什么数据提取到nc/db/csv, 然后对应wrf-cmaq网格化, 计算MDA8等等常规处理流程都搞起来

代码放在了这里, 仅供参考