Hwt2023

README.md

@@ -1,9 +1,63 @@
 # webplot plotting package
 
-This is a simplified version of the webplot plotting library used to create graphics for the NCAR ensemble system.
+### create graphics from MPAS ensemble
 
-For example, to create a 2-m ensemble mean temperature plot using ensemble data initialized at 2017061900, valid at 2017061912:
+For example, to plot max precipitation accumulation with contours of 500hPa mean height, and 
+mean 500hPa wind barbs from forecast hours [12, 18]:
 
-make_webplot.py -d=2017070500 -f=t2_mean -tr=12 -t='2-m Ensemble Mean Temperature (F)'
+```
+python webplot.py 20170518 --fill precipacc/max --fhr 12 18 --contour hgt500/mean \
+    --barb wind500/mean --title 'Max precip acc, mean 500hPa hgt [m] and wind barbs [kt]' \
+    --idir /glade/campaign/mmm/parc/schwartz/MPAS_JEDI/15-3km_mesh/cold_start \
+    --init_file /glade/scratch/schwartz/MPAS/15-3km_mesh/mpas_init/2019042518/init.nc \
+    --mesh 15-3km_mesh
+```
 
-The package uses basemap objects that are stored as pickle files. These should be created with the saveNewMap() function in webplot.py before attempting to plot.
+MPAS initialization time provided as first argument.
+
+```
+usage: webplot.py [-h] [--autolevels] [-b BARB] [-c CONTOUR] [-con] [-d] [--domain {CONUS,NA,SGP,NGP,CGP,SW,NW,SE,NE,MATL}] [--ENS_SIZE ENS_SIZE] [-f FILL] [--fhr FHR [FHR ...]] [--idir IDIR]
+                  [--init_file INIT_FILE] [--meshstr MESHSTR] [--nbarbs NBARBS] [--nlon_max NLON_MAX] [--nlat_max NLAT_MAX] [--sigma SIGMA] [-t TITLE]
+                  date
+
+Web plotting script for NCAR ensemble
+
+positional arguments:
+  date                  model initialization time
+
+options:
+  -h, --help            show this help message and exit
+  --autolevels          use min/max to determine levels for plot (default: False)
+  -b BARB, --barb BARB  barb field (FIELD_PRODUCT_THRESH) (default: None)
+  -c CONTOUR, --contour CONTOUR
+                        contour field (FIELD_PRODUCT_THRESH) (default: None)
+  -con, --convert       run final image through imagemagick (default: True)
+  -d, --debug           turn on debugging (default: False)
+  --domain {CONUS,NA,SGP,NGP,CGP,SW,NW,SE,NE,MATL}
+                        domain of plot (default: CONUS)
+  --ENS_SIZE ENS_SIZE   number of members in ensemble (default: 1)
+  -f FILL, --fill FILL  fill field (FIELD_PRODUCT_THRESH), FIELD options:precip,precip-24hr,precip-48hr,precipacc,sbcape,mlcape,mucape,sbcinh,mlcinh,pwat,t2,t2depart,t2-
+                        0c,mslp,td2,td2depart,thetae,thetapv,rh2m,pblh,hmuh,hmneguh,hmuh03,hmuh01,rvort1,sspf,hmup,hmdn,hmwind,hmgrp,cref,ref1km,srh3,srh1,shr06mag,shr01mag,zlfc,zlcl,ltg1,ltg2,ltg3,olrtoa,thck1000-
+                        500,thck1000-
+                        850,hgt200,hgt250,hgt300,hgt500,hgt700,hgt850,hgt925,speed200,speed250,speed300,speed500,speed700,speed850,speed925,temp200,temp250,temp300,temp500,temp700,temp850,temp925,td500,td700,td850,
+                        td925,rh300,rh500,rh700,rh850,rh925,pvort320k,speed10m,speed10m-
+                        tc,stp,uhratio,crefuh,wind10m,windsfc,wind1km,wind6km,windpv,shr06,shr01,bunkers,wind200,wind250,wind300,wind500,wind700,wind850,wind925,vort500,vort700,vort850,vortpv PRODUCT may be one of
+                        [max,maxstamp,min,mean,meanstamp,prob,neprob,problt,neproblt,paintball,stamp,spaghetti] (default: None)
+  --fhr FHR [FHR ...]   list of forecast hours (default: [12])
+  --idir IDIR           path to model output (default: /glade/campaign/mmm/parc/schwartz/MPAS_ensemble_paper)
+  --init_file INIT_FILE
+                        path to file with lat/lon/area of mesh (default: None)
+  --meshstr MESHSTR     mesh id. used to prefix output file and pickle file (default: None)
+  --nbarbs NBARBS       max barbs in one dimension (default: 32)
+  --nlon_max NLON_MAX   max pts in longitude dimension (default: 1500)
+  --nlat_max NLAT_MAX   max pts in latitude dimension (default: 1500)
+  --sigma SIGMA         smooth probabilities using gaussian smoother (default: 2)
+  -t TITLE, --title TITLE
+                        title for plot (default: None)
+```
+
+### Installation
+
+Create conda environment described in environment.yaml.
+
+Use f2py3 to install vert2cell.

environment.yaml

@@ -0,0 +1,12 @@
+name: webplot
+channels:
+  - conda-forge
+  - defaults
+dependencies:
+  - pygrib
+  - scipy
+  - xarray
+  - netcdf4
+  - dask
+  - metpy
+prefix: /glade/u/home/ahijevyc/miniconda3/envs/webplot

interp.py

@@ -0,0 +1,152 @@
+import argparse
+import logging
+import numpy as np
+import pandas as pd
+import pdb
+import pickle
+import pygrib
+import os
+from scipy.spatial import KDTree, Delaunay
+import sys
+import xarray
+
+logging.basicConfig(format='%(asctime)s - %(message)s', level=logging.INFO)
+
+class regrid:
+    '''
+    Interpolate MPAS mesh to same grid as an example grib file.
+    Output netCDF and grib.
+    '''
+    def __init__(self):
+        self.args = parseargs()
+        idate = pd.to_datetime(self.args.date)
+        idir = self.args.idir
+        var = self.args.var
+
+        # Source mesh
+        path = self.args.init_file
+        logging.info(f"get mpas mesh file {path}")
+        mpas_mesh = xarray.open_dataset(path)
+        lonCell = mpas_mesh['lonCell']
+        lonCell = np.degrees(lonCell) #radians to degrees
+        lonCell[lonCell>=180] -= 360
+        mpas_mesh["lonCell"] = lonCell
+        mpas_mesh["latCell"] = np.degrees(mpas_mesh["latCell"]) #radians to degrees
+
+        grb_with_dest_mesh = self.args.grb_with_dest_mesh
+        with pygrib.open(grb_with_dest_mesh) as f:
+            grb = f.readline()
+            dest_latlons = grb.latlons()
+        grb.dataDate = int(idate.strftime("%Y%m%d"))
+        grb.dataTime = int(idate.strftime("%H%M"))
+        # destination lats and lons
+        dest_lats, dest_lons = dest_latlons
+
+        # weights and vertices for regridding 
+        pk_file = os.path.join(os.getenv('TMPDIR'), f"{self.args.meshstr}.pk")
+        if not os.path.exists(pk_file):
+            saveNewMap(mpas_mesh, dest_latlons, pk_file)
+        logging.info(f"load {pk_file}")
+        (ibox, vtx, wts) = pickle.load(open(pk_file, 'rb'))
+
+        # Source data
+        valid_times = [idate + pd.Timedelta(fhr, unit="hour") for fhr in self.args.fhr]
+        ifiles = [  os.path.join(idir,
+                    idate.strftime("%Y%m%d%H"),
+                    f'diag.{valid_time.strftime("%Y-%m-%d_%H.%M.%S")}.nc')
+                    for valid_time in valid_times ]
+
+        nt = len(ifiles) # used to reshape output
+        logging.info(f"open {len(ifiles)} files")
+        data = xarray.open_mfdataset(ifiles, concat_dim="Time", combine="nested")
+        data = data.assign_coords(Time=valid_times)
+        data = data[var].isel(nCells=ibox)
+        # dot product of nCells x 3 arrays vtx and wtx, preserving time dimension
+        out = np.einsum('tnj,nj->tn', np.take(data.values, vtx, axis=1), wts)
+        out = np.reshape(out, (nt,*dest_lats.shape))
+        out = xarray.DataArray(data=out, name=var, 
+                coords = dict(
+                    Time=data.Time,
+                    lat=(["y","x"], dest_lats), 
+                    lon=(["y","x"], dest_lons)), 
+                dims=["Time","y","x"], attrs=data.attrs)
+
+        ogrb = f'{idate.strftime("%Y%m%d_%H")}.grb2'
+        logging.info(f"write {ogrb}")
+        with open(ogrb, 'wb') as f:
+            for i, fhr in enumerate(self.args.fhr):
+                grb.values = out.values[i]
+                grb["forecastTime"] = fhr
+                f.write(grb.tostring())
+
+        ocdf = f'{idate.strftime("%Y%m%d_%H")}.nc'
+        logging.info(f"to_netcdf {ocdf}")
+        xarray.concat(out, pd.Index(self.args.fhr, name="forecast_hour")).to_netcdf(ocdf)
+
+
+def parseargs():
+    '''Parse arguments and argparse Namespace'''
+
+    parser = argparse.ArgumentParser(description='regrid MPAS',
+            formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    parser.add_argument('date', help='model initialization time')
+    parser.add_argument('init_file', help='path to init.nc file with latCell, lonCell of mesh')
+    parser.add_argument('-d', '--debug', action='store_true', help='turn on debugging')
+    parser.add_argument('--ENS_SIZE', type=int, default=1, help='number of members in ensemble')
+    parser.add_argument('--fhr', nargs='+', type=float, default=[12], help='list of forecast hours')
+    today = pd.Timestamp.now().strftime('%Y%m%d00') 
+    parser.add_argument('--grb_with_dest_mesh', default = f'/glade/scratch/sobash/HRRR/{today}/hrrr.t00z.wrfsfcf00.grib2',
+            help = 'grib file with destination mesh')
+    parser.add_argument('--idir', default='/glade/campaign/mmm/parc/schwartz/MPAS_ensemble_paper',
+            help="path to model output")
+    parser.add_argument('--meshstr', default="HRRR", help=f'mesh id. used to prefix pickle file')
+    parser.add_argument('--var', default="refl10cm_1km_max", help=f'name of variable to regrid')
+
+    args = parser.parse_args()
+    if args.debug:
+        logging.getLogger().setLevel(logging.DEBUG)
+
+    assert os.path.exists(args.init_file), "--init_file must be path to file with latCell, lonCell, etc."
+
+    return args
+
+def saveNewMap(mpas_mesh, dest_latlons, pk_file):
+    logging.info(f"saveNewMap: pk_file={pk_file}")
+
+    # destination lats and lons
+    dest_lats, dest_lons = dest_latlons
+
+    # Source lats and lons
+    latCell = mpas_mesh.latCell
+    lonCell = mpas_mesh.lonCell
+    # To save time, triangulate near destination mesh
+    # ibox: subscripts. speed up triangulation in interp_weights
+    tree = KDTree(np.c_[dest_lons.ravel(), dest_lats.ravel()])
+    dd, ii = tree.query(np.c_[lonCell, latCell], distance_upper_bound=0.25)
+    ibox = np.isfinite(dd)
+    latCell = latCell[ibox]
+    lonCell = lonCell[ibox]
+
+    logging.info(f"saveNewMap: triangulate {len(lonCell)} pts")
+    vtx, wts = interp_weights(np.vstack((lonCell,latCell)).T,np.vstack((dest_lons.flatten(), dest_lats.flatten())).T)
+    assert wts.min() >= 0, 'got negative weight'
+
+    pickle.dump((ibox,vtx,wts), open(pk_file, 'wb'))
+
+# from https://stackoverflow.com/questions/20915502/speedup-scipy-griddata-for-multiple-interpolations-between-two-irregular-grids
+def interp_weights(xyz, uvw):
+    logging.info("interp_weights: Delaunay")
+    tri = Delaunay(xyz)
+    logging.info("interp_weights: find_simplex")
+    simplex = tri.find_simplex(uvw)
+    logging.info("interp_weights: vertices")
+    vertices = np.take(tri.simplices, simplex, axis=0)
+    temp = np.take(tri.transform, simplex, axis=0)
+    d = xyz.shape[1]
+    delta = uvw - temp[:, d]
+    logging.info("interp_weights: bary")
+    bary = np.einsum('njk,nk->nj', temp[:, :d, :], delta)
+    return vertices, np.hstack((bary, 1 - bary.sum(axis=1, keepdims=True)))
+
+if __name__ == "__main__":
+    regrid()

mpas.py

@@ -0,0 +1,11 @@
+from fieldinfo import fieldinfo,readNCLcm
+import numpy as np
+
+# fieldinfo should have been imported from fieldinfo module.
+# Copy fieldinfo dictionary for MPAS. Change some fnames.
+fieldinfo['td2']['fname'] = ['surface_dewpoint']
+fieldinfo['td2depart']['fname'] = ['surface_dewpoint']
+for plev in ['500', '700', '850']:
+    fieldinfo['vort'+plev] = {'levels' : np.array([0,9,12,15,18,21,24,27,30,33])*1e-5, 'cmap': readNCLcm('prcp_1'), 'fname': ['vorticity_'+plev+'hPa']}
+fieldinfo['vortpv']        = {'levels' : [0,9,12,15,18,21,24,27,30,33], 'cmap': readNCLcm('prcp_1'), 'fname': ['vort_pv']}
+fieldinfo['stp']['fname'] = ['sbcape','lcl','srh_0_1km','uzonal_6km','umeridional_6km','uzonal_surface','umeridional_surface']

mpas_vort_cell.f90

@@ -0,0 +1,25 @@
+subroutine vert2cell(nEdgesOnCell, verticesOnCell, maxEdges, nVertLevels, nCells, nVertices, fieldv, fieldc)
+implicit none
+
+integer, intent(in) :: maxEdges, nVertLevels, nCells, nVertices
+!f2py required :: maxEdges, nVertLevels, nCells, nVertices
+integer, intent(in) :: nEdgesOnCell(nCells)
+integer, intent(in) :: verticesOnCell(maxEdges,nCells)
+real*8, intent(in)  :: fieldv(nVertLevels,nVertices)
+real*8, intent(out) :: fieldc(nVertLevels,nCells)
+
+integer i,j,k
+real*8 factor
+
+do k=1,nVertLevels
+    do i=1,nCells
+        factor = 1./nEdgesOnCell(i)
+        fieldc(k,i) = 0.
+        do j=1,nEdgesOnCell(i)
+            fieldc(k,i)=fieldc(k,i) + fieldv(k,verticesOnCell(j,i))
+        end do
+            fieldc(k,i) = factor*fieldc(k,i)
+    end do
+end do
+return
+end subroutine vert2cell

ncar.png

@@ -0,0 +1 @@
+/glade/work/ahijevyc/share/rt_ensemble/python_scripts/ncar.png