1 | #!/usr/bin/env python3 |
---|
2 | # -*- coding: utf-8 -*- |
---|
3 | |
---|
4 | #------------------------------------------------------------------------------# |
---|
5 | # |
---|
6 | # Scripts for processing of WRF and CAMx files to PALM dynamic driver |
---|
7 | # |
---|
8 | # This program is free software: you can redistribute it and/or modify |
---|
9 | # it under the terms of the GNU General Public License as published by |
---|
10 | # the Free Software Foundation, either version 3 of the License, or |
---|
11 | # (at your option) any later version. |
---|
12 | # |
---|
13 | # This program is distributed in the hope that it will be useful, |
---|
14 | # but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
15 | # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
16 | # GNU General Public License for more details. |
---|
17 | # |
---|
18 | # You should have received a copy of the GNU General Public License |
---|
19 | # along with this program. If not, see <https://www.gnu.org/licenses/>. |
---|
20 | # |
---|
21 | # Copyright 2018-2021 Institute of Computer Science |
---|
22 | # of the Czech Academy of Sciences, Prague |
---|
23 | # Authors: Krystof Eben, Jaroslav Resler, Pavel Krc |
---|
24 | # |
---|
25 | #------------------------------------------------------------------------------# |
---|
26 | |
---|
27 | '''WRF (+CAMx) utility module for PALM dynamic driver generator''' |
---|
28 | |
---|
29 | import os |
---|
30 | import math |
---|
31 | import numpy as np |
---|
32 | import pyproj |
---|
33 | import scipy.ndimage as ndimage |
---|
34 | import netCDF4 |
---|
35 | |
---|
36 | import metpy.calc as mpcalc |
---|
37 | from metpy.interpolate import interpolate_1d, log_interpolate_1d |
---|
38 | from metpy.units import units |
---|
39 | |
---|
40 | from palm_dynamic_config import * |
---|
41 | |
---|
42 | # Constants directly equivalent to WRF code |
---|
43 | radius = 6370000.0 |
---|
44 | g = 9.81 #m/s2 |
---|
45 | rd = 287. #dry air gas constant (J/kg/K) |
---|
46 | rd_cp = 2./7. #from WRF v4 technote (R_d / c_p) |
---|
47 | wrf_base_temp = 300. #NOT wrfout T00 |
---|
48 | |
---|
49 | _ax = np.newaxis |
---|
50 | |
---|
51 | class WRFCoordTransform(object): |
---|
52 | 'Coordinate transformer for WRFOUT files' |
---|
53 | |
---|
54 | def __init__(self, ncf): |
---|
55 | attr = lambda a: getattr(ncf, a) |
---|
56 | |
---|
57 | # Define grids |
---|
58 | |
---|
59 | # see http://www.pkrc.net/wrf-lambert.html |
---|
60 | #latlon_wgs84 = pyproj.Proj(proj='latlong', |
---|
61 | # ellps='WGS84', datum='WGS84', |
---|
62 | # no_defs=True) #don't use - WRF datum misuse |
---|
63 | |
---|
64 | latlon_sphere = pyproj.Proj(proj='latlong', |
---|
65 | a=radius, b=radius, |
---|
66 | towgs84='0,0,0', no_defs=True) |
---|
67 | |
---|
68 | lambert_grid = pyproj.Proj(proj='lcc', |
---|
69 | lat_1=attr('TRUELAT1'), |
---|
70 | lat_2=attr('TRUELAT2'), |
---|
71 | lat_0=attr('MOAD_CEN_LAT'), |
---|
72 | lon_0=attr('STAND_LON'), |
---|
73 | a=radius, b=radius, |
---|
74 | towgs84='0,0,0', no_defs=True) |
---|
75 | |
---|
76 | # resoltion in m |
---|
77 | self.dx = dx = attr('DX') |
---|
78 | self.dy = dy = attr('DY') |
---|
79 | |
---|
80 | # number of mass grid points |
---|
81 | self.nx = nx = attr('WEST-EAST_GRID_DIMENSION') - 1 |
---|
82 | self.ny = ny = attr('SOUTH-NORTH_GRID_DIMENSION') - 1 |
---|
83 | |
---|
84 | # distance between centers of mass grid points at edges |
---|
85 | extent_x = (nx - 1) * dx |
---|
86 | extent_y = (ny - 1) * dy |
---|
87 | |
---|
88 | # grid center in lambert |
---|
89 | center_x, center_y = pyproj.transform(latlon_sphere, lambert_grid, |
---|
90 | attr('CEN_LON'), attr('CEN_LAT')) |
---|
91 | |
---|
92 | # grid origin coordinates in lambert |
---|
93 | i0_x = center_x - extent_x*.5 |
---|
94 | j0_y = center_y - extent_y*.5 |
---|
95 | |
---|
96 | # Define fast transformation methods |
---|
97 | |
---|
98 | def latlon_to_ji(lat, lon): |
---|
99 | x, y = pyproj.transform(latlon_sphere, lambert_grid, |
---|
100 | lon, lat) |
---|
101 | return (y-j0_y)/dy, (x-i0_x)/dx |
---|
102 | self.latlon_to_ji = latlon_to_ji |
---|
103 | |
---|
104 | def ji_to_latlon(j, i): |
---|
105 | lon, lat = pyproj.transform(lambert_grid, latlon_sphere, |
---|
106 | i*dx+i0_x, j*dy+j0_y) |
---|
107 | return lat, lon |
---|
108 | self.ji_to_latlon = ji_to_latlon |
---|
109 | |
---|
110 | def verify(self, ncf): |
---|
111 | lat = ncf.variables['XLAT'][0] |
---|
112 | lon = ncf.variables['XLONG'][0] |
---|
113 | j, i = np.mgrid[0:self.ny, 0:self.nx] |
---|
114 | |
---|
115 | jj, ii = self.latlon_to_ji(lat, lon) |
---|
116 | d = np.hypot(jj-j, ii-i) |
---|
117 | print('error for ll->ji: max {0} m, avg {1} m.'.format(d.max(), d.mean())) |
---|
118 | |
---|
119 | llat, llon = self.ji_to_latlon(j, i) |
---|
120 | d = np.hypot(llat - lat, llon - lon) |
---|
121 | print('error for ji->ll: max {0} deg, avg {1} deg.'.format(d.max(), d.mean())) |
---|
122 | |
---|
123 | lat = ncf.variables['XLAT_U'][0] |
---|
124 | lon = ncf.variables['XLONG_U'][0] |
---|
125 | j, i = np.mgrid[0:self.ny, 0:self.nx+1] |
---|
126 | jj, ii = self.latlon_to_ji(lat, lon) |
---|
127 | ii = ii + .5 |
---|
128 | d = np.hypot(jj-j, ii-i) |
---|
129 | print('error for U-staggered ll->ji: max {0} m, avg {1} m.'.format(d.max(), d.mean())) |
---|
130 | |
---|
131 | class CAMxCoordTransform(object): |
---|
132 | 'Coordinate transformer for CAMx files running from WRF' |
---|
133 | |
---|
134 | def __init__(self, ncf): |
---|
135 | attr = lambda a: getattr(ncf, a) |
---|
136 | |
---|
137 | # Define grids |
---|
138 | |
---|
139 | latlon_sphere = pyproj.Proj(proj='latlong', |
---|
140 | a=radius, b=radius, |
---|
141 | towgs84='0,0,0', no_defs=True) |
---|
142 | |
---|
143 | lambert_grid = pyproj.Proj(proj='lcc', |
---|
144 | lat_1=attr('P_ALP'), |
---|
145 | lat_2=attr('P_BET'), |
---|
146 | lat_0=attr('YCENT'), |
---|
147 | lon_0=attr('P_GAM'), |
---|
148 | a=radius, b=radius, |
---|
149 | towgs84='0,0,0', no_defs=True) |
---|
150 | |
---|
151 | # resoltion in m |
---|
152 | self.dx = dx = attr('XCELL') |
---|
153 | self.dy = dy = attr('YCELL') |
---|
154 | |
---|
155 | # number of mass grid points |
---|
156 | self.nx = nx = attr('NCOLS') |
---|
157 | self.ny = ny = attr('NROWS') |
---|
158 | |
---|
159 | # grid origin coordinates in lambert |
---|
160 | i0_x = attr('XORIG') |
---|
161 | j0_y = attr('YORIG') |
---|
162 | |
---|
163 | # Define fast transformation methods |
---|
164 | |
---|
165 | def latlon_to_ji(lat, lon): |
---|
166 | x, y = pyproj.transform(latlon_sphere, lambert_grid, |
---|
167 | lon, lat) |
---|
168 | return (y-j0_y)/dy, (x-i0_x)/dx |
---|
169 | self.latlon_to_ji = latlon_to_ji |
---|
170 | |
---|
171 | def ji_to_latlon(j, i): |
---|
172 | lon, lat = pyproj.transform(lambert_grid, latlon_sphere, |
---|
173 | i*dx+i0_x, j*dy+j0_y) |
---|
174 | return lat, lon |
---|
175 | self.ji_to_latlon = ji_to_latlon |
---|
176 | |
---|
177 | def verify(self, ncf): |
---|
178 | lat = ncf.variables['latitude'][:] |
---|
179 | lon = ncf.variables['longitude'][:] |
---|
180 | j, i = np.mgrid[0:self.ny, 0:self.nx] |
---|
181 | |
---|
182 | jj, ii = self.latlon_to_ji(lat, lon) |
---|
183 | d = np.hypot(jj-j, ii-i) |
---|
184 | print('error for ll->ji: max {0} m, avg {1} m.'.format(d.max(), d.mean())) |
---|
185 | |
---|
186 | llat, llon = self.ji_to_latlon(j, i) |
---|
187 | d = np.hypot(llat - lat, llon - lon) |
---|
188 | print('error for ji->ll: max {0} deg, avg {1} deg.'.format(d.max(), d.mean())) |
---|
189 | |
---|
190 | class BilinearRegridder(object): |
---|
191 | '''Bilinear regridder for multidimensional data. |
---|
192 | |
---|
193 | By standard, the last two dimensions are always Y,X in that order. |
---|
194 | ''' |
---|
195 | def __init__(self, projected_x, projected_y, preloaded=False): |
---|
196 | projected_x = np.asanyarray(projected_x) |
---|
197 | projected_y = np.asanyarray(projected_y) |
---|
198 | self.shape = projected_x.shape |
---|
199 | self.rank = len(self.shape) |
---|
200 | assert self.shape == projected_y.shape |
---|
201 | |
---|
202 | y0 = np.floor(projected_y) |
---|
203 | yd = projected_y - y0 |
---|
204 | ydd = 1. - yd |
---|
205 | self.y0 = y0.astype('i8') |
---|
206 | |
---|
207 | x0 = np.floor(projected_x) |
---|
208 | xd = projected_x - x0 |
---|
209 | xdd = 1. - xd |
---|
210 | self.x0 = x0.astype('i8') |
---|
211 | |
---|
212 | if preloaded: |
---|
213 | # Prepare slices for preloading from NetCDF files (for cases where |
---|
214 | # the range of loaded Y, X coordinates is much less than total |
---|
215 | # size. The regrid method then expects preloaded data. |
---|
216 | |
---|
217 | ybase = self.y0.min() |
---|
218 | self.ys = slice(ybase, self.y0.max()+2) |
---|
219 | self.y0 -= ybase |
---|
220 | |
---|
221 | xbase = self.x0.min() |
---|
222 | self.xs = slice(xbase, self.x0.max()+2) |
---|
223 | self.x0 -= xbase |
---|
224 | |
---|
225 | self.y1 = self.y0 + 1 |
---|
226 | self.x1 = self.x0 + 1 |
---|
227 | |
---|
228 | self.weights = np.array([ |
---|
229 | ydd * xdd, #wy0x0 |
---|
230 | ydd * xd , #wy0x1 |
---|
231 | yd * xdd, #wy1x0 |
---|
232 | yd * xd , #wy1x1 |
---|
233 | ]) |
---|
234 | |
---|
235 | def regrid(self, data): |
---|
236 | # data may contain additional dimensions (before Y,X) |
---|
237 | dshape = data.shape[:-2] |
---|
238 | drank = len(dshape) |
---|
239 | |
---|
240 | # Prepare array for selected data |
---|
241 | sel_shape = (4,) + dshape + self.shape |
---|
242 | selection = np.empty(sel_shape, dtype=data.dtype) |
---|
243 | |
---|
244 | selection[0, ...] = data[..., self.y0, self.x0] |
---|
245 | selection[1, ...] = data[..., self.y0, self.x1] |
---|
246 | selection[2, ...] = data[..., self.y1, self.x0] |
---|
247 | selection[3, ...] = data[..., self.y1, self.x1] |
---|
248 | |
---|
249 | # Slice weights to match the extra dimensions |
---|
250 | wslice = ((slice(None),) + #weights |
---|
251 | (np.newaxis,) * drank + #data minus Y,X |
---|
252 | (slice(None),) * self.rank) #regridded shape |
---|
253 | |
---|
254 | w = selection * self.weights[wslice] |
---|
255 | return w.sum(axis=0) |
---|
256 | |
---|
257 | def print_dstat(desc, delta): |
---|
258 | print('Delta stats for {0} ({1:8g} ~ {2:8g}): bias = {3:8g}, MAE = {4:8g}, RMSE = {5:8g}'.format( |
---|
259 | desc, delta.min(), delta.max(), delta.mean(), np.abs(delta).mean(), np.sqrt((delta**2).mean()))) |
---|
260 | |
---|
261 | def barom_pres(p0, gp, gp0, t0): |
---|
262 | barom = 1. / (rd * t0) |
---|
263 | return p0 * np.exp((gp0-gp)*barom) |
---|
264 | |
---|
265 | def barom_gp(gp0, p, p0, t0): |
---|
266 | baromi = rd * t0 |
---|
267 | return gp0 - np.log(p/p0) * baromi |
---|
268 | |
---|
269 | def calc_ph_hybrid(f, mu): |
---|
270 | pht = f.variables['P_TOP'][0] |
---|
271 | c3f = f.variables['C3F'][0] |
---|
272 | c4f = f.variables['C4F'][0] |
---|
273 | c3h = f.variables['C3H'][0] |
---|
274 | c4h = f.variables['C4H'][0] |
---|
275 | return (c3f[:,_ax,_ax]*mu[_ax,:,:] + (c4f[:,_ax,_ax] + pht), |
---|
276 | c3h[:,_ax,_ax]*mu[_ax,:,:] + (c4h[:,_ax,_ax] + pht)) |
---|
277 | |
---|
278 | def calc_ph_sigma(f, mu): |
---|
279 | pht = f.variables['P_TOP'][0] |
---|
280 | eta_f = f.variables['ZNW'][0] |
---|
281 | eta_h = f.variables['ZNU'][0] |
---|
282 | return (eta_f[:,_ax,_ax]*mu[_ax,:,:] + pht, |
---|
283 | eta_h[:,_ax,_ax]*mu[_ax,:,:] + pht) |
---|
284 | |
---|
285 | def wrf_t(f): |
---|
286 | p = f.variables['P'][0,:,:,:] + f.variables['PB'][0,:,:,:] |
---|
287 | return (f.variables['T'][0,:,:,:] + wrf_base_temp) * np.power(0.00001*p, rd_cp) |
---|
288 | |
---|
289 | def calc_gp(f, ph): |
---|
290 | terr = f.variables['HGT'][0,:,:] |
---|
291 | gp0 = terr * g |
---|
292 | gp = [gp0] |
---|
293 | t = wrf_t(f) |
---|
294 | for lev in range(1, ph.shape[0]): |
---|
295 | gp.append(barom_gp(gp[-1], ph[lev,:,:], ph[lev-1,:,:], t[lev-1,:,:])) |
---|
296 | return np.array(gp) |
---|
297 | |
---|
298 | def palm_wrf_vertical_interp(infile, outfile, wrffile, z_levels, z_levels_stag, |
---|
299 | z_soil_levels, origin_z, terrain, wrf_hybrid_levs, vinterp_terrain_smoothing): |
---|
300 | |
---|
301 | zdim = len(z_levels) |
---|
302 | zwdim = len(z_levels_stag) |
---|
303 | zsoildim = len(z_soil_levels) |
---|
304 | dimnames = ['z', 'zw', 'zsoil'] # dimnames of palm vertical dims |
---|
305 | dimensions = [zdim , zwdim, zsoildim] |
---|
306 | |
---|
307 | print("infile: " , infile) |
---|
308 | print("outfile: ", outfile) |
---|
309 | |
---|
310 | try: |
---|
311 | os.remove(outfile) |
---|
312 | os.remove(infile+'_vinterp.log') |
---|
313 | except: |
---|
314 | pass |
---|
315 | |
---|
316 | nc_infile = netCDF4.Dataset(infile, 'r') |
---|
317 | nc_wrf = netCDF4.Dataset(wrffile, 'r') |
---|
318 | nc_outfile = netCDF4.Dataset(outfile, "w", format="NETCDF4") |
---|
319 | nc_outfile.createDimension('Time', None) |
---|
320 | for dimname in ['west_east', 'south_north', 'soil_layers_stag']: |
---|
321 | nc_outfile.createDimension(dimname, len(nc_infile.dimensions[dimname])) |
---|
322 | for i in range (len(dimnames)): |
---|
323 | nc_outfile.createDimension(dimnames[i], dimensions[i]) |
---|
324 | |
---|
325 | # Use hybrid ETA levels in WRF and stretch them so that the WRF terrain |
---|
326 | # matches either PALM terrain or flat terrain at requested height |
---|
327 | gpf = nc_infile.variables['PH'][0,:,:,:] + nc_infile.variables['PHB'][0,:,:,:] |
---|
328 | wrfterr = gpf[0]*(1./g) |
---|
329 | |
---|
330 | if vinterp_terrain_smoothing is None: |
---|
331 | target_terrain = terrain |
---|
332 | else: |
---|
333 | print('Smoothing PALM terrain for the purpose of dynamic driver with sigma={0} grid points.'.format( |
---|
334 | vinterp_terrain_smoothing)) |
---|
335 | target_terrain = ndimage.gaussian_filter(terrain, sigma=vinterp_terrain_smoothing, order=0) |
---|
336 | print('Morphing WRF terrain ({0} ~ {1}) to PALM terrain ({2} ~ {3})'.format( |
---|
337 | wrfterr.min(), wrfterr.max(), target_terrain.min(), target_terrain.max())) |
---|
338 | print_dstat('terrain shift', wrfterr - target_terrain[:,:]) |
---|
339 | |
---|
340 | # Load original dry air column pressure |
---|
341 | mu = nc_infile.variables['MUB'][0,:,:] + nc_infile.variables['MU'][0,:,:] |
---|
342 | pht = nc_wrf.variables['P_TOP'][0] |
---|
343 | |
---|
344 | # Shift column pressure so that it matches PALM terrain |
---|
345 | t = wrf_t(nc_infile) |
---|
346 | mu2 = barom_pres(mu+pht, target_terrain*g, gpf[0,:,:], t[0,:,:])-pht |
---|
347 | |
---|
348 | # Calculate original and shifted 3D dry air pressure |
---|
349 | if wrf_hybrid_levs: |
---|
350 | phf, phh = calc_ph_hybrid(nc_wrf, mu) |
---|
351 | phf2, phh2 = calc_ph_hybrid(nc_wrf, mu2) |
---|
352 | else: |
---|
353 | phf, phh = calc_ph_sigma(nc_wrf, mu) |
---|
354 | phf2, phh2 = calc_ph_sigma(nc_wrf, mu2) |
---|
355 | |
---|
356 | # Shift 3D geopotential according to delta dry air pressure |
---|
357 | tf = np.concatenate((t, t[-1:,:,:]), axis=0) # repeat highest layer |
---|
358 | gpf2 = barom_gp(gpf, phf2, phf, tf) |
---|
359 | # For half-levs, originate from gp full levs rather than less accurate gp halving |
---|
360 | gph2 = barom_gp(gpf[:-1,:,:], phh2, phf[:-1,:,:], t) |
---|
361 | zf = gpf2 * (1./g) - origin_z |
---|
362 | zh = gph2 * (1./g) - origin_z |
---|
363 | |
---|
364 | # Report |
---|
365 | gpdelta = gpf2 - gpf |
---|
366 | print('GP deltas by level:') |
---|
367 | for k in range(gpf.shape[0]): |
---|
368 | print_dstat(k, gpdelta[k]) |
---|
369 | |
---|
370 | # Because we require levels below the lowest level from WRF, we will always |
---|
371 | # add one layer at zero level with repeated values from the lowest level. |
---|
372 | # WRF-python had some special treatment for theta in this case. |
---|
373 | height = np.zeros((zh.shape[0]+1,) + zh.shape[1:], dtype=zh.dtype) |
---|
374 | height[0,:,:] = -999. #always below terrain |
---|
375 | height[1:,:,:] = zh |
---|
376 | heightw = np.zeros((zf.shape[0]+1,) + zf.shape[1:], dtype=zf.dtype) |
---|
377 | heightw[0,:,:] = -999. #always below terrain |
---|
378 | heightw[1:,:,:] = zf |
---|
379 | |
---|
380 | # ======================== SPECIFIC HUMIDITY ============================== |
---|
381 | qv_raw = nc_infile.variables['SPECHUM'][0] |
---|
382 | qv_raw = np.r_[qv_raw[0:1], qv_raw] |
---|
383 | |
---|
384 | # Vertical interpolation to grid height levels (specified in km!) |
---|
385 | # Levels start at 50m (below that the interpolation looks very sketchy) |
---|
386 | init_atmosphere_qv = interpolate_1d(z_levels, height, qv_raw) |
---|
387 | vdata = nc_outfile.createVariable('init_atmosphere_qv', "f4", ("Time", "z","south_north","west_east")) |
---|
388 | vdata[0,:,:,:] = init_atmosphere_qv |
---|
389 | |
---|
390 | # ======================= POTENTIAL TEMPERATURE ========================== |
---|
391 | pt_raw = nc_infile.variables['T'][0] + 300. # from perturbation pt to standard |
---|
392 | pt_raw = np.r_[pt_raw[0:1], pt_raw] |
---|
393 | init_atmosphere_pt = interpolate_1d(z_levels, height, pt_raw) |
---|
394 | |
---|
395 | #plt.figure(); plt.contourf(pt[0]) ; plt.colorbar() ; plt.show() |
---|
396 | vdata = nc_outfile.createVariable('init_atmosphere_pt', "f4", ("Time", "z","south_north","west_east")) |
---|
397 | vdata[0,:,:,:] = init_atmosphere_pt |
---|
398 | |
---|
399 | # ======================= Wind ========================================== |
---|
400 | u_raw = nc_infile.variables['U'][0] |
---|
401 | u_raw = np.r_[u_raw[0:1], u_raw] |
---|
402 | init_atmosphere_u = interpolate_1d(z_levels, height, u_raw) |
---|
403 | |
---|
404 | vdata = nc_outfile.createVariable('init_atmosphere_u', "f4", ("Time", "z","south_north","west_east")) |
---|
405 | vdata[0,:,:,:] = init_atmosphere_u |
---|
406 | |
---|
407 | v_raw = nc_infile.variables['V'][0] |
---|
408 | v_raw = np.r_[v_raw[0:1], v_raw] |
---|
409 | init_atmosphere_v = interpolate_1d(z_levels, height, v_raw) |
---|
410 | |
---|
411 | vdata = nc_outfile.createVariable('init_atmosphere_v', "f4", ("Time", "z","south_north","west_east")) |
---|
412 | #vdata.coordinates = "XLONG_V XLAT_V XTIME" |
---|
413 | vdata[0,:,:,:] = init_atmosphere_v |
---|
414 | |
---|
415 | w_raw = nc_infile.variables['W'][0] |
---|
416 | w_raw = np.r_[w_raw[0:1], w_raw] |
---|
417 | init_atmosphere_w = interpolate_1d(z_levels_stag, heightw, w_raw) |
---|
418 | |
---|
419 | vdata = nc_outfile.createVariable('init_atmosphere_w', "f4", ("Time", "zw","south_north","west_east")) |
---|
420 | #vdata.coordinates = "XLONG XLAT XTIME" |
---|
421 | vdata[0,:,:,:] = init_atmosphere_w |
---|
422 | |
---|
423 | # ===================== SURFACE PRESSURE ================================== |
---|
424 | surface_forcing_surface_pressure = nc_infile.variables['PSFC'] |
---|
425 | vdata = nc_outfile.createVariable('surface_forcing_surface_pressure', "f4", ("Time", "south_north","west_east")) |
---|
426 | vdata[0,:,:] = surface_forcing_surface_pressure[0,:,:] |
---|
427 | |
---|
428 | # ======================== SOIL VARIABLES (without vertical interpolation) ============= |
---|
429 | # soil temperature |
---|
430 | init_soil_t = nc_infile.variables['TSLB'] |
---|
431 | vdata = nc_outfile.createVariable('init_soil_t', "f4", ("Time", "zsoil","south_north","west_east")) |
---|
432 | vdata[0,:,:,:] = init_soil_t[0,:,:,:] |
---|
433 | |
---|
434 | # soil moisture |
---|
435 | init_soil_m = nc_infile.variables['SMOIS'] |
---|
436 | vdata = nc_outfile.createVariable('init_soil_m', "f4", ("Time","zsoil","south_north","west_east")) |
---|
437 | vdata[0,:,:,:] = init_soil_m[0,:,:,:] |
---|
438 | |
---|
439 | # zsoil |
---|
440 | zsoil = nc_wrf.variables['ZS'] #ZS:description = "DEPTHS OF CENTERS OF SOIL LAYERS" ; |
---|
441 | vdata = nc_outfile.createVariable('zsoil', "f4", ("zsoil")) |
---|
442 | vdata[:] = zsoil[0,:] |
---|
443 | |
---|
444 | # coordinates z, zw |
---|
445 | vdata = nc_outfile.createVariable('z', "f4", ("z")) |
---|
446 | vdata[:] = list(z_levels) |
---|
447 | |
---|
448 | vdata = nc_outfile.createVariable('zw', "f4", ("zw")) |
---|
449 | vdata[:] = list (z_levels_stag) |
---|
450 | |
---|
451 | # zsoil is taken from wrf - not need to define it |
---|
452 | |
---|
453 | nc_infile.close() |
---|
454 | nc_wrf.close() |
---|
455 | nc_outfile.close() |
---|
456 | |
---|
457 | def palm_wrf_gw(f, lon, lat, levels, tidx=0): |
---|
458 | '''Calculate geostrophic wind from WRF using metpy''' |
---|
459 | |
---|
460 | hgts, ug, vg = calcgw_wrf(f, lat, lon, levels, tidx) |
---|
461 | |
---|
462 | # extrapolate at the bottom |
---|
463 | hgts = np.r_[np.array([0.]), hgts] |
---|
464 | ug = np.r_[ug[0], ug] |
---|
465 | vg = np.r_[vg[0], vg] |
---|
466 | |
---|
467 | return minterp(levels, hgts, ug, vg) |
---|
468 | |
---|
469 | |
---|
470 | def minterp(interp_heights, data_heights, u, v): |
---|
471 | '''Interpolate wind using power law for agl levels''' |
---|
472 | |
---|
473 | pdata = data_heights ** gw_alpha |
---|
474 | pinterp = interp_heights ** gw_alpha |
---|
475 | hindex = np.searchsorted(data_heights, interp_heights, side='right') |
---|
476 | lindex = hindex - 1 |
---|
477 | assert lindex[0] >= 0 |
---|
478 | assert hindex[-1] < len(data_heights) |
---|
479 | lbound = pdata[lindex] |
---|
480 | hcoef = (pinterp - lbound) / (pdata[hindex] - lbound) |
---|
481 | #print(data_heights) |
---|
482 | #print(lindex) |
---|
483 | #print(hcoef) |
---|
484 | lcoef = 1. - hcoef |
---|
485 | iu = u[lindex] * lcoef + u[hindex] * hcoef |
---|
486 | iv = v[lindex] * lcoef + v[hindex] * hcoef |
---|
487 | return iu, iv |
---|
488 | |
---|
489 | def get_wrf_dims(f, lat, lon, xlat, xlong): |
---|
490 | '''A crude method, yet satisfactory for approximate WRF surroundings''' |
---|
491 | |
---|
492 | sqdist = (xlat - lat)**2 + (xlong - lon)**2 |
---|
493 | coords = np.unravel_index(sqdist.argmin(), sqdist.shape) |
---|
494 | |
---|
495 | xmargin = int(math.ceil(gw_wrf_margin_km * 1000 / f.DX)) #py2 ceil produces float |
---|
496 | ymargin = int(math.ceil(gw_wrf_margin_km * 1000 / f.DY)) |
---|
497 | y0, y1 = coords[0] - ymargin, coords[0] + ymargin |
---|
498 | x0, x1 = coords[1] - xmargin, coords[1] + xmargin |
---|
499 | assert 0 <= y0 < y1 < sqdist.shape[0], "Point {0} + surroundings not inside domain".format(coords[0]) |
---|
500 | assert 0 <= x0 < x1 < sqdist.shape[1], "Point {0} + surroundings not inside domain".format(coords[1]) |
---|
501 | |
---|
502 | return coords, (slice(y0, y1+1), slice(x0, x1+1)), (ymargin, xmargin) |
---|
503 | |
---|
504 | def calcgw_wrf(f, lat, lon, levels, tidx=0): |
---|
505 | # MFDataset removes the time dimension from XLAT, XLONG |
---|
506 | xlat = f.variables['XLAT'] |
---|
507 | xlslice = (0,) * (len(xlat.shape)-2) + (slice(None), slice(None)) |
---|
508 | xlat = xlat[xlslice] |
---|
509 | xlong = f.variables['XLONG'][xlslice] |
---|
510 | |
---|
511 | (iy, ix), area, (iby, ibx) = get_wrf_dims(f, lat, lon, xlat, xlong) |
---|
512 | areat = (tidx,) + area |
---|
513 | areatz = (tidx, slice(None)) + area |
---|
514 | #print('wrf coords', lat, lon, xlat[iy,ix], xlong[iy,ix]) |
---|
515 | #print(xlat[area][iby,ibx], xlong[area][iby,ibx], areat) |
---|
516 | |
---|
517 | # load area |
---|
518 | hgt = (f.variables['PH'][areatz] + f.variables['PHB'][areatz]) / 9.81 |
---|
519 | hgtu = (hgt[:-1] + hgt[1:]) * .5 |
---|
520 | pres = f.variables['P'][areatz] + f.variables['PB'][areatz] |
---|
521 | terrain = f.variables['HGT'][areat] |
---|
522 | |
---|
523 | # find suitable pressure levels |
---|
524 | yminpres, xminpres = np.unravel_index(pres[0].argmin(), pres[0].shape) |
---|
525 | pres1 = pres[0, yminpres, xminpres] - 1. |
---|
526 | |
---|
527 | aglpt = hgtu[:,iby,ibx] - terrain[iby,ibx] |
---|
528 | pres0 = pres[np.searchsorted(aglpt, levels[-1]), iby, ibx] |
---|
529 | plevels = np.arange(pres1, min(pres0, pres1)-1, -1000.) |
---|
530 | |
---|
531 | # interpolate wrf into pressure levels |
---|
532 | phgt = log_interpolate_1d(plevels, pres, hgtu, axis=0) |
---|
533 | |
---|
534 | # Set up some constants based on our projection, including the Coriolis parameter and |
---|
535 | # grid spacing, converting lon/lat spacing to Cartesian |
---|
536 | coriol = mpcalc.coriolis_parameter(np.deg2rad(xlat[area])).to('1/s') |
---|
537 | |
---|
538 | # lat_lon_grid_deltas doesn't work under py2, but for WRF grid it is still |
---|
539 | # not very accurate, better use direct values. |
---|
540 | #dx, dy = mpcalc.lat_lon_grid_deltas(xlong[area], xlat[area]) |
---|
541 | dx = f.DX * units.m |
---|
542 | dy = f.DY * units.m |
---|
543 | |
---|
544 | # Smooth height data. Sigma=1.5 for gfs 0.5deg |
---|
545 | res_km = f.DX / 1000. |
---|
546 | |
---|
547 | ug = np.zeros(plevels.shape, 'f8') |
---|
548 | vg = np.zeros(plevels.shape, 'f8') |
---|
549 | for i in range(len(plevels)): |
---|
550 | sh = ndimage.gaussian_filter(phgt[i,:,:], sigma=1.5*50/res_km, order=0) |
---|
551 | ugl, vgl = mpcalc.geostrophic_wind(sh * units.m, coriol, dx, dy) |
---|
552 | ug[i] = ugl[iby, ibx].magnitude |
---|
553 | vg[i] = vgl[iby, ibx].magnitude |
---|
554 | |
---|
555 | return phgt[:,iby,ibx], ug, vg |
---|
556 | |
---|
557 | # The following two functions calculate GW from GFS files, although this |
---|
558 | # function is currently not implemented in PALM dynamic driver generation |
---|
559 | # script |
---|
560 | |
---|
561 | def calcgw_gfs(v, lat, lon): |
---|
562 | height, lats, lons = v.data(lat1=lat-gw_gfs_margin_deg ,lat2=lat+gw_gfs_margin_deg, |
---|
563 | lon1=lon-gw_gfs_margin_deg, lon2=lon+gw_gfs_margin_deg) |
---|
564 | i = np.searchsorted(lats[:,0], lat) |
---|
565 | if abs(lats[i+1,0] - lat) < abs(lats[i,0] - lat): |
---|
566 | i = i+1 |
---|
567 | j = np.searchsorted(lons[0,:], lon) |
---|
568 | if abs(lons[0,i+1] - lon) < abs(lons[0,i] - lon): |
---|
569 | j = j+1 |
---|
570 | #print('level', v.level, 'height', height[i,j], lats[i,j], lons[i,j]) |
---|
571 | |
---|
572 | # Set up some constants based on our projection, including the Coriolis parameter and |
---|
573 | # grid spacing, converting lon/lat spacing to Cartesian |
---|
574 | f = mpcalc.coriolis_parameter(np.deg2rad(lats)).to('1/s') |
---|
575 | dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats) |
---|
576 | res_km = (dx[i,j]+dy[i,j]).magnitude / 2000. |
---|
577 | |
---|
578 | # Smooth height data. Sigma=1.5 for gfs 0.5deg |
---|
579 | height = ndimage.gaussian_filter(height, sigma=1.5*50/res_km, order=0) |
---|
580 | |
---|
581 | # In MetPy 0.5, geostrophic_wind() assumes the order of the dimensions is (X, Y), |
---|
582 | # so we need to transpose from the input data, which are ordered lat (y), lon (x). |
---|
583 | # Once we get the components,transpose again so they match our original data. |
---|
584 | geo_wind_u, geo_wind_v = mpcalc.geostrophic_wind(height * units.m, f, dx, dy) |
---|
585 | |
---|
586 | return height[i,j], geo_wind_u[i,j], geo_wind_v[i,j] |
---|
587 | |
---|
588 | def combinegw_gfs(grbs, levels, lat, lon): |
---|
589 | heights = [] |
---|
590 | us = [] |
---|
591 | vs = [] |
---|
592 | for grb in grbs: |
---|
593 | h, u, v = calcgw_gfs(grb, lat, lon) |
---|
594 | heights.append(h) |
---|
595 | us.append(u.magnitude) |
---|
596 | vs.append(v.magnitude) |
---|
597 | heights = np.array(heights) |
---|
598 | us = np.array(us) |
---|
599 | vs = np.array(vs) |
---|
600 | |
---|
601 | ug, vg = minterp(np.asanyarray(levels), heights[::-1], us[::-1], vs[::-1]) |
---|
602 | return ug, vg |
---|
603 | |
---|
604 | if __name__ == '__main__': |
---|
605 | import argparse |
---|
606 | |
---|
607 | parser = argparse.ArgumentParser(description=__doc__) |
---|
608 | parser.add_argument('-w', '--wrfout', help='verify wrfout file') |
---|
609 | parser.add_argument('-c', '--camx', help='verify camx file') |
---|
610 | args = parser.parse_args() |
---|
611 | |
---|
612 | if args.wrfout: |
---|
613 | f = netCDF4.Dataset(args.wrfout) |
---|
614 | |
---|
615 | print('Verifying coord transform:') |
---|
616 | t = WRFCoordTransform(f) |
---|
617 | t.verify(f) |
---|
618 | |
---|
619 | print('\nVerifying vertical levels:') |
---|
620 | mu = f.variables['MUB'][0,:,:] + f.variables['MU'][0,:,:] |
---|
621 | gp = f.variables['PH'][0,:,:,:] + f.variables['PHB'][0,:,:,:] |
---|
622 | |
---|
623 | print('\nUsing sigma:') |
---|
624 | phf, phh = calc_ph_sigma(f, mu) |
---|
625 | gp_calc = calc_gp(f, phf) |
---|
626 | delta = gp_calc - gp |
---|
627 | for lev in range(delta.shape[0]): |
---|
628 | print_dstat(lev, delta[lev]) |
---|
629 | |
---|
630 | print('\nUsing hybrid:') |
---|
631 | phf, phh = calc_ph_hybrid(f, mu) |
---|
632 | gp_calc = calc_gp(f, phf) |
---|
633 | delta = gp_calc - gp |
---|
634 | for lev in range(delta.shape[0]): |
---|
635 | print_dstat(lev, delta[lev]) |
---|
636 | |
---|
637 | f.close() |
---|
638 | |
---|
639 | if args.camx: |
---|
640 | f = netCDF4.Dataset(args.camx) |
---|
641 | t = CAMxCoordTransform(f) |
---|
642 | t.verify(f) |
---|
643 | f.close() |
---|