-
Notifications
You must be signed in to change notification settings - Fork 11
/
interpolated_model_timeseries.py
97 lines (83 loc) · 3.41 KB
/
interpolated_model_timeseries.py
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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#!/usr/bin/env python3
import xarray as xr
import pandas as pd
import numpy as np
from scipy.interpolate import RegularGridInterpolator
locations = [(25, 100), (26, 101), (27, 102), (28, 103), (29, 101), (30, 101)]
wrfouts = ['/nfs/a68/eebjs/wrfoutput/p2run/202/ctl2016/feb/wrfout_feb2016_ctl_rgdd25.nc',
'/nfs/a68/eebjs/wrfoutput/p2run/202/ctl2016/mar/wrfout_mar2016_ctl_rgdd25.nc',
'/nfs/a68/eebjs/wrfoutput/p2run/202/ctl2016/apr/wrfout_apr2016_ctl_rgdd25.nc',
'/nfs/a68/eebjs/wrfoutput/p2run/202/ctl2016/may/wrfout_may2016_ctl_rgdd25.nc']
#wrfouts = ['/nfs/a336/earlacoa/paper_aia_china/test_file/wrfout_combined-domains_global_0.25deg_2015-01_PM2_5_DRY.nc']
pols = ['PM2_5_DRY', 'PM10']
level = 0
with xr.open_dataset(wrfouts[0]) as ds:
if 'time' in ds[pol].coords:
tdim = 'time'
elif 'Time' in ds[pol].coords:
tdim = 'Time'
def get_catda(wrfouts, pol, level, tdim):
"""
Description:
Returns a concatenated data array of the files contained within.
Args:
wrfouts (list): List of wrfout files.
pol (str): Pollutant.
level (int): Model level as an index location.
tdim (str): Name of the time dimension.
Returns:
catda (xarray DataArray): Concatenated data array.
"""
da_list = []
for filepath in wrfouts:
with xr.open_dataset(filepath) as ds:
da = ds[pol]
tdim_len = len(da[tdim])
if tdim_len > 1:
da = da.chunk({tdim:tdim_len//10})
if 'bottom_top' in da.dims:
if level is None:
print('wrfouts contain z dimension, please supply a level')
da = da.loc[{'bottom_top':level}]
da_list.append(da)
catda = xr.concat(da_list, dim=tdim)
return catda
def interpolate_model_timeseries(locations, tdim, times='all', level=None,
pols):
"""
Description:
Interpolates timeseries at given lat/lon locations.
The wrfout files must be regridded to a rectilinear grid.
Args:
locations (list of tuples): Observation lat lon locations.
tdim (str): Name of the time dimension.
times (str, optional): Times of interpolation, default = 'all'.
Returns:
df (pandas DataFrame): Interpolated timeseries.
"""
pol_dfs = {}
for pol in pols:
print(pol+':')
catda = get_catda(wrfouts, pol, level=level, tdim=tdim)
time_coord = np.arange(0, len(catda[tdim]))
print('creating interpolator...', end=' ')
f = RegularGridInterpolator((time_coord,
catda.coords['lat'].values,
catda.coords['lon'].values),
catda.values
)
print('done')
df = pd.DataFrame(index=catda[tdim].values)
print('sampling locations...', end=' ')
for lat, lon in locations:
indexer = np.column_stack([time_coord,
[lat]*len(time_coord),
[lon]*len(time_coord)])
series = f(indexer)
df[(lat,lon)] = series
print('done')
pol_dfs[pol] = df
df = pd.concat(pol_dfs, axis=1)
return df
if __name__ == '__main__':
df = interpolate_model_timeseries(locations, tdim, times='all', level=level)