source: palm/trunk/UTIL/WRF_interface/dynamic/palm_wrf_utils.py @ 4866

Last change on this file since 4866 was 4843, checked in by raasch, 4 years ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

File size: 23.2 KB
Line 
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
29import os
30import math
31import numpy as np
32import pyproj
33import scipy.ndimage as ndimage
34import netCDF4
35
36import metpy.calc as mpcalc
37from metpy.interpolate import interpolate_1d, log_interpolate_1d
38from metpy.units import units
39
40from palm_dynamic_config import *
41
42# Constants directly equivalent to WRF code
43radius = 6370000.0
44g = 9.81 #m/s2
45rd = 287. #dry air gas constant (J/kg/K)
46rd_cp = 2./7. #from WRF v4 technote (R_d / c_p)
47wrf_base_temp = 300. #NOT wrfout T00
48
49_ax = np.newaxis
50
51class 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
131class 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
190class 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
257def 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
261def barom_pres(p0, gp, gp0, t0):
262    barom = 1. / (rd * t0)
263    return p0 * np.exp((gp0-gp)*barom)
264
265def barom_gp(gp0, p, p0, t0):
266    baromi = rd * t0
267    return gp0 - np.log(p/p0) * baromi
268
269def 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
278def 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
285def 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
289def 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
298def 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
457def 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
470def 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
489def 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
504def 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
561def 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
588def 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
604if __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()
Note: See TracBrowser for help on using the repository browser.