3 |
michaesp |
1 |
# -*- coding:utf-8 -*-
|
|
|
2 |
import cartopy
|
|
|
3 |
import copy
|
|
|
4 |
import math
|
|
|
5 |
import numpy as np
|
|
|
6 |
from pyproj import Geod
|
|
|
7 |
from datetime import timedelta
|
|
|
8 |
|
|
|
9 |
from dypy.intergrid import Intergrid
|
|
|
10 |
from dypy.constants import rdv, L, R, cp, t_zero
|
|
|
11 |
|
|
|
12 |
__all__ = ['CrossSection', 'interpolate', 'rotate_points', 'interval',
|
|
|
13 |
'dewpoint', 'esat', 'moist_lapse']
|
|
|
14 |
|
|
|
15 |
|
|
|
16 |
class CrossSection(object):
|
|
|
17 |
"""
|
|
|
18 |
Create a cross-section
|
|
|
19 |
return the distance, the pressure levels,
|
|
|
20 |
and the variables.
|
|
|
21 |
"""
|
|
|
22 |
|
|
|
23 |
def __init__(self, variables, coo, pressure,
|
|
|
24 |
int2p=False, pollon=-170, pollat=43, nbre=1000, order=1):
|
|
|
25 |
"""
|
|
|
26 |
variables: dictionary of variables {'name': np.array},
|
|
|
27 |
need to contain rlat, rlon
|
|
|
28 |
coo: list of coordinates [(startlon, startlat), (endlon, endlat)]
|
|
|
29 |
pressure: pressure coordinate 1d np.array
|
|
|
30 |
int2p: True if variables need to be interpolated,
|
|
|
31 |
require P in the variables
|
|
|
32 |
pollon: pole_longitude of the rotated grid
|
|
|
33 |
pollat: pole_latitude of the rotated grid
|
|
|
34 |
nre: nbre of points along the cross-section
|
|
|
35 |
order: order of interpolation see for details
|
|
|
36 |
"""
|
|
|
37 |
variables = copy.deepcopy(variables)
|
|
|
38 |
# test the input
|
|
|
39 |
required = ['rlon', 'rlat']
|
|
|
40 |
|
|
|
41 |
if int2p:
|
|
|
42 |
required.append('p')
|
|
|
43 |
for var in required:
|
|
|
44 |
if var not in variables.keys():
|
|
|
45 |
raise Exception('{} not in variables dictionary'.format(var))
|
|
|
46 |
|
|
|
47 |
rlon = variables.pop('rlon')
|
|
|
48 |
rlat = variables.pop('rlat')
|
|
|
49 |
self.coo = coo
|
|
|
50 |
self.nbre = nbre
|
|
|
51 |
self.pressure = pressure
|
|
|
52 |
self.nz = pressure.size
|
|
|
53 |
if int2p:
|
|
|
54 |
p = variables['p']
|
|
|
55 |
|
|
|
56 |
# find the coordinate of the cross-section
|
|
|
57 |
self.lon, self.lat, crlon, crlat, self.distance = find_cross_points(
|
|
|
58 |
coo, nbre, pollon, pollat)
|
|
|
59 |
|
|
|
60 |
self.distances = np.linspace(0, self.distance, nbre)
|
|
|
61 |
|
|
|
62 |
# get the information for intergrid
|
|
|
63 |
self._lo = np.array([rlat.min(), rlon.min()])
|
|
|
64 |
self._hi = np.array([rlat.max(), rlon.max()])
|
|
|
65 |
self._query_points = [[lat, lon] for lat, lon in zip(crlat, crlon)]
|
|
|
66 |
|
|
|
67 |
if int2p:
|
|
|
68 |
for var, data in variables.items():
|
|
|
69 |
if data.ndim > 2:
|
|
|
70 |
dataint = interpolate(data, p, pressure)
|
|
|
71 |
variables[var] = dataint
|
|
|
72 |
|
|
|
73 |
variables = self._int2cross(variables)
|
|
|
74 |
self.__dict__.update(variables)
|
|
|
75 |
|
|
|
76 |
def _int2cross(self, variables):
|
|
|
77 |
for var, data in variables.items():
|
|
|
78 |
if data.squeeze().ndim > 2:
|
|
|
79 |
datacross = np.zeros((self.nz, self.nbre))
|
|
|
80 |
for i, datah in enumerate(data.squeeze()):
|
|
|
81 |
intfunc = Intergrid(datah, lo=self._lo, hi=self._hi,
|
|
|
82 |
verbose=False, order=1)
|
|
|
83 |
datacross[i, :] = intfunc.at(self._query_points)
|
|
|
84 |
else:
|
|
|
85 |
datacross = np.zeros((self.nbre, ))
|
|
|
86 |
intfunc = Intergrid(data.squeeze(), lo=self._lo, hi=self._hi,
|
|
|
87 |
verbose=False, order=1)
|
|
|
88 |
datacross[:] = intfunc.at(self._query_points)
|
|
|
89 |
variables[var] = datacross
|
|
|
90 |
return variables
|
|
|
91 |
|
|
|
92 |
|
|
|
93 |
def find_cross_points(coo, nbre, pole_longitude, pole_latitude):
|
|
|
94 |
""" give nbre points along a great circle between coo[0] and coo[1]
|
|
|
95 |
iand rotated the points
|
|
|
96 |
"""
|
|
|
97 |
|
|
|
98 |
g = Geod(ellps='WGS84')
|
|
|
99 |
cross_points = g.npts(coo[0][0], coo[0][1], coo[1][0], coo[1][1], nbre)
|
|
|
100 |
lat = np.array([point[1] for point in cross_points])
|
|
|
101 |
lon = np.array([point[0] for point in cross_points])
|
|
|
102 |
rlon, rlat = rotate_points(
|
|
|
103 |
pole_longitude, pole_latitude, lon, lat)
|
|
|
104 |
distance = great_circle_distance(coo[0][0], coo[0][1],
|
|
|
105 |
coo[1][0], coo[1][1])
|
|
|
106 |
return lon, lat, rlon, rlat, distance
|
|
|
107 |
|
|
|
108 |
|
|
|
109 |
def great_circle_distance(lon1, lat1, lon2, lat2):
|
|
|
110 |
"""
|
|
|
111 |
return the distance (km) between points following a great circle
|
|
|
112 |
|
|
|
113 |
based on : https://gist.github.com/gabesmed/1826175
|
|
|
114 |
|
|
|
115 |
>>> great_circle_distance(0, 55, 8, 45.5)
|
|
|
116 |
1199.13065955064
|
|
|
117 |
"""
|
|
|
118 |
|
|
|
119 |
EARTH_CIRCUMFERENCE = 6378137 # earth circumference in meters
|
|
|
120 |
|
|
|
121 |
dLat = math.radians(lat2 - lat1)
|
|
|
122 |
dLon = math.radians(lon2 - lon1)
|
|
|
123 |
a = (math.sin(dLat / 2) * math.sin(dLat / 2) +
|
|
|
124 |
math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) *
|
|
|
125 |
math.sin(dLon / 2) * math.sin(dLon / 2))
|
|
|
126 |
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
|
|
|
127 |
distance = EARTH_CIRCUMFERENCE * c
|
|
|
128 |
|
|
|
129 |
return distance/1000
|
|
|
130 |
|
|
|
131 |
|
|
|
132 |
def interpolate(data, grid, interplevels):
|
|
|
133 |
"""
|
|
|
134 |
intrerpolate data to grid at interplevels
|
|
|
135 |
data and grid are numpy array in the form (z,lat,lon)
|
|
|
136 |
"""
|
|
|
137 |
data = data.squeeze()
|
|
|
138 |
grid = grid.squeeze()
|
|
|
139 |
shape = list(data.shape)
|
|
|
140 |
if (data.ndim > 3) | (grid.ndim > 3):
|
|
|
141 |
message = "data and grid need to be 3d array"
|
|
|
142 |
raise IndexError(message)
|
|
|
143 |
|
|
|
144 |
try:
|
|
|
145 |
nintlev = len(interplevels)
|
|
|
146 |
except:
|
|
|
147 |
interplevels = [interplevels]
|
|
|
148 |
nintlev = len(interplevels)
|
|
|
149 |
shape[-3] = nintlev
|
|
|
150 |
|
|
|
151 |
outdata = np.ones(shape)*np.nan
|
|
|
152 |
if nintlev > 20:
|
|
|
153 |
for idx, _ in np.ndenumerate(data[0]):
|
|
|
154 |
column = grid[:, idx[0], idx[1]]
|
|
|
155 |
column_GRID = data[:, idx[0], idx[1]]
|
|
|
156 |
|
|
|
157 |
value = np.interp(
|
|
|
158 |
interplevels,
|
|
|
159 |
column,
|
|
|
160 |
column_GRID,
|
|
|
161 |
left=np.nan,
|
|
|
162 |
right=np.nan)
|
|
|
163 |
outdata[:, idx[0], idx[1]] = value[:]
|
|
|
164 |
else:
|
|
|
165 |
for j, intlevel in enumerate(interplevels):
|
|
|
166 |
for lev in range(grid.shape[0]):
|
|
|
167 |
cond1 = grid[lev, :, :] > intlevel
|
|
|
168 |
cond2 = grid[lev-1, :, :] < intlevel
|
|
|
169 |
right = np.where(cond1 & cond2)
|
|
|
170 |
if right[0].size > 0:
|
|
|
171 |
sabove = grid[lev, right[0], right[1]]
|
|
|
172 |
sbelow = grid[lev-1, right[0], right[1]]
|
|
|
173 |
dabove = data[lev, right[0], right[1]]
|
|
|
174 |
dbelow = data[lev-1, right[0], right[1]]
|
|
|
175 |
result = (intlevel-sbelow)/(sabove-sbelow) * \
|
|
|
176 |
(dabove-dbelow)+dbelow
|
|
|
177 |
outdata[j, right[0], right[1]] = result
|
|
|
178 |
return outdata
|
|
|
179 |
|
|
|
180 |
|
|
|
181 |
def rotate_points(pole_longitude, pole_latitude, lon, lat, direction='n2r'):
|
|
|
182 |
"""
|
|
|
183 |
rotate points from non-rotated system to rotated system
|
|
|
184 |
n2r : non-rotated to rotated (default)
|
|
|
185 |
r2n : rotated to non-rotated
|
|
|
186 |
return rlon,rlat,
|
|
|
187 |
"""
|
|
|
188 |
lon = np.array(lon)
|
|
|
189 |
lat = np.array(lat)
|
|
|
190 |
rotatedgrid = cartopy.crs.RotatedPole(
|
|
|
191 |
pole_longitude=pole_longitude,
|
|
|
192 |
pole_latitude=pole_latitude
|
|
|
193 |
)
|
|
|
194 |
standard_grid = cartopy.crs.Geodetic()
|
|
|
195 |
|
|
|
196 |
if direction == 'n2r':
|
|
|
197 |
rotated_points = rotatedgrid.transform_points(standard_grid, lon, lat)
|
|
|
198 |
elif direction == 'r2n':
|
|
|
199 |
rotated_points = standard_grid.transform_points(rotatedgrid, lon, lat)
|
|
|
200 |
|
|
|
201 |
rlon, rlat, _ = rotated_points.T
|
|
|
202 |
return rlon, rlat
|
|
|
203 |
|
|
|
204 |
|
|
|
205 |
def interval(starting_date, ending_date, delta=timedelta(hours=1)):
|
|
|
206 |
curr = starting_date
|
|
|
207 |
while curr < ending_date:
|
|
|
208 |
yield curr
|
|
|
209 |
curr += delta
|
|
|
210 |
|
|
|
211 |
|
|
|
212 |
def dewpoint(p, qv):
|
|
|
213 |
"""
|
|
|
214 |
Calculate the dew point temperature based on p and qv following (eq.8):
|
|
|
215 |
Lawrence, M.G., 2005. The Relationship between Relative Humidity
|
|
|
216 |
and the Dewpoint Temperature in Moist Air:
|
|
|
217 |
A Simple Conversion and Applications.
|
|
|
218 |
Bulletin of the American Meteorological Society 86,
|
|
|
219 |
225–233. doi:10.1175/BAMS-86-2-225
|
|
|
220 |
"""
|
|
|
221 |
B1 = 243.04 # °C
|
|
|
222 |
C1 = 610.94 # Pa
|
|
|
223 |
A1 = 17.625
|
|
|
224 |
|
|
|
225 |
e = p*qv/(rdv+qv)
|
|
|
226 |
td = (B1*np.log(e/C1))/(A1-np.log(e/C1))
|
|
|
227 |
|
|
|
228 |
return td
|
|
|
229 |
|
|
|
230 |
|
|
|
231 |
def esat(t):
|
|
|
232 |
"""
|
|
|
233 |
Calculate the saturation vapor pressure for t in °C
|
|
|
234 |
|
|
|
235 |
Following eq. 6 of
|
|
|
236 |
Lawrence, M.G., 2005. The Relationship between Relative Humidity
|
|
|
237 |
and the Dewpoint Temperature in Moist Air:
|
|
|
238 |
A Simple Conversion and Applications.
|
|
|
239 |
Bulletin of the American Meteorological Society 86,
|
|
|
240 |
225–233. doi:10.1175/BAMS-86-2-225
|
|
|
241 |
"""
|
|
|
242 |
C1 = 610.94 # Pa
|
|
|
243 |
A1 = 17.625 #
|
|
|
244 |
B1 = 243.03 # °C
|
|
|
245 |
|
|
|
246 |
return C1*np.exp((A1 * t) / (B1 + t))
|
|
|
247 |
|
|
|
248 |
|
|
|
249 |
def moist_lapse(t, p):
|
|
|
250 |
""" Calculates moist adiabatic lapse rate for T (Celsius) and p (Pa)
|
|
|
251 |
Note: We calculate dT/dp, not dT/dz
|
|
|
252 |
See formula 3.16 in Rogers&Yau for dT/dz, but this must be combined with
|
|
|
253 |
the dry adiabatic lapse rate (gamma = g/cp) and the
|
|
|
254 |
inverse of the hydrostatic equation (dz/dp = -RT/pg)"""
|
|
|
255 |
|
|
|
256 |
a = 2. / 7.
|
|
|
257 |
b = rdv * L * L / (R * cp)
|
|
|
258 |
c = a * L / R
|
|
|
259 |
e = esat(t)
|
|
|
260 |
wsat = rdv * e / (p - e) # Rogers&Yau 2.18
|
|
|
261 |
numer = a * (t + t_zero) + c*wsat
|
|
|
262 |
denom = p * (1 + b * wsat / ((t + t_zero)**2))
|
|
|
263 |
|
|
|
264 |
return numer / denom
|