1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
from math import radians,sin,cos,asin,sqrt
import xarray as xr
import pandas as pd
import numpy as np
year = 2016
stime = "2016-12-11 00:00:00"
etime = "2016-12-19 23:00:00"
GeoFile = "/r008/xdxie/data/wrfchemdata/beijing201612/presource/domain/geo_em.d02.nc"
OutFile = "/r008/xdxie/data/wrfchemdata/beijing201612/observation/obs.d02.timemean_20161211-1219.csv"
inpdatadir = '/r008/xdxie/data/obs/airquality'
f = xr.open_dataset(GeoFile)
lat2d = f['XLAT_M'].sel(Time=0)
lon2d = f['XLONG_M'].sel(Time=0)
allsite = pd.read_csv(f"{inpdatadir}/obssiteinfo.csv",encoding='gbk')
def haversine_dis(lon1, lat1, lon2, lat2):
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
d_lon = lon2 - lon1
d_lat = lat2 - lat1
aa = sin(d_lat/2)**2 + cos(lat1)*cos(lat2)*sin(d_lon/2)**2
c = 2 * asin(sqrt(aa))
r = 6371
return c*r*1000
def findij(latlon):
a = abs(lat2d-latlon['lat'])+abs(lon2d-latlon['lon'])
i,j = np.unravel_index(a.argmin(),a.shape)
if (haversine_dis(lon2d[i,j], lat2d[i,j], latlon['lon'], latlon['lat']) > f.DX):
return np.nan
else:
return [i,j]
def read_station_hourly_data(StationCode, year, stime, etime):
temp = allsite.query("StationCode == @StationCode")
ProvinceCode = temp.ProvinceCode.values[0]
CityCode = temp.CityCode.values[0]
df = pd.read_csv(f"{inpdatadir}/QC+PPB/hour/{year}/{ProvinceCode}_{CityCode}_{StationCode}_{year}.csv",
index_col="datetime")
df.index = pd.to_datetime(df.index,format="%Y-%m-%d %H:%M:%S")
return df.loc[stime:etime]
latlons = allsite.loc[:,['lat','lon']]
a = latlons.apply(findij,axis=1)
a.name = "ij"
subsite = pd.concat([allsite, a], axis=1).dropna(subset=['ij'])
dfs = pd.DataFrame(columns=["station",'ix','jx','pm25', 'pm10', 'o3', 'co', 'so2', 'no2'])
for StationCode,ij in zip(subsite.StationCode,subsite.ij):
temp = read_station_hourly_data(StationCode, year, stime, etime).mean()
dfs = pd.concat([dfs,pd.DataFrame.from_records([{"station":f"{ij[0]}_{ij[1]}",
'ix':ij[0],
'jx':ij[1],
'pm25':temp.pm25,
'pm10':temp.pm10,
'o3':temp.o3,
'co':temp.co,
'so2':temp.so2,
'no2':temp.no2
}])])
dfs.groupby("station").mean().to_csv(OutFile)
|