[4766] | 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-2020 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() |
---|