From 6564b3c0e3badb098bf7ebc45151c0f1ecb1a7ba Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 9 Mar 2018 11:40:33 -0700 Subject: [PATCH 1/9] Updated pyf file for new signutures --- fortran/wrffortran.pyf | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/fortran/wrffortran.pyf b/fortran/wrffortran.pyf index 7211ff6..c4d9183 100644 --- a/fortran/wrffortran.pyf +++ b/fortran/wrffortran.pyf @@ -231,7 +231,7 @@ python module _wrffortran ! in integer, optional,intent(in),check(shape(prs,0)==mkzh),depend(prs) :: mkzh=shape(prs,0) integer intent(in) :: ter_follow end subroutine dpfcalc - subroutine dcapecalc3d(prs,tmk,qvp,ght,ter,sfp,cape,cin,cmsg,mix,mjy,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 + subroutine dcapecalc3d(prs,tmk,qvp,ght,ter,sfp,cape,cin,prsf,prs_new,tmk_new,qvp_new,ght_new,cmsg,mix,mjy,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 threadsafe use wrf_constants, only: tlclc2,gamma,tlclc1,rgasmd,tlclc4,g,tlclc3,thtecon3,eps,rd,cpmd,celkel,gammamd,eslcon2,eslcon1,cp,thtecon1,ezero,thtecon2 real(kind=8) dimension(mix,mjy,mkzh),intent(in) :: prs @@ -242,6 +242,11 @@ python module _wrffortran ! in real(kind=8) dimension(mix,mjy),intent(in),depend(mix,mjy) :: sfp real(kind=8) dimension(mix,mjy,mkzh),intent(out,in),depend(mix,mjy,mkzh) :: cape real(kind=8) dimension(mix,mjy,mkzh),intent(out,in),depend(mix,mjy,mkzh) :: cin + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: prsf + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: prs_new + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: tmk_new + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: qvp_new + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: ght_new real(kind=8) intent(in) :: cmsg integer, optional,intent(in),check(shape(prs,0)==mix),depend(prs) :: mix=shape(prs,0) integer, optional,intent(in),check(shape(prs,1)==mjy),depend(prs) :: mjy=shape(prs,1) @@ -251,7 +256,7 @@ python module _wrffortran ! in integer intent(inout) :: errstat character*(*) intent(inout) :: errmsg end subroutine dcapecalc3d - subroutine dcapecalc2d(prs,tmk,qvp,ght,ter,sfp,cape,cin,cmsg,mix,mjy,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 + subroutine dcapecalc2d(prs,tmk,qvp,ght,ter,sfp,cape,cin,prsf,prs_new,tmk_new,qvp_new,ght_new,cmsg,mix,mjy,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 threadsafe use wrf_constants, only: tlclc2,gamma,tlclc1,rgasmd,tlclc4,g,tlclc3,thtecon3,eps,rd,cpmd,celkel,gammamd,eslcon2,eslcon1,cp,thtecon1,ezero,thtecon2 real(kind=8) dimension(mix,mjy,mkzh),intent(in) :: prs @@ -262,6 +267,11 @@ python module _wrffortran ! in real(kind=8) dimension(mix,mjy),intent(in),depend(mix,mjy) :: sfp real(kind=8) dimension(mix,mjy,mkzh),intent(out,in),depend(mix,mjy,mkzh) :: cape real(kind=8) dimension(mix,mjy,mkzh),intent(out,in),depend(mix,mjy,mkzh) :: cin + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: prsf + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: prs_new + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: tmk_new + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: qvp_new + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: ght_new real(kind=8) intent(in) :: cmsg integer, optional,intent(in),check(shape(prs,0)==mix),depend(prs) :: mix=shape(prs,0) integer, optional,intent(in),check(shape(prs,1)==mjy),depend(prs) :: mjy=shape(prs,1) @@ -349,7 +359,7 @@ python module _wrffortran ! in real(kind=8), parameter,optional,depend(cp,rd) :: gamma=0.285714285714 real(kind=8), parameter,optional,depend(expon,rd,ussalr,g) :: exponi=5.25864379523 end module wrf_constants - subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt,haveqci,nz,ns,ew) ! in :_wrffortran:wrf_fctt.f90 + subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt,pf,haveqci,nz,ns,ew) ! in :_wrffortran:wrf_fctt.f90 threadsafe use wrf_constants, only: g,eps,rd,ussalr,abscoefi,abscoef,celkel real(kind=8) dimension(ew,ns,nz),intent(in) :: prs @@ -360,6 +370,7 @@ python module _wrffortran ! in real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: ght real(kind=8) dimension(ew,ns),intent(in),depend(ew,ns) :: ter real(kind=8) dimension(ew,ns),intent(out,in),depend(ew,ns) :: ctt + real(kind=8) dimension(ew,ns,nz),intent(inout),depend(ew,ns,nz) :: pf integer intent(in) :: haveqci integer, optional,intent(in),check(shape(prs,2)==nz),depend(prs) :: nz=shape(prs,2) integer, optional,intent(in),check(shape(prs,1)==ns),depend(prs) :: ns=shape(prs,1) @@ -730,7 +741,7 @@ python module _wrffortran ! in integer intent(inout) :: errstat real(kind=8) :: wrf_intrp_value end function wrf_intrp_value - subroutine wrf_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,sfp,smsfp,vcarray,interp_levels,numlevels,icase,ew,ns,nz,extrap,vcor,logp,rmsg,errstat,errmsg) ! in :_wrffortran:wrf_vinterp.f90 + subroutine wrf_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,sfp,smsfp,vcarray,interp_levels,numlevels,icase,ew,ns,nz,extrap,vcor,logp,tempout,rmsg,errstat,errmsg) ! in :_wrffortran:wrf_vinterp.f90 threadsafe use wrf_constants, only: tlclc2,tlclc3,sclht,tlclc4,thtecon3,eps,tlclc1,celkel,exponi,gammamd,ussalr,expon,thtecon1,algerr,gamma,thtecon2 real(kind=8) dimension(ew,ns,nz),intent(in) :: datain @@ -752,6 +763,7 @@ python module _wrffortran ! in integer intent(in) :: extrap integer intent(in) :: vcor integer intent(in) :: logp + real(kind=8) dimension(ew,ns),intent(inout),depend(ew,ns) :: tempout real(kind=8) intent(in) :: rmsg integer intent(inout) :: errstat character*(*) intent(inout) :: errmsg From 71fe678e916ab632ad8fb2b99ddc19a34bf7aa9d Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 9 Mar 2018 11:42:19 -0700 Subject: [PATCH 2/9] Moved work arrays out of fortran. Fixes #47 --- fortran/rip_cape.f90 | 28 ++++++++++++++++--------- fortran/wrf_fctt.f90 | 6 ++++-- fortran/wrf_vinterp.f90 | 6 ++++-- src/wrf/extension.py | 45 +++++++++++++++++++++++++++++------------ 4 files changed, 58 insertions(+), 27 deletions(-) diff --git a/fortran/rip_cape.f90 b/fortran/rip_cape.f90 index 5e9cdac..1fc3b23 100644 --- a/fortran/rip_cape.f90 +++ b/fortran/rip_cape.f90 @@ -272,6 +272,7 @@ END SUBROUTINE DPFCALC ! NCLFORTSTART SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& + prsf, prs_new, tmk_new, qvp_new, ght_new,& cmsg,mix,mjy,mkzh,ter_follow,& psafile, errstat, errmsg) USE wrf_constants, ONLY : CELKEL, G, EZERO, ESLCON1, ESLCON2, & @@ -293,6 +294,14 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) ::sfp REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cape REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cin + + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prsf + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prs_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: tmk_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: qvp_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: ght_new + + REAL(KIND=8), INTENT(IN) :: cmsg CHARACTER(LEN=*), INTENT(IN) :: psafile INTEGER, INTENT(INOUT) :: errstat @@ -308,15 +317,10 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat REAL(KIND=8) :: benamin, dz REAL(KIND=8), DIMENSION(150) :: buoy, zrel, benaccum - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prsf REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs REAL(KIND=8), DIMENSION(150,150) :: psaditmk LOGICAL :: elfound - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prs_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: tmk_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: qvp_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: ght_new ! To remove compiler warnings tmkpari = 0 @@ -597,6 +601,7 @@ END SUBROUTINE DCAPECALC3D ! NCLFORTSTART SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& + prsf, prs_new, tmk_new, qvp_new, ght_new,& cmsg,mix,mjy,mkzh,ter_follow,& psafile, errstat, errmsg) USE wrf_constants, ONLY : CELKEL, G, EZERO, ESLCON1, ESLCON2, & @@ -618,6 +623,13 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) ::sfp REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cape REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cin + + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prsf + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prs_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: tmk_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: qvp_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: ght_new + REAL(KIND=8), INTENT(IN) :: cmsg CHARACTER(LEN=*), INTENT(IN) :: psafile INTEGER, INTENT(INOUT) :: errstat @@ -635,15 +647,11 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat REAL(KIND=8) :: benamin, dz, pup, pdn REAL(KIND=8), DIMENSION(150) :: buoy, zrel, benaccum - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prsf REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs REAL(KIND=8), DIMENSION(150,150) :: psaditmk LOGICAL :: elfound REAL(KIND=8), DIMENSION(mkzh) :: eth_temp - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prs_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: tmk_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: qvp_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: ght_new + ! To remove compiler warnings errstat = 0 diff --git a/fortran/wrf_fctt.f90 b/fortran/wrf_fctt.f90 index 2ef7ecb..0be18ae 100644 --- a/fortran/wrf_fctt.f90 +++ b/fortran/wrf_fctt.f90 @@ -1,5 +1,5 @@ !NCLFORTSTART -SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew) +SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns, ew) USE wrf_constants, ONLY : EPS, USSALR, RD, G, ABSCOEFI, ABSCOEF, CELKEL IMPLICIT NONE @@ -12,6 +12,8 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: ter REAL(KIND=8), DIMENSION(ew,ns), INTENT(OUT) :: ctt + REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(INOUT) :: pf + !NCLEND ! REAL(KIND=8) :: znfac(nz) @@ -20,7 +22,7 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew INTEGER i,j,k,ripk REAL(KIND=8) :: opdepthu, opdepthd, dp, arg1, fac, prsctt, ratmix REAL(KIND=8) :: arg2, agl_hgt, vt - REAL(KIND=8), DIMENSION(ew,ns,nz) :: pf + REAL(KIND=8) :: p1, p2 !$OMP PARALLEL diff --git a/fortran/wrf_vinterp.f90 b/fortran/wrf_vinterp.f90 index c34efd7..d81270b 100644 --- a/fortran/wrf_vinterp.f90 +++ b/fortran/wrf_vinterp.f90 @@ -127,7 +127,7 @@ END FUNCTION wrf_intrp_value !NCLFORTSTART SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& sfp, smsfp, vcarray, interp_levels, numlevels,& - icase, ew, ns, nz, extrap, vcor, logp, rmsg,& + icase, ew, ns, nz, extrap, vcor, logp, tempout, rmsg,& errstat, errmsg) USE wrf_constants, ONLY : ALGERR, SCLHT, EXPON, EXPONI, GAMMA, GAMMAMD, TLCLC1, & TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3, & @@ -146,6 +146,9 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& REAL(KIND=8), DIMENSION(ew,ns,numlevels), INTENT(OUT) :: dataout REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: vcarray REAL(KIND=8), DIMENSION(numlevels), INTENT(IN) :: interp_levels + + REAL(KIND=8), DIMENSION(ew,ns), INTENT(INOUT) :: tempout + REAL(KIND=8), INTENT(IN) :: rmsg INTEGER, INTENT(INOUT) :: errstat CHARACTER(LEN=*), INTENT(INOUT) :: errmsg @@ -156,7 +159,6 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& INTEGER :: i, j, k, kupper !itriv INTEGER :: ifound, isign !miy,mjx INTEGER :: log_errcnt, interp_errcnt, interp_errstat - REAL(KIND=8), DIMENSION(ew,ns) :: tempout REAL(KIND=8) :: rlevel, vlev, diff REAL(KIND=8) :: tmpvlev REAL(KIND=8) :: vcp1, vcp0, valp0, valp1 diff --git a/src/wrf/extension.py b/src/wrf/extension.py index 167d27c..5bcccfa 100755 --- a/src/wrf/extension.py +++ b/src/wrf/extension.py @@ -618,20 +618,33 @@ def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow, else: cape_routine = dcapecalc2d + # Work arrays + k_left_shape = (p_hpa.shape[2], p_hpa.shape[0], p_hpa.shape[1]) + prsf = np.empty(k_left_shape, np.float64, order="F") + prs_new = np.empty(k_left_shape, np.float64, order="F") + tmk_new = np.empty(k_left_shape, np.float64, order="F") + qvp_new = np.empty(k_left_shape, np.float64, order="F") + ght_new = np.empty(k_left_shape, np.float64, order="F") + # note that p_hpa, tk, qv, and ht have the vertical flipped result = cape_routine(p_hpa, - tk, - qv, - ht, - ter, - sfp, - capeview, - cinview, - missing, - ter_follow, - psafile, - errstat, - errmsg) + tk, + qv, + ht, + ter, + sfp, + capeview, + cinview, + prsf, + prs_new, + tmk_new, + qvp_new, + ght_new, + missing, + ter_follow, + psafile, + errstat, + errmsg) if int(errstat) != 0: raise DiagnosticError("".join(npbytes_to_str(errmsg)).strip()) @@ -766,6 +779,8 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): if outview is None: outview = np.empty_like(ter) + pf = np.empty(p_hpa.shape[0:3], np.float64, order="F") + result = wrfcttcalc(p_hpa, tk, qice, @@ -774,6 +789,7 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): ght, ter, outview, + pf, haveqci) return result @@ -858,7 +874,9 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, if outview is None: outdims = field.shape[0:2] + interp_levels.shape outview = np.empty(outdims, field.dtype, order="F") - + + tempout = np.zeros(field.shape[0:2], np.float64, order="F") + errstat = np.array(0) errmsg = np.zeros(Constants.ERRLEN, "c") @@ -877,6 +895,7 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, extrap, vcor, logp, + tempout, missing, errstat, errmsg) From e6b4f7382b4cb485aab9e8182444e78af543418e Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Wed, 14 Mar 2018 15:54:20 -0600 Subject: [PATCH 3/9] Added checks for _key=None in extract_vars. Fix typo with updraft helicity cache key. Fixes #51. --- src/wrf/g_helicity.py | 2 +- src/wrf/util.py | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/wrf/g_helicity.py b/src/wrf/g_helicity.py index 02ae780..cb3beb5 100755 --- a/src/wrf/g_helicity.py +++ b/src/wrf/g_helicity.py @@ -176,7 +176,7 @@ def get_uh(wrfin, timeidx=0, method="cat", squeeze=True, """ ncvars = extract_vars(wrfin, timeidx, ("W", "PH", "PHB", "MAPFAC_M"), - method, squeeze, cache, meta=False, _key=None) + method, squeeze, cache, meta=False, _key=_key) wstag = ncvars["W"] ph = ncvars["PH"] diff --git a/src/wrf/util.py b/src/wrf/util.py index 3aa9583..7bf9927 100644 --- a/src/wrf/util.py +++ b/src/wrf/util.py @@ -973,10 +973,12 @@ def _combine_dict(wrfdict, varname, timeidx, method, meta, _key): is_moving = is_moving_domain(wrfdict, varname, _key=_key) + _cache_key = _key[first_key] if _key is not None else None + first_array = _extract_var(wrfdict[first_key], varname, timeidx, is_moving=is_moving, method=method, squeeze=False, cache=None, meta=meta, - _key=_key[first_key]) + _key=_cache_key) # Create the output data numpy array based on the first array outdims = [numkeys] @@ -992,10 +994,11 @@ def _combine_dict(wrfdict, varname, timeidx, method, meta, _key): break else: keynames.append(key) + _cache_key = _key[key] if _key is not None else None vardata = _extract_var(wrfdict[key], varname, timeidx, is_moving=is_moving, method=method, squeeze=False, cache=None, meta=meta, - _key=_key[key]) + _key=_cache_key) if outdata.shape[1:] != vardata.shape: raise ValueError("data sequences must have the " From 196e28760045e5763f034eeb7b5196f9006083c9 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 16 Mar 2018 10:05:18 -0600 Subject: [PATCH 4/9] Added th alias for theta. Fixes #52. --- src/wrf/routines.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/wrf/routines.py b/src/wrf/routines.py index 5442956..f3665fe 100644 --- a/src/wrf/routines.py +++ b/src/wrf/routines.py @@ -112,7 +112,7 @@ _VALID_KARGS = {"cape2d" : ["missing"], "wspd_wdir10" : ["units"], "uvmet_wspd_wdir" : ["units"], "uvmet10_wspd_wdir" : ["units"], - "ctt" : [], + "ctt" : ["fill_nocloud", "missing", "units"], "cloudfrac" : ["vert_type", "low_thresh", "mid_thresh", "high_thresh"], "geopt_stag" : [], @@ -138,7 +138,8 @@ _ALIASES = {"cape_2d" : "cape2d", "td2" : "dp2m", "cfrac" : "cloudfrac", "wspd_wdir_uvmet" : "uvmet_wspd_wdir", - "wspd_wdir_uvmet10" : "uvmet10_wspd_wdir" + "wspd_wdir_uvmet10" : "uvmet10_wspd_wdir", + "th" : "theta" } class ArgumentError(Exception): From cc86d68a35efb5de440426e595c4286c875d1108 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 16 Mar 2018 16:34:44 -0600 Subject: [PATCH 5/9] Fix issues related to cloud top temperature. Fixed indexing bug. Fixed incorrect computation of optical depth when cloud ice is not available. Users can use fill values for cloud free areas. Users can now specify the optical depth threshold that triggers the calculation. Fixes #45. --- doc/source/_templates/product_table.txt | 6 +- doc/source/new.rst | 20 +++ fortran/wrf_fctt.f90 | 54 +++--- fortran/wrffortran.pyf | 5 +- src/wrf/computation.py | 34 +++- src/wrf/extension.py | 10 +- src/wrf/g_ctt.py | 38 ++++- src/wrf/metadecorators.py | 4 +- src/wrf/routines.py | 2 +- src/wrf/version.py | 2 +- test/ci_tests/make_test_file.py | 2 +- test/comp_utest.py | 4 +- test/ipynb/Test_CTT.ipynb | 214 ++++++++++++++++++++++++ test/ipynb/Test_CTT_Prod.ipynb | 214 ++++++++++++++++++++++++ 14 files changed, 565 insertions(+), 44 deletions(-) create mode 100644 test/ipynb/Test_CTT.ipynb create mode 100644 test/ipynb/Test_CTT_Prod.ipynb diff --git a/doc/source/_templates/product_table.txt b/doc/source/_templates/product_table.txt index c49c3ac..b2ee337 100644 --- a/doc/source/_templates/product_table.txt +++ b/doc/source/_templates/product_table.txt @@ -13,11 +13,13 @@ +--------------------+---------------------------------------------------------------+-----------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------+ | cape_3d | 3D cape and cin | J kg-1 | **missing** (float): Fill value for output only | +--------------------+---------------------------------------------------------------+-----------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------+ -| ctt | Cloud Top Temperature | degC | **units** (str) : Set to desired units. Default is *'degC'*. | +| ctt | Cloud Top Temperature | degC | **fill_nocloud** (boolean): Set to True to use fill values for cloud free regions rather than surface temperature. Default is *False*. | | | | | | -| | | K | | +| | | K | **missing** (float): The fill value to use when *fill_nocloud* is True. | | | | | | +| | | | **opt_thresh** (float): The optical depth required to trigger the cloud top temperature calculation. Default is 1.0. | | | | degF | | +| | | | **units** (str) : Set to desired units. Default is *'degC'*. | +--------------------+---------------------------------------------------------------+-----------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------+ | cloudfrac | Cloud Fraction | % | **vert_type** (str): The vertical coordinate type for the cloud thresholds. Must be 'height_agl', 'height_msl', or 'pres'. Default is 'height_agl'. | | | | | | diff --git a/doc/source/new.rst b/doc/source/new.rst index 2962534..c07e5a5 100644 --- a/doc/source/new.rst +++ b/doc/source/new.rst @@ -4,6 +4,26 @@ What's New Releases ------------- +v1.1.3 +^^^^^^^^^^^^^^ + +- Release 1.1.3 +- Fixed/Enhanced the cloud top temperature diagnostic. + - Optical depth was not being calculated correctly when + cloud ice mixing ratio was not available. + - Fixed an indexing bug that caused crashes on Windows, but should have been + crashing on all platforms. + - Users can now specify if they want cloud free regions to use fill values, + rather than the default behavior of using the surface temperature. + - Users can now specify the optical depth required to trigger the cloud + top temperature calculation. However, the default value of 1.0 should be + sufficient for most users. +- Added 'th' alias for the theta product. +- Fixed a crash issue related to updraft helicity when a dictionary is + used as the input. +- The cape_2d diagnostic can now work with a single column of data, just like + cape_3d. + v1.1.2 ^^^^^^^^^^^^^^ diff --git a/fortran/wrf_fctt.f90 b/fortran/wrf_fctt.f90 index 0be18ae..2e61632 100644 --- a/fortran/wrf_fctt.f90 +++ b/fortran/wrf_fctt.f90 @@ -1,5 +1,6 @@ !NCLFORTSTART -SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns, ew) +SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci,& + fill_nocloud, missing, opt_thresh, nz, ns, ew) USE wrf_constants, ONLY : EPS, USSALR, RD, G, ABSCOEFI, ABSCOEF, CELKEL IMPLICIT NONE @@ -7,13 +8,15 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns !f2py threadsafe !f2py intent(in,out) :: ctt - INTEGER, INTENT(IN) :: nz, ns, ew, haveqci + INTEGER, INTENT(IN) :: nz, ns, ew, haveqci, fill_nocloud REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: ght, prs, tk, qci, qcw, qvp REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: ter REAL(KIND=8), DIMENSION(ew,ns), INTENT(OUT) :: ctt - + REAL(KIND=8), INTENT(IN) :: missing + REAL(KIND=8), INTENT(IN) :: opt_thresh REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(INOUT) :: pf + !NCLEND ! REAL(KIND=8) :: znfac(nz) @@ -60,12 +63,12 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns DO i=1,ew opdepthd = 0.D0 k = 0 - prsctt = 0 + prsctt = -1 ! Integrate downward from model top, calculating path at full ! model vertical levels. - DO k=1, nz + DO k=2,nz opdepthu = opdepthd ripk = nz - k + 1 @@ -76,21 +79,23 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns END IF IF (haveqci .EQ. 0) then - IF (tk(i,j,k) .LT. CELKEL) then + IF (tk(i,j,ripk) .LT. CELKEL) then ! Note: abscoefi is m**2/g, qcw is g/kg, so no convrsion needed - opdepthd = opdepthu + ABSCOEFI*qcw(i,j,k) * dp/G + opdepthd = opdepthu + ABSCOEFI*qcw(i,j,ripk) * dp/G ELSE - opdepthd = opdepthu + ABSCOEF*qcw(i,j,k) * dp/G + opdepthd = opdepthu + ABSCOEF*qcw(i,j,ripk) * dp/G END IF ELSE opdepthd = opdepthd + (ABSCOEF*qcw(i,j,ripk) + ABSCOEFI*qci(i,j,ripk))*dp/G END IF - IF (opdepthd .LT. 1. .AND. k .LT. nz) THEN + IF (opdepthd .LT. opt_thresh .AND. k .LT. nz) THEN CYCLE - ELSE IF (opdepthd .LT. 1. .AND. k .EQ. nz) THEN - prsctt = prs(i,j,1) + ELSE IF (opdepthd .LT. opt_thresh .AND. k .EQ. nz) THEN + IF (fill_nocloud .EQ. 0) THEN + prsctt = prs(i,j,1) + ENDIF EXIT ELSE fac = (1. - opdepthu)/(opdepthd - opdepthu) @@ -100,17 +105,22 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns END IF END DO - DO k=2,nz - ripk = nz - k + 1 - p1 = prs(i,j,ripk+1) - p2 = prs(i,j,ripk) - IF (prsctt .GE. p1 .AND. prsctt .LE. p2) THEN - fac = (prsctt - p1)/(p2 - p1) - arg1 = fac*(tk(i,j,ripk) - tk(i,j,ripk+1)) - CELKEL - ctt(i,j) = tk(i,j,ripk+1) + arg1 - EXIT - END IF - END DO + ! prsctt should only be 0 if fill values are used + IF (prsctt .GT. -1) THEN + DO k=2,nz + ripk = nz - k + 1 + p1 = prs(i,j,ripk+1) + p2 = prs(i,j,ripk) + IF (prsctt .GE. p1 .AND. prsctt .LE. p2) THEN + fac = (prsctt - p1)/(p2 - p1) + arg1 = fac*(tk(i,j,ripk) - tk(i,j,ripk+1)) - CELKEL + ctt(i,j) = tk(i,j,ripk+1) + arg1 + EXIT + END IF + END DO + ELSE + ctt(i,j) = missing + END IF END DO END DO !$OMP END DO diff --git a/fortran/wrffortran.pyf b/fortran/wrffortran.pyf index c4d9183..be27452 100644 --- a/fortran/wrffortran.pyf +++ b/fortran/wrffortran.pyf @@ -359,7 +359,7 @@ python module _wrffortran ! in real(kind=8), parameter,optional,depend(cp,rd) :: gamma=0.285714285714 real(kind=8), parameter,optional,depend(expon,rd,ussalr,g) :: exponi=5.25864379523 end module wrf_constants - subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt,pf,haveqci,nz,ns,ew) ! in :_wrffortran:wrf_fctt.f90 + subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt,pf,haveqci,fill_nocloud,missing,opt_thresh,nz,ns,ew) ! in :_wrffortran:wrf_fctt.f90 threadsafe use wrf_constants, only: g,eps,rd,ussalr,abscoefi,abscoef,celkel real(kind=8) dimension(ew,ns,nz),intent(in) :: prs @@ -372,6 +372,9 @@ python module _wrffortran ! in real(kind=8) dimension(ew,ns),intent(out,in),depend(ew,ns) :: ctt real(kind=8) dimension(ew,ns,nz),intent(inout),depend(ew,ns,nz) :: pf integer intent(in) :: haveqci + integer intent(in) :: fill_nocloud + real(kind=8) intent(in) :: missing + real(kind=8) intent(in) :: opt_thresh integer, optional,intent(in),check(shape(prs,2)==nz),depend(prs) :: nz=shape(prs,2) integer, optional,intent(in),check(shape(prs,1)==ns),depend(prs) :: ns=shape(prs,1) integer, optional,intent(in),check(shape(prs,0)==ew),depend(prs) :: ew=shape(prs,0) diff --git a/src/wrf/computation.py b/src/wrf/computation.py index 3973709..164d905 100644 --- a/src/wrf/computation.py +++ b/src/wrf/computation.py @@ -1014,6 +1014,11 @@ def cloudfrac(vert, relh, vert_inc_w_height, low_thresh, mid_thresh, high_thresh (:obj:`float`): The bottom vertical threshold for what is considered a high cloud. + + missing (:obj:`float:`, optional): The fill value to use for areas + where the surface is higher than the cloud threshold level + (e.g. mountains). Default is + :data:`wrf.default_fill(numpy.float64)`. meta (:obj:`bool`): Set to False to disable metadata and return :class:`numpy.ndarray` instead of @@ -1047,8 +1052,9 @@ def cloudfrac(vert, relh, vert_inc_w_height, low_thresh, mid_thresh, @set_alg_metadata(2, "pres_hpa", refvarndims=3, description="cloud top temperature") @convert_units("temp", "c") -def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True, - units="degC"): +def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, + fill_nocloud=False, missing=default_fill(np.float64), + opt_thresh=1.0, meta=True, units="degC"): """Return the cloud top temperature. This is the raw computational algorithm and does not extract any variables @@ -1094,6 +1100,23 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True, qice (:class:`xarray.DataArray` or :class:`numpy.ndarray`, optional): Ice mixing ratio in [kg/kg] with the same dimensionality as *pres_hpa*. + + fill_nocloud (:obj:`bool`, optional): Set to True to use fill values in + regions where clouds are not detected (optical depth less than 1). + Otherwise, the output will contain the surface temperature for + areas without clouds. Default is False. + + missing (:obj:`float`, optional): The fill value to use for areas + where no clouds are detected. Only used if *fill_nocloud* is + True. Default is + :data:`wrf.default_fill(numpy.float64)`. + + opt_thresh (:obj:`float`, optional): The required amount of optical + depth (looking from top down) required to trigger a cloud top + temperature calculation. Values less than this threshold will be + treated as no cloud areas. For thin cirrus, a value of .1 might be + appropriate. For large CB clouds, 1000.0 might be appropriate. + Default is 1.0. meta (:obj:`bool`): Set to False to disable metadata and return :class:`numpy.ndarray` instead of @@ -1128,8 +1151,13 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True, haveqci = 0 else: haveqci = 1 if qice.any() else 0 + + _fill_nocloud = 1 if fill_nocloud else 0 + + ctt = _ctt(pres_hpa, tkel, qice, qcld, qv, height, terrain, haveqci, + _fill_nocloud, missing, opt_thresh) - return _ctt(pres_hpa, tkel, qice, qcld, qv, height, terrain, haveqci) + return ma.masked_values(ctt, missing) diff --git a/src/wrf/extension.py b/src/wrf/extension.py index 5bcccfa..42c4734 100755 --- a/src/wrf/extension.py +++ b/src/wrf/extension.py @@ -767,10 +767,11 @@ def _xytoll(map_proj, truelat1, truelat2, stdlon, lat1, lon1, @check_args(0, 3, (3,3,3,3,3,3,2)) -@left_iteration(3, 2, ref_var_idx=0, ignore_args=(7,)) +@left_iteration(3, 2, ref_var_idx=0, ignore_args=(7,8,9,10)) @cast_type(arg_idxs=(0,1,2,3,4,5,6)) @extract_and_transpose() -def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): +def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, fill_nocloud, + missing, opt_thresh, outview=None): """Wrapper for wrfcttcalc. Located in wrf_fctt.f90. @@ -790,7 +791,10 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): ter, outview, pf, - haveqci) + haveqci, + fill_nocloud, + missing, + opt_thresh) return result diff --git a/src/wrf/g_ctt.py b/src/wrf/g_ctt.py index c3dff09..1c59c01 100644 --- a/src/wrf/g_ctt.py +++ b/src/wrf/g_ctt.py @@ -1,11 +1,12 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) -import numpy as n +import numpy as np +import numpy.ma as ma #from .extension import computectt, computetk from .extension import _ctt, _tk -from .constants import Constants, ConversionFactors +from .constants import Constants, ConversionFactors, default_fill from .destag import destagger from .decorators import convert_units from .metadecorators import copy_and_set_metadata @@ -19,7 +20,8 @@ from .util import extract_vars @convert_units("temp", "c") def get_ctt(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, meta=True, _key=None, - units="degC"): + fill_nocloud=False, missing=default_fill(np.float64), + opt_thresh=1.0, units="degC"): """Return the cloud top temperature. This functions extracts the necessary variables from the NetCDF file @@ -64,6 +66,27 @@ def get_ctt(wrfin, timeidx=0, method="cat", _key (:obj:`int`, optional): A caching key. This is used for internal purposes only. Default is None. + fill_nocloud (:obj:`bool`, optional): Set to True to use fill values in + regions where clouds are not detected (optical depth less than 1). + Otherwise, the output will contain the surface temperature for + areas without clouds. Default is False. + + missing (:obj:`float`, optional): The fill value to use for areas + where no clouds are detected. Only used if *fill_nocloud* is + True. Default is + :data:`wrf.default_fill(numpy.float64)`. + + opt_thresh (:obj:`float`, optional): The amount of optical + depth (integrated from top down) required to trigger a cloud top + temperature calculation. The cloud top temperature is calculated at + the vertical level where this threshold is met. Vertical columns + with less than this threshold will be treated as cloud free areas. + In general, the larger the value is for this + threshold, the lower the altitude will be for the cloud top + temperature calculation, and therefore higher cloud top + temperature values. Default is 1.0, which should be sufficient for + most users. + units (:obj:`str`): The desired units. Refer to the :meth:`getvar` product table for a list of available units for 'ctt'. Default is 'degC'. @@ -93,7 +116,7 @@ def get_ctt(wrfin, timeidx=0, method="cat", method, squeeze, cache, meta=False, _key=_key) except KeyError: - qice = n.zeros(qv.shape, qv.dtype) + qice = np.zeros(qv.shape, qv.dtype) haveqci = 0 else: qice = icevars["QICE"] * 1000.0 #g/kg @@ -116,6 +139,9 @@ def get_ctt(wrfin, timeidx=0, method="cat", geopt_unstag = destagger(geopt, -3) ght = geopt_unstag / Constants.G - ctt = _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci) + _fill_nocloud = 1 if fill_nocloud else 0 + + ctt = _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, _fill_nocloud, + missing, opt_thresh) - return ctt + return ma.masked_values(ctt, missing) diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index 27425c1..b35d57d 100644 --- a/src/wrf/metadecorators.py +++ b/src/wrf/metadecorators.py @@ -2034,8 +2034,8 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): outcoords = {} # Left-most is always cape_cin or cape_cin_lcl_lfc if is2d: - outdims[0] = "cape_cin_lcl_lfc" - outcoords["cape_cin_lcl_lfc"] = ["cape", "cin", "lcl", "lfc"] + outdims[0] = "mcape_mcin_lcl_lfc" + outcoords["mcape_mcin_lcl_lfc"] = ["mcape", "mcin", "lcl", "lfc"] else: outdims[0] = "cape_cin" outcoords["cape_cin"] = ["cape", "cin"] diff --git a/src/wrf/routines.py b/src/wrf/routines.py index f3665fe..df6d533 100644 --- a/src/wrf/routines.py +++ b/src/wrf/routines.py @@ -112,7 +112,7 @@ _VALID_KARGS = {"cape2d" : ["missing"], "wspd_wdir10" : ["units"], "uvmet_wspd_wdir" : ["units"], "uvmet10_wspd_wdir" : ["units"], - "ctt" : ["fill_nocloud", "missing", "units"], + "ctt" : ["fill_nocloud", "missing", "opt_thresh", "units"], "cloudfrac" : ["vert_type", "low_thresh", "mid_thresh", "high_thresh"], "geopt_stag" : [], diff --git a/src/wrf/version.py b/src/wrf/version.py index e39ded1..18b2dba 100644 --- a/src/wrf/version.py +++ b/src/wrf/version.py @@ -1,2 +1,2 @@ -__version__ = "1.1.2" +__version__ = "1.1.3" diff --git a/test/ci_tests/make_test_file.py b/test/ci_tests/make_test_file.py index e6ec29b..fcdfac3 100644 --- a/test/ci_tests/make_test_file.py +++ b/test/ci_tests/make_test_file.py @@ -12,7 +12,7 @@ from wrf import (getvar, interplevel, interpline, vertcross, vinterp, py2round, VARS_TO_KEEP = ("Times", "XLAT", "XLONG", "XLAT_U", "XLAT_V", "XLONG_U", "XLONG_V", "U", "V", "W", "PH", "PHB", "T", "P", "PB", "Q2", "T2", "PSFC", "U10", "V10", "XTIME", "QVAPOR", "QCLOUD", - "QGRAUP", "QRAIN", "QSNOW", "MAPFAC_M", "MAPFAC_U", + "QGRAUP", "QRAIN", "QSNOW", "QICE", "MAPFAC_M", "MAPFAC_U", "MAPFAC_V", "F", "HGT", "RAINC", "RAINSH", "RAINNC") DIMS_TO_TRIM = ("west_east", "south_north", "bottom_top", "bottom_top_stag", diff --git a/test/comp_utest.py b/test/comp_utest.py index 55448a8..958cd75 100644 --- a/test/comp_utest.py +++ b/test/comp_utest.py @@ -609,14 +609,14 @@ if __name__ == "__main__": #varnames = ["helicity"] varnames=["avo", "pvo", "eth", "dbz", "helicity", "updraft_helicity", "omg", "pw", "rh", "slp", "td", "tk", "tv", "twb", "uvmet", - "cloudfrac"] + "cloudfrac", "ctt"] omp_set_num_threads(omp_get_num_procs()-1) omp_set_schedule(OMP_SCHED_STATIC, 0) omp_set_dynamic(False) # Turn this one off when not needed, since it's slow - #varnames += ["cape_2d", "cape_3d"] + varnames += ["cape_2d", "cape_3d"] for varname in varnames: for i,wrfnc in enumerate((NCFILE, NCGROUP)): diff --git a/test/ipynb/Test_CTT.ipynb b/test/ipynb/Test_CTT.ipynb new file mode 100644 index 0000000..533ea13 --- /dev/null +++ b/test/ipynb/Test_CTT.ipynb @@ -0,0 +1,214 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import numpy as np\n", + "from netCDF4 import Dataset\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.cm import get_cmap\n", + "import cartopy.crs as crs\n", + "from cartopy.feature import NaturalEarthFeature\n", + "\n", + "from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords\n", + "\n", + "# Open the NetCDF file\n", + "ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d01_2005-08-28_12:00:00\")\n", + "\n", + "# Get the sea level pressure\n", + "ctt = getvar(ncfile, \"ctt\", fill_nocloud=False)\n", + "slp = getvar(ncfile, \"slp\")\n", + "\n", + "\n", + "# Get the latitude and longitude points\n", + "lats, lons = latlon_coords(ctt)\n", + "\n", + "# Get the cartopy mapping object\n", + "cart_proj = get_cartopy(ctt)\n", + "\n", + "# Create a figure\n", + "fig = plt.figure(figsize=(12,9))\n", + "# Set the GeoAxes to the projection used by WRF\n", + "ax = plt.axes(projection=cart_proj)\n", + "\n", + "# Download and add the states and coastlines\n", + "states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n", + " name='admin_1_states_provinces_shp')\n", + "ax.add_feature(states, linewidth=.5)\n", + "ax.coastlines('50m', linewidth=0.8)\n", + "\n", + "# Make the contour outlines and filled contours for the smoothed sea level pressure.\n", + "levels = np.arange(-80, 20, 5)\n", + "plt.contour(to_np(lons), to_np(lats), to_np(slp), 10, colors=\"black\",\n", + " transform=crs.PlateCarree())\n", + "plt.contourf(to_np(lons), to_np(lats), to_np(ctt), levels=levels, transform=crs.PlateCarree(),\n", + " cmap=get_cmap(\"Greys\"))\n", + "\n", + "# Add a color bar\n", + "plt.colorbar(ax=ax, shrink=.88)\n", + "\n", + "# Set the map limits. Not really necessary, but used for demonstration.\n", + "ax.set_xlim(cartopy_xlim(ctt))\n", + "ax.set_ylim(cartopy_ylim(ctt))\n", + "\n", + "# Add the gridlines\n", + "ax.gridlines(color=\"black\", linestyle=\"dotted\")\n", + "\n", + "plt.title(\"Cloud Top Temperature (degC)\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import numpy as np\n", + "from netCDF4 import Dataset\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.cm import get_cmap\n", + "import cartopy.crs as crs\n", + "from cartopy.feature import NaturalEarthFeature\n", + "\n", + "from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords\n", + "\n", + "# Open the NetCDF file\n", + "ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d01_2005-08-28_12:00:00\")\n", + "\n", + "# Get the sea level pressure\n", + "ctt = getvar(ncfile, \"ctt\", fill_nocloud=True, opt_thresh=1.0)\n", + "slp = getvar(ncfile, \"slp\")\n", + "\n", + "\n", + "# Get the latitude and longitude points\n", + "lats, lons = latlon_coords(ctt)\n", + "\n", + "# Get the cartopy mapping object\n", + "cart_proj = get_cartopy(ctt)\n", + "\n", + "# Create a figure\n", + "fig = plt.figure(figsize=(12,9))\n", + "# Set the GeoAxes to the projection used by WRF\n", + "ax = plt.axes(projection=cart_proj)\n", + "\n", + "# Download and add the states and coastlines\n", + "states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n", + " name='admin_1_states_provinces_shp')\n", + "ax.add_feature(states, linewidth=.5)\n", + "ax.coastlines('50m', linewidth=0.8)\n", + "\n", + "# Make the contour outlines and filled contours for the smoothed sea level pressure.\n", + "levels = np.arange(-80, 20, 5)\n", + "plt.contour(to_np(lons), to_np(lats), to_np(slp), 10, colors=\"black\",\n", + " transform=crs.PlateCarree())\n", + "plt.contourf(to_np(lons), to_np(lats), to_np(ctt), levels=levels, transform=crs.PlateCarree(),\n", + " cmap=get_cmap(\"Greys\"))\n", + "\n", + "# Add a color bar\n", + "plt.colorbar(ax=ax, shrink=.88)\n", + "\n", + "# Set the map limits. Not really necessary, but used for demonstration.\n", + "ax.set_xlim(cartopy_xlim(ctt))\n", + "ax.set_ylim(cartopy_ylim(ctt))\n", + "\n", + "# Add the gridlines\n", + "ax.gridlines(color=\"black\", linestyle=\"dotted\")\n", + "\n", + "plt.title(\"Cloud Top Temperature (degC)\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import numpy as np\n", + "from netCDF4 import Dataset\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.cm import get_cmap\n", + "import cartopy.crs as crs\n", + "from cartopy.feature import NaturalEarthFeature\n", + "\n", + "from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords\n", + "\n", + "# Open the NetCDF file\n", + "ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d01_2005-08-28_12:00:00\")\n", + "\n", + "# Get the sea level pressure\n", + "cfrac = getvar(ncfile, \"cfrac\")[2, :]\n", + "slp = getvar(ncfile, \"slp\")\n", + "\n", + "\n", + "# Get the latitude and longitude points\n", + "lats, lons = latlon_coords(ctt)\n", + "\n", + "# Get the cartopy mapping object\n", + "cart_proj = get_cartopy(ctt)\n", + "\n", + "# Create a figure\n", + "fig = plt.figure(figsize=(12,9))\n", + "# Set the GeoAxes to the projection used by WRF\n", + "ax = plt.axes(projection=cart_proj)\n", + "\n", + "# Download and add the states and coastlines\n", + "states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n", + " name='admin_1_states_provinces_shp')\n", + "ax.add_feature(states, linewidth=.5)\n", + "ax.coastlines('50m', linewidth=0.8)\n", + "\n", + "# Make the contour outlines and filled contours for the smoothed sea level pressure.\n", + "levels = np.arange(0, 1.1, .1)\n", + "plt.contour(to_np(lons), to_np(lats), to_np(slp), 10, colors=\"black\",\n", + " transform=crs.PlateCarree())\n", + "plt.contourf(to_np(lons), to_np(lats), to_np(cfrac), levels=levels, transform=crs.PlateCarree(),\n", + " cmap=get_cmap(\"Greys_r\"))\n", + "\n", + "# Add a color bar\n", + "plt.colorbar(ax=ax, shrink=.88)\n", + "\n", + "# Set the map limits. Not really necessary, but used for demonstration.\n", + "ax.set_xlim(cartopy_xlim(ctt))\n", + "ax.set_ylim(cartopy_ylim(ctt))\n", + "\n", + "# Add the gridlines\n", + "ax.gridlines(color=\"black\", linestyle=\"dotted\")\n", + "\n", + "plt.title(\"Cloud Fraction\")\n", + "\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.13" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/test/ipynb/Test_CTT_Prod.ipynb b/test/ipynb/Test_CTT_Prod.ipynb new file mode 100644 index 0000000..1bf8707 --- /dev/null +++ b/test/ipynb/Test_CTT_Prod.ipynb @@ -0,0 +1,214 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import numpy as np\n", + "from netCDF4 import Dataset\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.cm import get_cmap\n", + "import cartopy.crs as crs\n", + "from cartopy.feature import NaturalEarthFeature\n", + "\n", + "from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords\n", + "\n", + "# Open the NetCDF file\n", + "ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d01_2005-08-28_12:00:00\")\n", + "\n", + "# Get the sea level pressure\n", + "ctt = getvar(ncfile, \"ctt\")\n", + "slp = getvar(ncfile, \"slp\")\n", + "\n", + "\n", + "# Get the latitude and longitude points\n", + "lats, lons = latlon_coords(ctt)\n", + "\n", + "# Get the cartopy mapping object\n", + "cart_proj = get_cartopy(ctt)\n", + "\n", + "# Create a figure\n", + "fig = plt.figure(figsize=(12,9))\n", + "# Set the GeoAxes to the projection used by WRF\n", + "ax = plt.axes(projection=cart_proj)\n", + "\n", + "# Download and add the states and coastlines\n", + "states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n", + " name='admin_1_states_provinces_shp')\n", + "ax.add_feature(states, linewidth=.5)\n", + "ax.coastlines('50m', linewidth=0.8)\n", + "\n", + "# Make the contour outlines and filled contours for the smoothed sea level pressure.\n", + "levels = np.arange(-80, 20, 5)\n", + "plt.contour(to_np(lons), to_np(lats), to_np(slp), 10, colors=\"black\",\n", + " transform=crs.PlateCarree())\n", + "plt.contourf(to_np(lons), to_np(lats), to_np(ctt), levels=levels, transform=crs.PlateCarree(),\n", + " cmap=get_cmap(\"Greys\"))\n", + "\n", + "# Add a color bar\n", + "plt.colorbar(ax=ax, shrink=.88)\n", + "\n", + "# Set the map limits. Not really necessary, but used for demonstration.\n", + "ax.set_xlim(cartopy_xlim(ctt))\n", + "ax.set_ylim(cartopy_ylim(ctt))\n", + "\n", + "# Add the gridlines\n", + "ax.gridlines(color=\"black\", linestyle=\"dotted\")\n", + "\n", + "plt.title(\"Cloud Top Temperature (degC)\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import numpy as np\n", + "from netCDF4 import Dataset\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.cm import get_cmap\n", + "import cartopy.crs as crs\n", + "from cartopy.feature import NaturalEarthFeature\n", + "\n", + "from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords\n", + "\n", + "# Open the NetCDF file\n", + "ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d01_2005-08-28_12:00:00\")\n", + "\n", + "# Get the sea level pressure\n", + "ctt = getvar(ncfile, \"ctt\")\n", + "slp = getvar(ncfile, \"slp\")\n", + "\n", + "\n", + "# Get the latitude and longitude points\n", + "lats, lons = latlon_coords(ctt)\n", + "\n", + "# Get the cartopy mapping object\n", + "cart_proj = get_cartopy(ctt)\n", + "\n", + "# Create a figure\n", + "fig = plt.figure(figsize=(12,9))\n", + "# Set the GeoAxes to the projection used by WRF\n", + "ax = plt.axes(projection=cart_proj)\n", + "\n", + "# Download and add the states and coastlines\n", + "states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n", + " name='admin_1_states_provinces_shp')\n", + "ax.add_feature(states, linewidth=.5)\n", + "ax.coastlines('50m', linewidth=0.8)\n", + "\n", + "# Make the contour outlines and filled contours for the smoothed sea level pressure.\n", + "levels = np.arange(-80, 20, 5)\n", + "plt.contour(to_np(lons), to_np(lats), to_np(slp), 10, colors=\"black\",\n", + " transform=crs.PlateCarree())\n", + "plt.contourf(to_np(lons), to_np(lats), to_np(ctt), levels=levels, transform=crs.PlateCarree(),\n", + " cmap=get_cmap(\"Greys\"))\n", + "\n", + "# Add a color bar\n", + "plt.colorbar(ax=ax, shrink=.88)\n", + "\n", + "# Set the map limits. Not really necessary, but used for demonstration.\n", + "ax.set_xlim(cartopy_xlim(ctt))\n", + "ax.set_ylim(cartopy_ylim(ctt))\n", + "\n", + "# Add the gridlines\n", + "ax.gridlines(color=\"black\", linestyle=\"dotted\")\n", + "\n", + "plt.title(\"Cloud Top Temperature (degC)\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import numpy as np\n", + "from netCDF4 import Dataset\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.cm import get_cmap\n", + "import cartopy.crs as crs\n", + "from cartopy.feature import NaturalEarthFeature\n", + "\n", + "from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords\n", + "\n", + "# Open the NetCDF file\n", + "ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d01_2005-08-28_12:00:00\")\n", + "\n", + "# Get the sea level pressure\n", + "cfrac = getvar(ncfile, \"cfrac\")[2, :]\n", + "slp = getvar(ncfile, \"slp\")\n", + "\n", + "\n", + "# Get the latitude and longitude points\n", + "lats, lons = latlon_coords(ctt)\n", + "\n", + "# Get the cartopy mapping object\n", + "cart_proj = get_cartopy(ctt)\n", + "\n", + "# Create a figure\n", + "fig = plt.figure(figsize=(12,9))\n", + "# Set the GeoAxes to the projection used by WRF\n", + "ax = plt.axes(projection=cart_proj)\n", + "\n", + "# Download and add the states and coastlines\n", + "states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n", + " name='admin_1_states_provinces_shp')\n", + "ax.add_feature(states, linewidth=.5)\n", + "ax.coastlines('50m', linewidth=0.8)\n", + "\n", + "# Make the contour outlines and filled contours for the smoothed sea level pressure.\n", + "levels = np.arange(0, 1.1, .1)\n", + "plt.contour(to_np(lons), to_np(lats), to_np(slp), 10, colors=\"black\",\n", + " transform=crs.PlateCarree())\n", + "plt.contourf(to_np(lons), to_np(lats), to_np(cfrac), levels=levels, transform=crs.PlateCarree(),\n", + " cmap=get_cmap(\"Greys_r\"))\n", + "\n", + "# Add a color bar\n", + "plt.colorbar(ax=ax, shrink=.88)\n", + "\n", + "# Set the map limits. Not really necessary, but used for demonstration.\n", + "ax.set_xlim(cartopy_xlim(ctt))\n", + "ax.set_ylim(cartopy_ylim(ctt))\n", + "\n", + "# Add the gridlines\n", + "ax.gridlines(color=\"black\", linestyle=\"dotted\")\n", + "\n", + "plt.title(\"Cloud Fraction\")\n", + "\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From 78de49faeab9cdf8f46f9aa6dcb4c57b02b81c5c Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Mon, 19 Mar 2018 11:10:49 -0600 Subject: [PATCH 6/9] Updated CI test file for new CTT changes. Updated docs for raw ctt routine. --- src/wrf/computation.py | 16 ++++++++++------ test/ci_tests/ci_result_file.nc | Bin 2148317 -> 2148537 bytes 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/wrf/computation.py b/src/wrf/computation.py index 164d905..641b164 100644 --- a/src/wrf/computation.py +++ b/src/wrf/computation.py @@ -1111,12 +1111,16 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, True. Default is :data:`wrf.default_fill(numpy.float64)`. - opt_thresh (:obj:`float`, optional): The required amount of optical - depth (looking from top down) required to trigger a cloud top - temperature calculation. Values less than this threshold will be - treated as no cloud areas. For thin cirrus, a value of .1 might be - appropriate. For large CB clouds, 1000.0 might be appropriate. - Default is 1.0. + opt_thresh (:obj:`float`, optional): The amount of optical + depth (integrated from top down) required to trigger a cloud top + temperature calculation. The cloud top temperature is calculated at + the vertical level where this threshold is met. Vertical columns + with less than this threshold will be treated as cloud free areas. + In general, the larger the value is for this + threshold, the lower the altitude will be for the cloud top + temperature calculation, and therefore higher cloud top + temperature values. Default is 1.0, which should be sufficient for + most users. meta (:obj:`bool`): Set to False to disable metadata and return :class:`numpy.ndarray` instead of diff --git a/test/ci_tests/ci_result_file.nc b/test/ci_tests/ci_result_file.nc index 2cf9ffb9fafc9c006f3fcf9369eb985a5f187729..e15d100a4b5f5d5a4ec7e9526024de31d652a014 100644 GIT binary patch delta 11726 zcmb_i2{cvT_kSY~&ofWYT&55Tk&J~5nVL!FF*BJ`9wMPa$|aPckO+yFM21A>B4h}e zDl-kr^ncHLr(fUiTL1rl{nz^6wYc}Q_de&IyZ1hO?|04_j=v1AOuP(FkU&&6mrDA> z5iUhZS?~@NJTA;P$qz9ApG`7G?nZe$i9C};K=Sa#^yzwUoK5$gZpZwYs=+Ondnd%I5AWP zG1!3Y&Ot+3bZB`fei`9}Qu2g$CE0L$6J1vjh8C!*ry`ZaXGTa~ID(c2dQVa^;g&U$ zG8_~mwQr1zqr0oU>-pdaTFPkPD9X!WFr>Tf`41RCWui7fx_f!f(gP|ppn7>Yo3{We zMIA60+5+kIfXeI5P${Cyp-bV3Ox$TRRK#90$08V5N~2rxq>^yn;P=NgLq zSv*PSq&%RkNdbLm7>mY0lGaFel-%C65RT;1ir+?yOw1+5e0Q3(84VNOCwtGgtrllo+@_y(#1DXRBMJw;MJsKp7E$)wn zk3y4J(s1hdgJ2G_TSYSE4R+pj4Uiv5o&9Wb*ioe zR~PrKtqP$nA2?ij?P5HMxiThaHN_EoT=ZH~5Sq*JwkdneYi|P)QOhh0;mSCB%&43} zrrcf5QT|pXv$ZrBWKeh3OHV$1 zG9f?pw*mT-vry1&7C8Z%^m;{!Yub_pLanF=c4<}d2IJa^*)A5PgW5LtrOlogFci+p z*swCG>am*~eWv*qWfYTbVB2S|2hHHnf`lb8+nbu7Qlv+f#acv+ckzjjBZ_0O#&0E~ zcS$zxjW|-uhg9fFZ)@i6y521!B4zVbGch|`GFv23QrK6d(@YYbYD15b9=+L8fGUcW zHZC<~yfGzr!_nMa=NzBhV}l)rnER5ZJH+)hm0O#Y7h~-uMf161lJ*|@q-;`<-qK8< z?;shipjKFlWsxI^Cwv;-s%hq`&*7>QKE>>!FQ#nn6h?E6MTgG5jwd&Py(d18BU$wH zK?esGLz$Dox(^&eP{!JEQ?+(~RDGC@akhwg7V5?o z>=@ekW-N-ZCu)kD^>=KZ=ruEm&#Dn)Xv^b!Tn;DK3AreN+>6wyzbm4T84CV;N>^6BUx-drdsld(>?mANiZJ9UH$k1J} z2_kwir}kJfM0_$odOJAiu)fHv-LVSodf5ldp1~YBf(`8Cu=ZK5=D2KIWyH7$r;yA1 z=^ID`&9t2^0!_7->syt@;Ec~|V+jgTx>|MK*BW_r9OUT4Zn6hkqI4`wgZc$g9S=1X zCw7jC7)$WoVyXxep`$mfj|dY#qU%8;EFo7+e@)SlMO7chDrWhbo3kN{Q(D(P7TZ=g zdz1bZjkAbpr#_x7e?5iXJ&aSQC4V!aq99vGwB=UBBW0asQ|TLtE_aXVggX#i!-Y{S z7N$`=lXQsw=uqLH2o^nEvCK%hysln8F^Q-3i>&SWg|&6Cx*mE~-57*|gj@$JpV3u% zxSL*>{J~FhGy>cwl7sB#@{ghec=Qi?=$&K_!qVp2XuHV8Yk8KO$(*#2HNDv^)%|66 zR51#nd~{ceW5h&|#2)S`(<*NQqV-)#f{^Q-i1MAazL#G|J`^z*pfSnywsN%P6Arzb zWFEu%w!+4nB`4e+J+7AslM(W|k|fH_WpjMY-Zoo=Zb#jTlfpK8b*@BLYe|$anwMPl z8@=hEkGNc}=~Ev@rxI^_S>oXF;5(3>4W9d!jzv@kT%L9fei9)nDIHr$^1p_z#Skds z+&zq~zAjo>Kh6|(ns(HZ({{h6I4cT`q#+)19lhHeprFq3l08G<-Z#$a`JYTv!( zn&nWLy{I|TpHn8w_+`pinsAJp+Xd09d<|QAEi457qfmPQ0nb&jRl8WyzGi4bCx8=`_V`6WVsre+-X?zR{)3VNoPJZ5`9_y;=&>)6}2 z2IFMgcgN3NIP_UOdxvFcx9^Kj*(&!^z39yx!{|3q5}nG=U^aQRUG-Bri{T?RkJ%?~ z?hBbwdR*sknVK*VpC>klD|daj__GSj8D>!t_LE}L9BWHcPZYab z4V|XZUj?_Ol*QWv*01s$&W*lsZ%BU7a;@pxoUvQCVIe(VoFK3Jor(&-=fCHU%6wio zJ?j!MT^?=j?x5XM_A03_uhEKtNm0mcyl`MJ1G&?>O7QuYbI+3Qu<}+o9zv-U;pJw=k~>34R%9pG}YzQ>Tl~ zf0i`VIvug_bIlL_=&3pXw8i&WJu1bLIKQ`0zA9$O6z@^Ufda$o+M zYoZ(+bbtRfUpM|To23^M#^D?uDuiA`e}(g;ms&iV^Y=7#`i|Qg^|!q_HnyR0wg1

cTa_oN3;*Un3;4Xwh{>)*33b=I_PeV685ejT0Ud;c)QFnZwBnaT1e{5GKniu&#j1R4f@BU}oZ zvAVA0QBe`F%M8c(d&Bz2aMiPuFV6MYPrMj;T;uWXt)X*T@8GYqN1}?mJIat&SAT`- z4qj2S_B~#CVssrTTh^!SB_{P@7fiyl%gBe0#byaisiiY z`tF&$<83pF7pI7CS>_hJFL~|ru)5H+`d-Oy?G?Yij8AnSo_~GckeUDPh6A5kE3#5Xy6)@WE}5mmA-6)0-cUm#;gtG!^=vDpuEPi;;?O^W^=za=hSJ#*@4-=h%C&kX_YEpbywPQm?8?Y*u%7j|7% zbT4-~V7I5N#%y^gyF-{s#purcweglKUy?2prX!^_3WFBfJ!S`@@dpUwH7h>7ieIaT zzj=B#TCHW3^`!@txARnn7JEDW3~Ku1=4{fGt1dUN>K?Lq{bx#L;+M%{&$Gt6`3L-3 zMi=&0yIeH>)_&yY&j)eIs05jpg4rdX4l~VjReOz>|D1TbTP;HTJKxH+&sTqrypyn* zsaYl9(hTF8?nxC5O7;4Ayy`ec+f;2a`5`~CZlB+x#fVl&jkRBudh+d$dkkD#@#3nO z`O2T;ZLvN@vVnM&{Ja$(Zi!2XlzS=n#oO5E?vLSOblV(*gOTH43Wh5P$`a_(`a zw|{f$pAKhvS>ccqg>>5|+RI%%H>FbJ-Z9fo%bKiaqn^1>~$Z5f&SQ=z$6Z7h&iC3e&FC(1- z`~(tr5Go6-WPaaG>@*O}sI>THA8f8Yx%5b0Hkr5}pW1n`I)7mz_>AxC@5z0*n6L6p zP9_5jo9d^&p0-{Y>zKE{v(_Paar)@%cXF4$r622_D<1yv^84t0iLEH&X7i7q!{-rF z>)e&+l&0&{&z-FFcAJoU#B^A$+^;Lc)@kAWlF9Q$0v9p)%Y53o+@e{Z!Jj45xB7ac z)1y7g4`__`wS^@v?MY;uKi~gzrSOzUs6VC;VuK9X*N~-kR zP8a!$<8){UdBo|8J}=n;PxJV--E>OV?p8bxkSe}CbEc|BN~5d$@o>+S+mphO?|nKG zZ``V^@e0>Eyy-vmdyK^$@qF#6BBB`jeo5y|QQvoqjN!4llodGJ#{7Z6nUZ;-?jKKQ ze5yQ7NuHr?*)1NqnlW%ZarN!o!69pZ28NE zu+-OJwq5O4Uivcd4OCb9=M*uhJ|W8by1wZ9u2?x)b(V#uPFr~-ebYVZ(qd0~PP^xa z=X}9s0|UnCRp+jBc5U8$UMTNe-9Qtl=P|OW8~M2Wy4_xnmmJ6a`flQ;IC*Rvw#wy$ z4qB!(H8;^-#1r_^HmmF>d%^;D%cRn%t7L?ZFGbHMRf|tlEYf{mJ9BtOU_t&$IU~c% z**?Q1-r24#i>kQ*QNh-Z)i(`v>x_qtr}Bay&3DV_0V`I{e=HzlgRcno4zN1eU9P6O-W6^{Jyns z{8f3@b?xy8n#y23Uq4@84YfBy`v-_amfu7VNXf5!=@vJ+tBQWp+DiL2^Ruvq$z`Dz zlWOleb@83jc3S3y`8QxG~d13D_`D+arvdm#K~x=O}LVhaY;>4@%w;y zgZ`wZ5%!Su+T z)~6;xJFmj^!{LsJk92I_b|u`Ybu(1PTQ=}AEI+!btMKN-)klPM6Hi~i#6V512p7EW z?(YSKiks!}Uw^{{VG-7AYORkr!Cv0_U&yF(Mq z0?%66TC~?eLT<-kvw|=7Wiv)u5g$HeY2Hn9cO{`{{o-{z@}kac1MddMEiB4qUz#HE z4%hCTBEPP(?e^3ou={u^n3w3qCZx1S-HAPUAyT%h#`UCTCuUTZqfNIc7Vn?@G~sSC zV+27AFHi}qve<8?DIrkQ7Ip7P+N1F&tD5iMw&XNkR@P?EHn(e*L?3qNx-~~^yWXKV!t*LcO)#n97{at-(t~uB7=b2 zioB8%OK(<@?710nham51)^#MuQ`dWL@*X|@MT<5^4fpNVn^WS)IJrL4vP57{9Hz(L zxzv5nGwu;bn^=yhajV=3w@e3kyyc1DbvCQ(#qZv-cggIz!;I4^*U+oc+CnC zwdj)1lBwb7ZfjmF;kxVY9fq3ut(&oZO^uCBTv<2z&bd6yjXHEp_=&YCe~p-^5dMXM zkn>!3HlyC>)q~pW){7?=O~QGv(^>SnjAx~8VZFW-6lCr`{#(Xs>ZjW!g_FVWX9@K? zpIzz2<2Nh(-3Dun=9+7dedBfsQU2{#ZV;w#$ez3@D>?tQZ_tyVnfO5=Wcca~aUl|~ z(abgO&2M)6=-$Th=e^&# za2!fyg26LuiujKR*}aW9@{I8p5AOBdagM25eJ%e)27h|148L-5qkVXq5nYU!hr7a+ z@B+J@9Xe8yrK{h*I1)U(X$L>$alYtPoO-nyS)d-SYIqNF`F6BC4eOh#L`tfP#TR90SAK=h~?|X8k6A zr9~&0X?a$-X)#^ktnzKFmsF+yam4Puil;Ph;@a$@9oC#F)(Y&n{ zq<6$jXurJOlht2wTo|iU(Y+R-=i>IX(=aacE37s(8Zr;Z<42Z2Q=QhC_C~)Bc_%woQA5x#npU?C=6CE+1f`oGT|U zx#4FSD}V=DaMd(;90XmqokhTucm>P1w?TSR<*$fS*_ymlgfOusWT$C3eCa zwfSj6gU&U4#-G4c+!vb)oA|t`EAdf0bO58hz5SGcXpYUJ>g zh2!ic#P$Q6x{eH%l_^vh?9W<0mH%fIpNjvi-&66769XUQeH zJ>ERqbVTU_RM(1rjIzz;TBR&lFYnu4nBI>-c0#2?7&i5@rr1qNe_+v{+B@t@CuI`= z;GWih;wLVifu!A%lMqt}I-338t<7cErzs2mKRVEX2-i`_s}p?); zDUGX(20|fYM`(y<`*2%y#FCFRJ$y(T8rnP(rYA!G~7-*~S1|&tOM5MoUK#nV9+K=Jp0h=pulo*wZ_7CuIJih^j_hZ|JpH`|A1jf`8(F9QvCQmEz~)Hd8{GiV&qQk4SXOK!WPph zS^(IN18iqeD6Vb1hbdFw0TuPQ|(G;rMMcNuK<1h_Vfpa{TXYsP!}~J1W)4 zO`Z}o@XBnIfjehlMNs1@b|c%{Q(GBW3FK2}B1On{j)@|4=^jZ)Fex`;dr~7!7|~Fo zHS#=i2OZ0HoI@@nbciLH*jOF}XCqlPn{-5#SJ)&@#KY<6&H*g3bAXX>1|g|MG~I{0 z&yVP`B0nJLiK}OE{s@G=9~u6=(Px>?O?h*z8%~7h)F6-eKc_4T24k@({se8f;kfX9 zZj{*wNA03&18wi*YipmL56QUWgb-p)kf}S46Jgv8xwzv5_rG~V-QY#-rpjXi9DUEb z+S|G>wE&h3l|@Hp+539=`Ke7pweC1R_!Knkj>97kk3zrPal8btG0G%^qxMppU|l^) zbK2Y5-7&z?T~I}6n+!+EQ`xiu6#n4@fTKX=AcL%fGt~f9kxE4b1@kBY;yx;ombB4p zP1@E)a1j7Y3CO#`aL@v?Mw? zIRQ-eroS;A34gcno&<4&1j*O(8<#4oOGR=}Z1&0 zZa%m`Lj=qH>M&DvmKaI@##YCNj zO{+}Ard8&bOh`L^mVX6)M1uij!ol(w6_`+T4T5vB2rv$FQ$m4&q)C$=QJRzG9ty8E z3F0E<8=UvH$G~7u<+8OYGpEs#; zw}_Zx=(iI`vK+7gS^IxmfE)|{#^iGW3ytTeDZAR4h1_(|uro)I8t_QU$H8Dz9cIqk z&mfY6paOz8zz-OPqoL1V5J!;+*d)Z~%F&6Uabvb48;%n~pX1nvbi0e)p80fR78I{R zKt7pRRg}+nY#wp#9LGDjNP0fVY?SwxzB{Wzej<;Yd||KR7%5>S+Y3htBd}-u6n99n zIi*Z?8X}t~hZ_gdB7uXfI5?L0>jhp$Bf6x6gILnR8KaEkVCTvsaj?;`LbH9uCclVV z;6F!fsM4So2aQ234lt<20S1}{?**C#7-$w?pjj#=z(6$sIt0}Kz(Dn&15gbBKBxu& z26ZIBz*L|;uoMLoU@%}RAOK7SFetD91_c(tz~{jWf$0Drm=0iIx@}B)KMXiM5CEnG zlK_|!@IjOT7^F9Vr>6cH;{POD>IgqTZRFEUq$|iOy1tm3SYOO@xd0&?NIYzT|A9sH zxj`DvcuqP{TVTnD&meziJo6tCZkzPL881PuFgnT7Voj7J%eu`%yO5*-iG7=-=z^Dk zFF{r=crHAsH9)weRM{}tn>SFr3*HJrDlisY@SJ#1L4Xeea)K(5petSmLAsl`;&t($ zasUN6kV6oN#U#EC@@+MqIB$a7U?hgdbL`cK7GKAc{&wIceNX<5WYW#F);Kwtn zjZ@yobM5QLkHokJ{u~4!$eTdP)o&%YiyFj~W(8#Mn*TWBl_!+>n!f`j&2++nTorIa zSDW}t5&qoJfoA@QKM|PxLyO!Ah^yqESbt;kv<+SyzYond^B+fQzUGIrUh_B8DJYZZ zCq{O?gg9DA^kF;;g{~+^qIaxyg+O&CNO36a7LGRBt(8B|tyKVnIQ*ID(JgQki%55Y z>P7@QP_I_I2ce`@9mh<-{y72#R{=}igvBflOOBG)(MvwH5^UeaXye28@7j_5Wm;3GuJl3U&}JNSy7f zJPyZCE{u*Nk0^t;ldJuXSkY8+RcJ<3Qx~G+MDDN_qAXtMTM0F>!)J*senN1%KbaZW z@Be!(2;?8oKgd5+O!%{!0|Y<<0vO~UplpzT00#L7$bto8UZn17AireCkPQ|h7lVK#M*eFMsb8~ XEisxXq=rFcHs_HZFwpBJA+G-dFRTi{ delta 11922 zcmbt)2{@GB`?uK`Ok=IFgq*9i$hpZuF-<4fi`pOitR)l)8BugP{3Zv{>wk%PW zgruTuAyJZd#?0y0?|Z%P^?%>%eXfgVKKFf}^PJ~i&V4`UoToR#VPo;(umo}VSXc;RLbWmRhzgAKdk`;R7hDRrjz|X^2lk^e$kjG0jg_4K} zZobk%G#%xdIfvc8Z_zFfx7M{MD0{sFIH#UFJ8!Q2+jGVl}DEYSAjy&+E8zE$AoVs ziAX8)&G}GYUy{ZK_x8Q-&BH^bqC%x$q);jNXftXpo7@F7>57+y5|-i4$fkeJrX|@C zd`Rvqa9vB}@-c@~q$)iGFAQFRf?kt#O;oi-P)CK_LF*mu=Ir6_=zcy}hJiXB7?S#U zIy!`h{rM;~pt8^!AUwRiQkeji6;QoBT`zS2Dpf0VbTA*Tc0gk%OK$gp_VIJ}_rKuh zT+jn07$DBT>_~G42ox4-<7F#DUD`nY^-e}mh%Ew#l+J&B9=ei=ISIY9MFb*-s!>p* zA669$n@TgL4j})S*I-nE9fDL~hg3>~v*T~72p}qUprA#6Aw0>hihqXjrX32(Ou;n4 z+k~M$dQ=JWg(toeYUE`jp#omfgNOQ{QO?NtH3f6Xpc_R(o!p>t_CwDZQM|}&P3LEl z3?*@poxBhu1ZP4OD=oTHeUPRJx}7fnjVpcnc(-{#`9RS9wn`86DPQOrsL$UIjWCm+ zvZZ?a;r7Vj2koFT3@R1rX}cqcVp^zv2TGS{8e_{6ydKEm(sk21AgbKFs6)@m_-VX$ znKt1=UF6v^4Ydam7*Xr9gKS5%ZgeXdm6_}1i@y3+TsKiSQBtMe$4#{30$L3 zx~4E2pC%*j`2ajblZky#bXdzWUk~ej$z?fF^oYCco(>T`CS>mlPGx_2CR2KZKk~{H zJW7|<($4)R=YE^}IKN>Q#Ij8S)Qd&oEsn$B9e0>wWTXPj1~>TF^F(jsMG7AIIPLlB zdby?(O<_ct1ZsR$TkOPBdSKf3l8>=p9wN`%z2O-k&XQFDcTDOAPzHRNnL$sT+Mw(&4~p zzG6z z)RC!szsL)`fi5A(P$n4J%v_c4Aag@99Ok~R&${P+M#}+0RHKcMP0h*UYe!TfkKb2q zz&S8UYfHY;cv2J{li2x^uxly#1$&K9^Gwq)tKqTmT&^|^m0NDGqWxoW=?E7{uKsk^ti}C*I8M+bd4EQ6`@gjoIC(nJ63|_wEEU6U$n7!x8-b$cxDI7RIaL z-y87>A}aRU@~mH(l3w0{p)__9P=l!t&>13hjJ#5gn7GM4X3m`b8uTCS0<}?9je#O3 zdbLy8^EdA>oxRMZ(V4$xSXq#*Cfb=CRw%Er+$VYCrd!4-jVnYa_bWS*Y?_8w`9Ctk zUquD)3<_gAp}9Nr>Z3fFLB8GMk6xY^$+ zsJ@q>w~;B3T&lS~Sp(;~DUx+qM1n~%GTiu)uv|o2=Ao*6GJZZ>;@ThZxhikHw64%y zH(?4$5jrAZiMet_jwqGh@K}ydL_9V3-5lEIZHv>HVV%T}_DfDo)25t*Z0jrOFLaWl zOHko)gdGkN1g%)48|=b?bjikNm$R7!k-VbX4an0LS1)(txCIvIvlmR8%bPo0|PZD=@6V^x>j=C5i9(jNe^wx-@4wB*E02 zc#X*?kuwYCgO&dnE7|0GSz83n#EHqk+SIB}y^AnDxxkoYdZ=WtVnkT+`*&(ymIuXO~VcYD4x z)HJATCT5jgYwrwEQyAtQRI)Sn8qu(x-Q&Ng&QtesPewNK4p9`Z>+P$1OC+e@S2Es+ zuj}ep!bZQhZ(fd^RK~US&SHc1_Qt(d43Xht(u`3lj;8Y8*sX+NElvAVABTv(zu#W= zu#(Gc*6J7YXgfR|1tUqDd{x*PLD#D~!@pxnk}U~28Ji~5A>U@F23-_jWKwWi+tE9f zN#E`Gq9uoEX9Mwoxv?|$kai3Af*t3#EAL&CosfbGC%%d0{# zSG*weRzm)T%yLx~y61?vs6&CO>H~~kNZ{wf+wlhMv(x3ft3Hn##9zO73n9jFBX-qe zpsP>Rd+%w^f9@X}cAoNYcBq_w9Io4x9xqlRckss^vF8D&^qn`2(_$HZKKE2OGvacI zw_Q(Z>J!0i*mZyi<4ohMSbX_ls$N=W!*ls9_k?HJ<7h#SjT+nKeF~SI1IR<0sDaosHIVH=lT(lD9av;qH9) zS;Ebyp;krPgCoyuJVjmMEBXcMJUa6 z&6OwjP*;louLikOMETdGm8^rXQ0`hKGWYebXV*iIKkNQ(n|*Kgs4VZ46_*z_erPLr zao5eti|Rdv=LDX=8W7_Rice1ss2}YK3Am{9wqZPQ^lE(?>+h@~&FOn_Z`Q{rjy?K) zn!})%RIBWVACvS5Vc7X_tZ#`~Nta%96jS?E%O>j5E<~(??k3T!cQ&BmXEs9V>4eE_ z>%LD?LGped4m0FP{p|Iqa$m@be|p3I32ve(R%O?2oYoD2^aCGzy9)KIbEB@yua*zj zH-GNdz9%*Ewx<96`u6|>&mmTRmQ|wf)8oI&UHmi#9yJqRD-;agUx=&!EHE`Hxl?h` z1)D)($sLyQKe0|QbAy$kJ@X*)KJOJ4QZ~2m1EVb9*M7s537<|-3AIZ)Y#MB zK8)pE`VwEgY77-xA(@oX<*Ks$9y` zJjSNGn;=>DnAmCL{pFvc@r`n!SC#W3&}5JP4XK<^%bPw=XDpKFO9vJsy91uOp6MG2 zbqgP{d0NRkr>j(by`0!(xn}A+3v0bp8^+{2r}Qbb_g&ZX+}bwm3ZHm)pg`Zw zDUa~A$lxphdiK~p}~^*&1mgQ3<)>X==HE0PZZLA+~1tJ9Nyup9sFfb zVx`2MSA)5i&Lz${(mR-V&-7{SbkBufK^JBdgPAXxdCF(Xi+#2n-I3!`GX8xw#sB)3 zgT#u>VsE{~iqXzK6NR{J_cMtRN!6~SCRPf!QUjER`i4c1uJKND*OeXd&e-(*aq(8# z=C2rQL4}u92A#!io8zWGO4CjTY)GUu?Om82mmsdyp6XxHxEapXQ$f796g_h_KDIo) zW@TtaDbU39g|hW)j}vdc>KQGGeeb=M z)^7K{xd(R_^+5krUw^W=+ds%C@OicI29?SzBFJ z9I_fxTin5qRWz~LK~jK{9>4wW%ypNlULkUl$nlAcm%QI?W~<&J{;F8qE!b)Pjq7ynT7Z7x%dRzR#SO_~-_PE&*;413 z9tHVG_H{WV_KstB#aVf?q}QJBx>86=llomTXU-58d3tX4`q>IUhr<&KV=EIC$`Ylc zPozcz%qJ=i|1_;x_94!k&(>J1!2yw+_jqkH{rrNNf96bZ zz#HVTwY5WQVN8JrdPSr?1?BIeOOh6bixm6{4hy{6RN-jru<&-U3BB037lEB1#<84q zhcYMF=f%^l*s5jn7`O(qqukh@J}}jcLT`T4;7vR4$28tQvM@gs^WAI=^2|qDec*V! z`iXVHQEhE-;e5t?fhU3(ynN@M4XrrSq24II)!vz!;P&TD8KtRrQ|87(*B%yCb#Lr* z8M^ni-EnEr!ic5yOH`&6-F z;O6VKrcoV(fT4E^)1iu>d6(=rFS=f6TSRhph<>?YL`73%nZ zJsLB~3w;>RyM^4?#r@px92;a8u$u-w+LxRk?8W53j4 zmg^~HnJVgadZWL)ynaZ{e7JOY=0O_KYgIU|=C-%h_ZntPYsIZEJ9aGm$?v>thOR9e&^_0_c`-eQPNdrZjd0(yIrzN?mQ?ls3$k@@*U5(p{4N|xkPO#bxZsu{hWcU zxyw?B0Y!}E=L*wvKhF`r{6qvbBy(Taj2MD&yLMEhM~n6kTvP8!HGdH8ox*lidY_8+ z!44H1PK(8chyB68Yx;hO$It`EBN{1)LniWIgi7=WYgDEkaE15PlQdbpv&@r^5oPG` z2&CTbq}yg)*`^UiNcZfeNY$I)>2tW_iuh7GOmd&P6>F^<&)?-l{ld`wE)DOyyII5 zLgI$}r@r>R=iSmjV>E3i%kx69yUS$yN-RoI$24upJBYt-P|N=~sl`rR9;s1}6(ZJe zWo`u=_qWWwdicz8PrloX=!{wRzR2A2=j|uadLze{VqU(icIbTBa!yGuaG)rw4AUxI zoKtEcYVMIZsOitOT5-BV5vJJ?o_Zp`eo6dUN@Pk^BWteS^-2NDb>gI+^FNr>xID--g`f zed}jD%s>3*(1;{XI=>0+zOpqH|m&1sV#pA){_jnq1l6@+8FyAM$`Thcwb!nJQ`{m$fWM~@<+ z(HUMn@;cQuTa92wzR%aMTv-!*?)9d;@7hpgH?I0XZh(GfPr(-*z46%h*0-3uh!q25MY-b z@&8_9_=;OZLnGKY9p0LfkcfY6kRc>*vg1v1J2&qzYh>z&4}z@G?{b%snNL=<@jnId zXI7C#&)!J7*w+#Jo@e%H=r$i(i2bIaYiJ{I#>UhLA3u{xGG%JGRp95C{ZP^C9~B=r zCv%^^_wHFQ|FH{VwI`%8HMF^Z$QnG`Ao2QR?pWE!7=?iLXAE3#p9S0%d=T!6u-xB( z5HxHFf$g1bg2Y5jZfkzgiF@qFRb8#wk~3@Z=~rvrFxQ5I?OfcYGU@MIFMc!BL$ZBP z*&-7Gs@{bXM=(gENNgMoTg;@7(x~0%5z@tfF9DO7OHead1Xtg-&?YlwNqD^-Rfp(k zM{lnvBUjZ78=ywT&Gu?AobwQ6%~#HXb$fALtcSXS?nUc&`_|$k zFF?g$fgNxQXyO2`cT$@JJTJbYOWvR`gJyeBQMiA%7Vv3dpf0rc_o4z3*C!aMY7z5l z@rL%-!tA3~$huqUqlb`)l5zUKR@X_f6ZC_8hzJzJ-$Xllh71@&n5=;RtizLjvoWy3 z;IaZytQ-s^Ru1&*N(AE=nqoGpMJS8(i0f}KROQGIp+%5Z2~VC;_D;myWD`za(mbKA zN~g(9dnK}go>@yV!jK#^0>UT28@x6Qi34Ijfu+OL{UI0+FeDfQ;=Ge=F0_9bjU(3I0Zi(uKOG(XyyM;sK;)zK zgF5f#eBS=SJwT`K7Ey%H3HnlZU|>j^jmg9STsn{x{4+hT7O+HUEW`!dfHpTsas*9~20I)u9e8(u!CM0?%*PqM9Xc>LWH^Ra;C|Zuzlj-n?TtqiA6AL_ zNcOgmV~j0uME5Ksb-IWuoMR*v&M|2s5cmzIJhCrhqP}e~k-lv(=OrQPGckWTFxW^! z9FN9QUtl2>BSI{E=%`K`5ePyscK~k$_>UvfCKO+MULb3a57v4G{Xrv45!&>a)QjNm44%Tc9GA9_0@bvqe_($`X(+JcN zrh8yTbSkj@g#Vrs$a{!SThIkCq&Q7O7!U72TTgcI|a?R1A%=s0d#<~ z-+6aOdk;KB^u+SQ?LR>2o>)E$rccz^h9PBuBb@_r^AjQj^M z#XJATjAX187_uA8Q)EcF@RRF{G)E1`Q!Sk7ixol2FF7_+=6^2es}~*(P5WY5k!lP} zqtslb<%bP~p`m&|tU5AR>cIobXjMtt{#Xu9#N}ixHJWX2VM(^PSX9jst?Mi_x6Jv$ zLdyBUTF-)rFhC<#*pG1J#Lg;AbN_@Uu5DAWq1z zQ$~TrL`m#tCrRw*kY+%Hs&PN0SY5W=rvC|9xS6ETv z1n!FOXDlwEngb4%4WbW09Bz0E8DJ8MIsP^YWwQAjQ%(x(sqVQ#F?**Ha6IY@($T#U zU?B$?`J#XG-*Xj6C1CCf^FE*rEp_s9LdP!)aY23VcrjQKwC;}YMgGV_h)^uZ2~D2G zA40_6#Xxp0@0`UWLgVo?^IIk0NmdD*PXiH0IJjsuE*uw$3rDyb zjc6z({KqA-@(E3h5(88-BPgcif>=KZctT$v614xQlya!zZxvFK>fe}h{vo8|n;c^t^^bm65=gw#VJ zbaWX+i*W3RoEM!3+d2$bbNYgiaaI_FVNR1=EIj3Wq9*y8Vbiy8Vbdt^lqP zK-zVN=R0O6C`OoRa@5(BjtqyQI7;X+s3c*J*ts)d@K{&N-%(Al>; z0v}1>F@THP?ugBuM_d-)Nv&eI{xu>%2U9WVV!j;DWu?TFBJPasEv zcI0AZ!<|P8PQFW=<&$Gk+J@t9q0@C#SB-a-H95^4y)4*sL296xc&^oPenNV|- zVjI6QV#X5-33c+fGwwe~X71nOAEneJoKP!4kN}0W^QSPTa8bX)Sg1+eu!H|UNr0kZ z$^*dEPWb@+8&hTuqB{Et>H{Y__?6*Kt&m_Re-HAaheij5Y=VICux$?kVtv7nLed4> z@Kb69v8dxO_({iK2r!`G=cY;8uLO)Shz3_YMYJ0^g?7w0V9O$m_Tj3`N=?J>$iMp+mRO0KOG0~)0`=HD8nIHa0#&{8YRs_T|K$} zkm{|`kf9&e1@hrzL!+=JMm{n~(wY(Ql7WxxiXssNNhE@h0}igF zOyb=q^neL&N`gi&3teOd!`MZB?5We~!Qeb#B@PtUWXATD7vd@Thbuv@qiJ;N&A2Q1 zF8LRjIe-=B9Nr61edB~1o)r=&3P~Y8DF{(LQixK|&n7MJr+TEdf>0U!ayZ1QC=~Xm z3K1fU{H>my|H8s#+$(oQlRUCc9GIrawC(FldJPx#>7ZN_B~R@>hpEqp>Xn61$dAfX zaujzoqJF;d9}_}x!ZmH61%Dw<>RC&Yu&q!h4#rBF_ZNaOQjVW~4G|LfbNm$88{|Fk z>LBj{4DuepAn$?pAnyST(jLGd?a?ri4&*)1A;^0GgS-bi0C^AaLEZxxJC$0{KuG`yuPkNg7-Ww#4){V&*U(}w^6 From 8523dea59ab9f3b8b72feaab3754c7d0bc0eedbe Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Tue, 20 Mar 2018 16:22:07 -0600 Subject: [PATCH 7/9] Added support added for calculating cape_2d with a single vertical column. Fixes #46. --- doc/source/tutorial.rst | 3 +- doc/source/tutorials/wrf_workshop_2017.rst | 4 +- doc/source/tutorials/wrf_workshop_2018.rst | 410 +++++++++++++++++++++ src/wrf/computation.py | 35 +- src/wrf/g_cape.py | 6 +- src/wrf/metadecorators.py | 28 +- src/wrf/specialdec.py | 45 +-- test/comp_utest.py | 65 ++++ 8 files changed, 544 insertions(+), 52 deletions(-) create mode 100644 doc/source/tutorials/wrf_workshop_2018.rst diff --git a/doc/source/tutorial.rst b/doc/source/tutorial.rst index 29b68b5..dc05e15 100644 --- a/doc/source/tutorial.rst +++ b/doc/source/tutorial.rst @@ -12,7 +12,7 @@ Upcoming Tutorials .. toctree:: :maxdepth: 1 - tutorials/tutorial_03_2018.rst + tutorials/wrf_workshop_2018.rst Past Tutorials @@ -22,4 +22,5 @@ Past Tutorials :maxdepth: 1 tutorials/wrf_workshop_2017.rst + tutorials/tutorial_03_2018.rst diff --git a/doc/source/tutorials/wrf_workshop_2017.rst b/doc/source/tutorials/wrf_workshop_2017.rst index 44ecbef..d01ed3b 100644 --- a/doc/source/tutorials/wrf_workshop_2017.rst +++ b/doc/source/tutorials/wrf_workshop_2017.rst @@ -1,5 +1,5 @@ -WRF Workshop 2017 -===================== +WRF Users' Workshop 2017 +========================== Welcome wrf-python tutorial attendees! diff --git a/doc/source/tutorials/wrf_workshop_2018.rst b/doc/source/tutorials/wrf_workshop_2018.rst new file mode 100644 index 0000000..9a8f8e8 --- /dev/null +++ b/doc/source/tutorials/wrf_workshop_2018.rst @@ -0,0 +1,410 @@ +WRF Users' Workshop 2018 +========================= + +Welcome WRF-Python tutorial attendees! + +The instructions below should be completed prior to arriving at the tutorial. +There will not be enough time to do this during the tutorial. + +Prerequisites +--------------- + +This tutorial assumes that you have basic knowledge of how to type commands +in to a command terminal using your preferred operating system. You +should know some basic directory commands like *cd*, *mkdir*, *cp*, *mv*. + +Regarding Python, to understand the examples in this tutorial, you +should have some experience with Python basics. Below is a list of some +Python concepts that you will see in the examples, but don't worry if you aren't +familiar with everything. + +- Opening a Python interpreter and entering commands. +- Importing packages via the import statement. +- Familiarity with some of the basic Python types: str, list, tuple, dict, bool, float, int, None. +- Creating a list, tuple, or dict with "[ ]", "( )", "{ }" syntax (e.g. my_list = [1,2,3,4,5]). +- Accessing dict/list/tuple items with the "x[ ]" syntax (e.g. my_list_item = my_list[0]). +- Slicing str/list/tuple with the ":" syntax (e.g. my_slice = my_list[1:3]). +- Using object methods and attributes with the "x.y" syntax (e.g. my_list.append(6)). +- Calling functions (e.g. result = some_function(x, y)) +- Familiarity with numpy would be helpful, as only a very brief introduction + is provided. +- Familiarity with matplotlib would be helpful, as only a very brief + introduction is provided. + +If you are completely new to Python, that shouldn't be a problem, since +most of the examples consist of basic container types and function calls. It +would be helpful to look at some introductory material before arriving at the +tutorial. If you've programmed before, picking up Python isn't too difficult. + +Here are some links: + +https://www.learnpython.org/ + +https://developers.google.com/edu/python/ + + +Step 1: Open a Command Terminal +-------------------------------- + +To begin, you will first need to know how to open a command line terminal for +your operating system. + +For Windows: + +.. code-block:: none + + WINDOWS + r + type cmd in the run window + +For Mac: + +.. code-block:: none + + Finder -> Applications -> Utilities -> Terminal + +For Linux: + +.. code-block:: none + + Try one of the following: + + CTRL + ALT + T + CTRL + ALT + F2 + + +Step 2: Download Miniconda +---------------------------- + +For this tutorial, you will need to download and install Miniconda. We are +going to use Python 2.7, but it will also work with Python 3.5+. + +Please use the appropriate link below to download Miniconda for your operating +system. + +.. note:: + + 64-bit OS recommended + +`Win64 `_ + +`Mac `_ + +`Linux `_ + +For more information, see: https://conda.io/miniconda.html + +.. note:: + + **What is Miniconda?** + + If you have used the Anaconda distribution for Python before, then you will be + familiar with Miniconda. The Anaconda Python distribution includes numerous + scientific packages out of the box, which can be difficult for users to build and + install. More importantly, Anaconda includes the conda package manager. + + The conda package manager is a utility (similar to yum or apt-get) that installs + packages from a repository of pre-compiled Python packages. These repositories + are called channels. Conda makes it easy for Python users to install and + uninstall packages, and also can be used to create isolated Python environments + (more on that later). + + Miniconda is a bare bones implementation of Anaconda and only includes the + conda package manager. Since we are going to use the conda-forge channel to + install our scientific packages, Miniconda avoids any complications between + packages provided by Anaconda and conda-forge. + + +Step 3: Install Miniconda +---------------------------- + +Windows: + + 1. Browse to the directory where you downloaded Miniconda2-latest-Windows-x86_64.exe. + + 2. Double click on Miniconda2-latest-Windows-x86_64.exe. + + 3. Follow the instructions. + + 4. Restart your command terminal. + +Mac and Linux: + + For Mac and Linux, the installer is a bash script. + + 1. Using a terminal, you need to execute the bash shell script that you downloaded by + doing:: + + bash /path/to/Miniconda2-latest-MacOSX-x86_64.sh [Mac] + + bash /path/to/Miniconda2-latest-Linux-x86_64.sh [Linux] + + 2. Follow the instructions. + + 3. At the end of the installation, it will ask if you want to add the + miniconda2 path to your bash environment. If you are unsure what to do, + you should say "yes". If you say "no", we're going to assume you know + what you are doing. + + If you said "yes", then once you restart your shell, the miniconda2 Python + will be found instead of the system Python when you type the "python" + command. If you want to undo this later, then you can edit + either ~/.bash_profile or ~/.bashrc (depending on OS used) and + comment out the line that looks similar to:: + + # added by Miniconda2 4.1.11 installer + export PATH="/path/to/miniconda2/bin:$PATH" + + 4. Restart your command terminal. + + 5. [Linux and Mac Users Only] Miniconda only works with bash. If bash is + not your default shell, then you need to activate the bash shell by typing + the following in to your command terminal:: + + bash + + 6. Verify that your system is using the correct Python interpreter by typing + the following in to your command terminal:: + + which python + + You should see the path to your miniconda installation. If not, see the + note below. + + .. note:: + + If you have already installed another Python distribution, like Enthought + Canopy, you will need to comment out any PATH entries for that distribution + in your .bashrc or .bash_profile. Otherwise, your shell environment may + pick to wrong Python installation. + + If bash is not your default shell type, and the PATH variable has been + set in .bash_profile by the miniconda installer, try executing + "bash -l" instead of the "bash" command in step 5. + + +Step 4: Set Up the Conda Environment +-------------------------------------- + +If you are new to the conda package manager, one of the nice features of conda +is that you can create isolated Python environments that prevent package +incompatibilities. This is similar to the *virtualenv* package that some +Python users may be familiar with. However, conda is not compatible with +virtualenv, so only use conda environments when working with conda. + +The name of our conda environment for this tutorial is: **tutorial_2018**. + +Follow the instructions below to create the tutorial_2018 environment. + + 1. Open a command terminal if you haven't done so. + + 2. [Linux and Mac Users Only] The conda package manager only works with bash, + so if bash is not your current shell, type:: + + bash + + 3. Add the conda-forge channel to your conda package manager. + + Type or copy the command below in to your command terminal. You should + run this command even if you have already done it in the past. + This will ensure that conda-forge is set as the highest priority channel. + + :: + + conda config --add channels conda-forge + + .. note:: + + Conda-forge is a community driven collection of packages that are + continually tested to ensure compatibility. We highly recommend using + conda-forge when working with conda. See https://conda-forge.github.io/ + for more details on this excellent project. + + 4. Create the conda environment for the tutorial. + + Type or copy this command in to your command terminal:: + + conda create -n tutorial_2018 python=2.7 matplotlib cartopy netcdf4 jupyter git ffmpeg wrf-python + + Type "y" when prompted. It will take several minutes to install everything. + + This command creates an isolated Python environment named *tutorial_2018*, and installs + the python interpreter, matplotlib, cartopy, netcdf4, jupyter, git, ffmpeg, and wrf-python + packages. + + .. note:: + + When the installation completes, your command terminal might post a message similar to: + + .. code-block:: none + + If this is your first install of dbus, automatically load on login with: + + mkdir -p ~/Library/LaunchAgents + cp /path/to/miniconda2/envs/tutorial_test/org.freedesktop.dbus-session.plist ~/Library/LaunchAgents/ + launchctl load -w ~/Library/LaunchAgents/org.freedesktop.dbus-session.plist + + This is indicating that the dbus package can be set up to automatically load on login. You + can either ignore this message or type in the commands as indicated on your command terminal. + The tutorial should work fine in either case. + + 5. Activate the conda environment. + + To activate the tutorial_2018 Python environment, type the following + in to the command terminal: + + For Linux and Mac (using bash):: + + source activate tutorial_2018 + + For Windows:: + + activate tutorial_2018 + + You should see (tutorial_2018) on your command prompt. + + To deactivate your conda environment, type the following in to the + command terminal: + + For Linux and Mac:: + + source deactivate + + For Windows:: + + deactivate tutorial_2018 + + +Step 5: Download the Student Workbook +--------------------------------------- + +The student workbook for the tutorial is available on GitHub. The tutorial_2018 +conda environment includes the git application needed to download the repository. + +These instructions download the tutorial in to your home directory. If you want +to place the tutorial in to another directory, we're going to assume you know +how to do this yourself. + +To download the student workbook, follow these instructions: + + 1. Activate the tutorial_2018 conda environment following the instructions + in the previous step (*source activate tutorial_2018* or + *activate tutorial_2018*). + + 2. Change your working directory to the home directory by typing the + following command in to the command terminal: + + For Linux and Mac:: + + cd ~ + + For Windows:: + + cd %HOMEPATH% + + 3. Download the git repository for the tutorial by typing the following + in to the command terminal:: + + git clone https://github.com/NCAR/wrf_python_tutorial.git + + 4. There may be additional changes to the tutorial after you have downloaded + it. To pull down the latest changes, type the following in to the + command terminal: + + For Linux and Mac:: + + source activate tutorial_2018 + + cd ~/wrf_python_tutorial/wrf_workshop_2018 + + git pull + + For Windows:: + + activate tutorial_2018 + + cd %HOMEPATH%\wrf_python_tutorial\wrf_workshop_2018 + + git pull + + .. note:: + + If you try the "git pull" command and it returns an error indicating + that you have made changes to the workbook, this is probably because + you ran the workbook and it contains the cell output. To fix this, + first do a checkout of the workbook, then do the pull. + + .. code-block:: none + + git checkout -- . + git pull + + +Step 6: Verify Your Environment +---------------------------------- + +Verifying that your environment is correct involves importing a few +packages and checking for errors (you may see some warnings for matplotlib +or xarray, but you can safely ignore these). + + 1. Activate the tutorial_2018 conda environment if it isn't already active + (see instructions above). + + 2. Open a python terminal by typing the following in to the command + terminal:: + + python + + 3. Now type the following in to the Python interpreter:: + + >>> import netCDF4 + >>> import matplotlib + >>> import xarray + >>> import wrf + + 4. You can exit the Python interpreter using **CTRL + D** + + +Step 7: Obtain WRF Output Files +---------------------------------- + +For this tutorial, we strongly recommend that you use your own WRF output files. +The tutorial includes an easy way to point to your own data files. The WRF +output files should all be from the same WRF run and use the same domain. +If your files are located on another system (e.g. yellowstone), then copy 2 or +3 of these files to your local computer prior to the tutorial. + +If you do not have any of your own WRF output files, then you can download the +instructor data files from a link that should have been provided to you in an +email prior to the tutorial. + +If you are using the link provided in the email for your data, you can follow +the instructions below to place your data in the default location for your +workbook. + + 1. The link in the email should take you to a location on an Amazon cloud + drive. + + 2. If you hover your mouse over the wrf_tutorial_data.zip file, you'll see + an empty check box appear next to the file name. Click this check + box. + + 3. At the bottom of the screen, you'll see a Download button next to a + cloud icon. Click this button to start the download. + + 4. The download was most likely placed in to your ~/Downloads folder + [%HOMEPATH%\\Downloads for Windows]. Using your preferred method of choice + for unzipping files, unzip this file in to your home directory. Your data + should now be in ~/wrf_tutorial_data + [%HOMEPATH%\\wrf_tutorial_data for Windows]. + + 5. Verify that you have three WRF output files in that directory. + + +Getting Help +---------------- + +If you experience problems during this installation, please send a question +to the :ref:`google-group` support mailing list. + + +We look forward to seeing you at the tutorial! diff --git a/src/wrf/computation.py b/src/wrf/computation.py index 641b164..f181285 100644 --- a/src/wrf/computation.py +++ b/src/wrf/computation.py @@ -714,10 +714,10 @@ def smooth2d(field, passes, meta=True): @set_cape_alg_metadata(is2d=True, copyarg="pres_hpa") def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow, missing=default_fill(np.float64), meta=True): - """Return the two-dimensional CAPE, CIN, LCL, and LFC. + """Return the two-dimensional MCAPE, MCIN, LCL, and LFC. This function calculates the maximum convective available potential - energy (CAPE), maximum convective inhibition (CIN), + energy (MCAPE), maximum convective inhibition (MCIN), lifted condensation level (LCL), and level of free convection (LFC). This function uses the RIP [Read/Interpolate/plot] code to calculate potential energy (CAPE) and convective inhibition @@ -725,18 +725,28 @@ def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow, in the column (i.e. something akin to Colman's MCAPE). CAPE is defined as the accumulated buoyant energy from the level of free convection (LFC) to the equilibrium level (EL). CIN is defined as the accumulated negative - buoyant energy from the parcel starting point to the LFC. The word 'parcel' - here refers to a 500 meter deep parcel, with actual temperature and - moisture averaged over that depth. + buoyant energy from the parcel starting point to the LFC. + + The cape_2d algorithm works by first finding the maximum theta-e height + level in the lowest 3000 m. A parcel with a depth of 500 m is then + calculated and centered over this maximum theta-e height level. The + parcel's moisture and temperature characteristics are calculated by + averaging over the depth of this 500 m parcel. This 'maximum' parcel + is then used to compute MCAPE, MCIN, LCL and LFC. The leftmost dimension of the returned array represents four different quantities: - - return_val[0,...] will contain CAPE [J kg-1] - - return_val[1,...] will contain CIN [J kg-1] + - return_val[0,...] will contain MCAPE [J kg-1] + - return_val[1,...] will contain MCIN [J kg-1] - return_val[2,...] will contain LCL [m] - return_val[3,...] will contain LFC [m] + This function also supports computing MCAPE along a single vertical + column. In this mode, the *pres_hpa*, *tkel*, *qv* and *height* arguments + must be one-dimensional vertical columns, and the *terrain* and + *psfc_hpa* arguments must be scalar values + (:obj:`float`, :class:`numpy.float32` or :class:`numpy.float64`). This is the raw computational algorithm and does not extract any variables from WRF output files. Use :meth:`wrf.getvar` to both extract and compute @@ -749,6 +759,9 @@ def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow, least three dimensions. The rightmost dimensions can be top_bottom x south_north x west_east or bottom_top x south_north x west_east. + When operating on only a single column of values, the vertical + column can be bottom_top or top_bottom. In this case, *terrain* + and *psfc_hpa* must be scalars. Note: @@ -774,12 +787,16 @@ def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow, terrain (:class:`xarray.DataArray` or :class:`numpy.ndarray`): Terrain height in [m]. This is at least a two-dimensional array with the same dimensionality as *pres_hpa*, excluding the vertical - (bottom_top/top_bottom) dimension. + (bottom_top/top_bottom) dimension. When operating on a single + vertical column, this argument must be a scalar (:obj:`float`, + :class:`numpy.float32`, or :class:`numpy.float64`). psfc_hpa (:class:`xarray.DataArray` or :class:`numpy.ndarray`): The surface pressure in [hPa]. This is at least a two-dimensional array with the same dimensionality as *pres_hpa*, excluding the - vertical (bottom_top/top_bottom) dimension. + vertical (bottom_top/top_bottom) dimension. When operating on a + singlevertical column, this argument must be a scalar + (:obj:`float`, :class:`numpy.float32`, or :class:`numpy.float64`). Note: diff --git a/src/wrf/g_cape.py b/src/wrf/g_cape.py index 5f5ef10..0a43c7a 100755 --- a/src/wrf/g_cape.py +++ b/src/wrf/g_cape.py @@ -13,13 +13,13 @@ from .metadecorators import set_cape_metadata @set_cape_metadata(is2d=True) def get_2dcape(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, meta=True, _key=None, missing=default_fill(np.float64)): - """Return the 2d fields of CAPE, CIN, LCL, and LFC. + """Return the 2d fields of MCAPE, MCIN, LCL, and LFC. The leftmost dimension of the returned array represents four different quantities: - - return_val[0,...] will contain CAPE [J kg-1] - - return_val[1,...] will contain CIN [J kg-1] + - return_val[0,...] will contain MCAPE [J kg-1] + - return_val[1,...] will contain MCIN [J kg-1] - return_val[2,...] will contain LCL [m] - return_val[3,...] will contain LFC [m] diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index b35d57d..40e8e49 100644 --- a/src/wrf/metadecorators.py +++ b/src/wrf/metadecorators.py @@ -1985,7 +1985,7 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): p = argvals[copyarg] missing = argvals["missing"] - # Note: cape_3d supports using only a single column of data + # Note: 2D/3D cape supports using only a single column of data is1d = p.ndim == 1 # Need to squeeze the right dimensions for 1D cape @@ -1997,31 +1997,34 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): outattrs = OrderedDict() + if is2d: - outname = "cape_2d" + if is1d: + outname = "cape_2d" + else: + outname = "cape_2d_column" outattrs["description"] = "mcape ; mcin ; lcl ; lfc" outattrs["units"] = "J kg-1 ; J kg-1 ; m ; m" outattrs["MemoryOrder"] = "XY" + outattrs["MemoryOrder"] = "" else: if not is1d: outname = "cape_3d" - outattrs["description"] = "cape; cin" - outattrs["units"] = "J kg-1 ; J kg-1" outattrs["MemoryOrder"] = "XYZ" else: - outname = "cape_column" - outattrs["description"] = "cape; cin" - outattrs["units"] = "J kg-1 ; J kg-1" + outname = "cape_3d_column" outattrs["MemoryOrder"] = "Z" + outattrs["description"] = "cape; cin" + outattrs["units"] = "J kg-1 ; J kg-1" if isinstance(p, DataArray): if is2d: - # Right dims - outdims[-2:] = p.dims[-2:] - # Left dims - outdims[1:-2] = p.dims[0:-3] - + if not is1d: + # Right dims + outdims[-2:] = p.dims[-2:] + # Left dims + outdims[1:-2] = p.dims[0:-3] else: if not is1d: # Right dims @@ -2030,6 +2033,7 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): outdims[1:-3] = p.dims[0:-3] else: outdims[1] = p.dims[0] + outcoords = {} # Left-most is always cape_cin or cape_cin_lcl_lfc diff --git a/src/wrf/specialdec.py b/src/wrf/specialdec.py index f72899b..1795f25 100644 --- a/src/wrf/specialdec.py +++ b/src/wrf/specialdec.py @@ -229,7 +229,14 @@ def cape_left_iter(alg_dtype=np.float64): ter_follow = args[8] is2d = i3dflag == 0 - is1d = np.isscalar(sfp) + # Note: This should still work with DataArrays + is1d = np.isscalar(sfp) or np.size(sfp) == 1 + + # Make sure sfp and terrain are regular floats for 1D case + # This should also work with DataArrays + if is1d: + ter = float(ter) + sfp = float(sfp) orig_dtype = p_hpa.dtype @@ -542,39 +549,27 @@ def check_cape_args(): ter_follow = args[8] is2d = False if i3dflag != 0 else True + is1d = ((np.isscalar(sfp) or np.size(sfp) == 1) or + (np.isscalar(ter) or np.size(ter) == 1)) if not (p_hpa.shape == tk.shape == qv.shape == ht.shape): raise ValueError("arguments 0, 1, 2, 3 must be the same shape") # 2D CAPE does not allow for scalars - if is2d: - if np.isscalar(ter) or np.isscalar(sfp): - raise ValueError("arguments 4 and 5 cannot be scalars with " - "cape_2d routine") - + if not is1d: if ter.ndim != p_hpa.ndim-1 or sfp.ndim != p_hpa.ndim-1: raise ValueError("arguments 4 and 5 must have " "{} dimensions".format(p_hpa.ndim-1)) - - # 3D cape is allowed to be just a vertical column with scalars - # for terrain and psfc_hpa. else: - if np.isscalar(ter) and np.isscalar(sfp): - if p_hpa.ndim != 1: - raise ValueError("arguments 0-3 " - "must be 1-dimensional when " - "arguments 4 and 5 are scalars") - if is2d: - raise ValueError("argument 7 must be 0 when using 1D data") - else: - if ((np.isscalar(ter) and not np.isscalar(sfp)) or - (not np.isscalar(ter) and np.isscalar(sfp))): - raise ValueError("arguments 4 and 5 must both be scalars") - else: - if ter.ndim != p_hpa.ndim-1 or sfp.ndim != p_hpa.ndim-1: - raise ValueError("arguments 4 and 5 " - "must have {} dimensions".format( - p_hpa.ndim-1)) + if np.size(ter) != np.size(sfp): + raise ValueError("arguments 4 and 5 must both be scalars or " + "both be arrays") + + # Only need to test p_hpa since we assured args 0-3 have same ndim + if p_hpa.ndim != 1: + raise ValueError("arguments 0-3 " + "must be 1-dimensional when " + "arguments 4 and 5 are scalars") return wrapped(*args, **kwargs) diff --git a/test/comp_utest.py b/test/comp_utest.py index 958cd75..864bb7a 100644 --- a/test/comp_utest.py +++ b/test/comp_utest.py @@ -597,9 +597,71 @@ def test_cape3d_1d(wrfnc): nt.assert_allclose(to_np(result), to_np(ref)) return func + + +def test_cape2d_1d(wrfnc): + + def func(self): + varnames = ("T", "P", "PB", "QVAPOR", "PH", "PHB", "HGT", "PSFC") + ncvars = extract_vars(wrfnc, 0, varnames, method="cat", squeeze=True, + cache=None, meta=True) + + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + qv = ncvars["QVAPOR"] + ph = ncvars["PH"] + phb = ncvars["PHB"] + ter = ncvars["HGT"] + psfc = ncvars["PSFC"] + + col_idxs = (slice(None), t.shape[-2]//2, t.shape[-1]//2) + t = t[col_idxs] + p = p[col_idxs] + pb = pb[col_idxs] + qv = qv[col_idxs] + ph = ph[col_idxs] + phb = phb[col_idxs] + ter = float(ter[col_idxs[1:]]) + psfc = float(psfc[col_idxs[1:]]) + + full_t = t + Constants.T_BASE + full_p = p + pb + tkel = tk(full_p, full_t, meta=False) + + geopt = ph + phb + geopt_unstag = destagger(geopt, -1) + z = geopt_unstag/Constants.G + + # Convert pressure to hPa + p_hpa = ConversionFactors.PA_TO_HPA * full_p + psfc_hpa = ConversionFactors.PA_TO_HPA * psfc + + i3dflag = 0 + ter_follow = 1 + + result = cape_2d(p_hpa, tkel, qv, z, ter, psfc_hpa, ter_follow) + + print ("RESULT", result) + + ref = getvar(wrfnc, "cape_2d") + + ref = ref[(slice(None),) + col_idxs[1:]] + + print ("REF", ref) + + nt.assert_allclose(to_np(result), to_np(ref)) + + return func if __name__ == "__main__": + from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule, + omp_set_dynamic, omp_get_num_procs, OMP_SCHED_STATIC) + omp_set_num_threads(omp_get_num_procs()-1) + omp_set_schedule(OMP_SCHED_STATIC, 0) + omp_set_dynamic(False) + varnames = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz", "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", @@ -639,6 +701,9 @@ if __name__ == "__main__": func = test_cape3d_1d(wrfnc) setattr(WRFVarsTest, "test_cape3d_1d", func) + func = test_cape2d_1d(wrfnc) + setattr(WRFVarsTest, "test_cape2d_1d", func) + ut.main() From 5799ccf6d0f84a53009d44a5e1c462738c2967a2 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Thu, 29 Mar 2018 15:12:18 -0600 Subject: [PATCH 8/9] Fixed bugs related to input failures with dictionaries. Modified routines that use qvapor to make copies so that scipy.io.netcdf works in the newer release. Added input tests. --- doc/source/new.rst | 2 + src/wrf/g_dewpoint.py | 8 +- src/wrf/g_latlon.py | 131 ++++++++++++++++++++-- src/wrf/g_rh.py | 8 +- src/wrf/g_slp.py | 4 +- src/wrf/interp.py | 22 ++-- src/wrf/latlonutils.py | 27 ++--- src/wrf/metadecorators.py | 8 +- src/wrf/util.py | 3 +- test/comp_utest.py | 6 +- test/ipynb/Doc_Examples.ipynb | 2 +- test/ipynb/Test_CTT.ipynb | 2 +- test/ipynb/WRF_python_demo.ipynb | 37 +++++-- test/test_inputs.py | 181 +++++++++++++++++++++++++++++++ test/utests.py | 2 +- 15 files changed, 381 insertions(+), 62 deletions(-) create mode 100644 test/test_inputs.py diff --git a/doc/source/new.rst b/doc/source/new.rst index c07e5a5..6103c43 100644 --- a/doc/source/new.rst +++ b/doc/source/new.rst @@ -21,8 +21,10 @@ v1.1.3 - Added 'th' alias for the theta product. - Fixed a crash issue related to updraft helicity when a dictionary is used as the input. +- Dictionary inputs now work correctly with xy_to_ll and ll_to_xy. - The cape_2d diagnostic can now work with a single column of data, just like cape_3d. + v1.1.2 ^^^^^^^^^^^^^^ diff --git a/src/wrf/g_dewpoint.py b/src/wrf/g_dewpoint.py index ad4a245..e9778d2 100755 --- a/src/wrf/g_dewpoint.py +++ b/src/wrf/g_dewpoint.py @@ -76,7 +76,9 @@ def get_dp(wrfin, timeidx=0, method="cat", squeeze=True, p = ncvars["P"] pb = ncvars["PB"] - qvapor = ncvars["QVAPOR"] + # Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to + # break with every release + qvapor = ncvars["QVAPOR"].copy() # Algorithm requires hPa full_p = .01*(p + pb) @@ -152,7 +154,9 @@ def get_dp_2m(wrfin, timeidx=0, method="cat", squeeze=True, # Algorithm requires hPa psfc = .01*(ncvars["PSFC"]) - q2 = ncvars["Q2"] + # Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to + # break with every release + q2 = ncvars["Q2"].copy() q2[q2 < 0] = 0 td = _td(psfc, q2) diff --git a/src/wrf/g_latlon.py b/src/wrf/g_latlon.py index d25f64a..44363c9 100755 --- a/src/wrf/g_latlon.py +++ b/src/wrf/g_latlon.py @@ -1,9 +1,16 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) -from .util import extract_vars, get_id -from .latlonutils import (_lat_varname, _lon_varname, _ll_to_xy, _xy_to_ll) +from collections import OrderedDict + +import numpy as np +from xarray import DataArray + +from .util import extract_vars, get_id, get_iterable, is_mapping, to_np +from .py3compat import viewkeys +from .latlonutils import _lat_varname, _lon_varname, _ll_to_xy, _xy_to_ll from .metadecorators import set_latlon_metadata +from .config import xarray_enabled def get_lat(wrfin, timeidx=0, method="cat", squeeze=True, @@ -151,7 +158,101 @@ def get_lon(wrfin, timeidx=0, method="cat", squeeze=True, return lon_var[varname] -# TODO: Do we need the user to know about method, squeeze, cache for this? + +def _llxy_mapping(wrfin, x_or_lat, y_or_lon, func, timeidx, stagger, + squeeze, meta, as_int=None): + + keynames = [] + # This might not work once mapping iterators are implemented + numkeys = len(wrfin) + + key_iter = iter(viewkeys(wrfin)) + first_key = next(key_iter) + keynames.append(first_key) + + first_args = [wrfin[first_key], x_or_lat, y_or_lon, timeidx, squeeze, + meta, stagger] + if as_int is not None: + first_args.append(as_int) + + first_array = func(*first_args) + + # Create the output data numpy array based on the first array + outdims = [numkeys] + outdims += first_array.shape + outdata = np.empty(outdims, first_array.dtype) + outdata[0,:] = first_array[:] + + idx = 1 + while True: + try: + key = next(key_iter) + except StopIteration: + break + else: + keynames.append(key) + + args = [wrfin[first_key], x_or_lat, y_or_lon, timeidx, squeeze, + meta, stagger] + if as_int is not None: + args.append(as_int) + + result_array = func(*args) + + if outdata.shape[1:] != result_array.shape: + raise ValueError("data sequences must have the " + "same size for all dictionary keys") + outdata[idx,:] = to_np(result_array)[:] + idx += 1 + + if xarray_enabled() and meta: + outname = str(first_array.name) + # Note: assumes that all entries in dict have same coords + outcoords = OrderedDict(first_array.coords) + + # First find and store all the existing key coord names/values + # This is applicable only if there are nested dictionaries. + key_coordnames = [] + coord_vals = [] + existing_cnt = 0 + while True: + key_coord_name = "key_{}".format(existing_cnt) + + if key_coord_name not in first_array.dims: + break + + key_coordnames.append(key_coord_name) + coord_vals.append(to_np(first_array.coords[key_coord_name])) + + existing_cnt += 1 + + # Now add the key coord name and values for THIS dictionary. + # Put the new key_n name at the bottom, but the new values will + # be at the top to be associated with key_0 (left most). This + # effectively shifts the existing 'key_n' coordinate values to the + # right one dimension so *this* dicionary's key coordinate values + # are at 'key_0'. + key_coordnames.append(key_coord_name) + coord_vals.insert(0, keynames) + + # make it so that key_0 is leftmost + outdims = key_coordnames + list(first_array.dims[existing_cnt:]) + + + # Create the new 'key_n', value pairs + for coordname, coordval in zip(key_coordnames, coord_vals): + outcoords[coordname] = coordval + + + outattrs = OrderedDict(first_array.attrs) + + outarr = DataArray(outdata, name=outname, coords=outcoords, + dims=outdims, attrs=outattrs) + else: + outarr = outdata + + return outarr + @set_latlon_metadata(xy=True) def ll_to_xy(wrfin, latitude, longitude, timeidx=0, @@ -215,9 +316,15 @@ def ll_to_xy(wrfin, latitude, longitude, timeidx=0, be a :class:`numpy.ndarray` object with no metadata. """ + if is_mapping(wrfin): + return _llxy_mapping(wrfin, latitude, longitude, ll_to_xy, + timeidx, stagger, squeeze, meta, as_int) _key = get_id(wrfin) - return _ll_to_xy(latitude, longitude, wrfin, timeidx, stagger, "cat", + _wrfin = get_iterable(wrfin) + + return _ll_to_xy(latitude, longitude, _wrfin, timeidx, stagger, "cat", squeeze, None, _key, as_int, **{}) + @set_latlon_metadata(xy=True) @@ -319,10 +426,10 @@ def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True, return _ll_to_xy(latitude, longitude, None, 0, True, "cat", squeeze, None, None, as_int, **projparams) - - + + @set_latlon_metadata(xy=False) -def xy_to_ll(wrfin, x, y, timeidx=0, stagger=None, squeeze=True, meta=True): +def xy_to_ll(wrfin, x, y, timeidx=0, squeeze=True, meta=True, stagger=None): """Return the latitude and longitude for specified x,y coordinates. The *x* and *y* arguments can be a single value or a sequence of values. @@ -370,9 +477,6 @@ def xy_to_ll(wrfin, x, y, timeidx=0, stagger=None, squeeze=True, meta=True): - 'v': Use the same staggered grid as the v wind component, which has a staggered south_north (y) dimension. - as_int (:obj:`bool`): Set to True to return the x,y values as - :obj:`int`, otherwise they will be returned as :obj:`float`. - Returns: :class:`xarray.DataArray` or :class:`numpy.ndarray`: The latitude and longitude values whose leftmost dimension is 2 @@ -382,8 +486,13 @@ def xy_to_ll(wrfin, x, y, timeidx=0, stagger=None, squeeze=True, meta=True): be a :class:`numpy.ndarray` object with no metadata. """ + if is_mapping(wrfin): + return _llxy_mapping(wrfin, x, y, xy_to_ll, + timeidx, stagger, squeeze, meta) + _key = get_id(wrfin) - return _xy_to_ll(x, y, wrfin, timeidx, stagger, "cat", True, None, + _wrfin = get_iterable(wrfin) + return _xy_to_ll(x, y, _wrfin, timeidx, stagger, "cat", True, None, _key, **{}) diff --git a/src/wrf/g_rh.py b/src/wrf/g_rh.py index 5af8ae5..6e4bb6b 100755 --- a/src/wrf/g_rh.py +++ b/src/wrf/g_rh.py @@ -71,7 +71,9 @@ def get_rh(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, t = ncvars["T"] p = ncvars["P"] pb = ncvars["PB"] - qvapor = ncvars["QVAPOR"] + # Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to + # break with every release + qvapor = ncvars["QVAPOR"].copy() full_t = t + Constants.T_BASE full_p = p + pb @@ -144,7 +146,9 @@ def get_rh_2m(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, meta=False, _key=_key) t2 = ncvars["T2"] psfc = ncvars["PSFC"] - q2 = ncvars["Q2"] + # Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to + # break with every release + q2 = ncvars["Q2"].copy() q2[q2 < 0] = 0 rh = _rh(q2, psfc, t2) diff --git a/src/wrf/g_slp.py b/src/wrf/g_slp.py index a1f5a96..77bac48 100755 --- a/src/wrf/g_slp.py +++ b/src/wrf/g_slp.py @@ -83,7 +83,9 @@ def get_slp(wrfin, timeidx=0, method="cat", squeeze=True, t = ncvars["T"] p = ncvars["P"] pb = ncvars["PB"] - qvapor = ncvars["QVAPOR"] + # Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to + # break with every release + qvapor = ncvars["QVAPOR"].copy() ph = ncvars["PH"] phb = ncvars["PHB"] diff --git a/src/wrf/interp.py b/src/wrf/interp.py index 1892947..34d7e6f 100755 --- a/src/wrf/interp.py +++ b/src/wrf/interp.py @@ -8,7 +8,7 @@ from .extension import (_interpz3d, _vertcross, _interpline, _smooth2d, _monotonic, _vintrp) from .metadecorators import set_interp_metadata -from .util import extract_vars, is_staggered, get_id, to_np +from .util import extract_vars, is_staggered, get_id, to_np, get_iterable from .py3compat import py3range from .interputils import get_xy, get_xy_z_params, to_xy_coords from .constants import Constants, default_fill, ConversionFactors @@ -548,6 +548,8 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, """ _key = get_id(wrfin) + _wrfin = get_iterable(wrfin) + # Remove case sensitivity field_type = field_type.lower() if field_type is not None else "none" vert_coord = vert_coord.lower() if vert_coord is not None else "none" @@ -583,7 +585,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, interp_levels = np.asarray(interp_levels, np.float64) # TODO: Check if field is staggered - if is_staggered(wrfin, field): + if is_staggered(_wrfin, field): raise RuntimeError("Please unstagger field in the vertical") # Check for valid coord @@ -606,23 +608,23 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, # Extract vriables #timeidx = -1 # Should this be an argument? - ncvars = extract_vars(wrfin, timeidx, ("PSFC", "QVAPOR", "F"), + ncvars = extract_vars(_wrfin, timeidx, ("PSFC", "QVAPOR", "F"), method, squeeze, cache, meta=False, _key=_key) sfp = ncvars["PSFC"] * ConversionFactors.PA_TO_HPA qv = ncvars["QVAPOR"] coriolis = ncvars["F"] - terht = get_terrain(wrfin, timeidx, units="m", + terht = get_terrain(_wrfin, timeidx, units="m", method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) - tk = get_temp(wrfin, timeidx, units="k", + tk = get_temp(_wrfin, timeidx, units="k", method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) - p = get_pressure(wrfin, timeidx, units="pa", + p = get_pressure(_wrfin, timeidx, units="pa", method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) - ght = get_height(wrfin, timeidx, msl=True, units="m", + ght = get_height(_wrfin, timeidx, msl=True, units="m", method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) @@ -639,7 +641,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, vcord_array = np.exp(-ght/sclht) elif vert_coord == "ght_agl": - ht_agl = get_height(wrfin, timeidx, msl=False, units="m", + ht_agl = get_height(_wrfin, timeidx, msl=False, units="m", method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) @@ -647,7 +649,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, vcord_array = np.exp(-ht_agl/sclht) elif vert_coord in ("theta", "th"): - t = get_theta(wrfin, timeidx, units="k", + t = get_theta(_wrfin, timeidx, units="k", method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) @@ -671,7 +673,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, idir = 1 delta = 0.01 - eth = get_eth(wrfin, timeidx, method=method, squeeze=squeeze, + eth = get_eth(_wrfin, timeidx, method=method, squeeze=squeeze, cache=cache, meta=False, _key=_key) p_hpa = p * ConversionFactors.PA_TO_HPA diff --git a/src/wrf/latlonutils.py b/src/wrf/latlonutils.py index ccd5221..00940ef 100644 --- a/src/wrf/latlonutils.py +++ b/src/wrf/latlonutils.py @@ -131,8 +131,9 @@ def _get_proj_params(wrfin, timeidx, stagger, method, squeeze, cache, _key): """ - if timeidx < 0: - raise ValueError("'timeidx' must be greater than 0") + if timeidx is not None: + if timeidx < 0: + raise ValueError("'timeidx' must be greater than 0") attrs = extract_global_attrs(wrfin, attrs=("MAP_PROJ", "TRUELAT1", "TRUELAT2", "STAND_LON", @@ -383,12 +384,10 @@ def _ll_to_xy(latitude, longitude, wrfin=None, timeidx=0, if ref_lat.size == 1: outdim = [2, lats.size] - #outdim = [lats.size, 2] extra_dims = [outdim[1]] else: # Moving domain will have moving ref_lats/ref_lons outdim = [2, ref_lat.size, lats.size] - #outdim = [lats.size, ref_lat.size, 2] extra_dims = outdim[1:] result = np.empty(outdim, np.float64) @@ -402,11 +401,11 @@ def _ll_to_xy(latitude, longitude, wrfin=None, timeidx=0, ref_lat_val = ref_lat[0] ref_lon_val = ref_lon[0] else: - ref_lat_val = ref_lat[left_idxs[-1]] - ref_lon_val = ref_lon[left_idxs[-1]] + ref_lat_val = ref_lat[left_idxs[-2]] + ref_lon_val = ref_lon[left_idxs[-2]] - lat = lats[left_idxs[0]] - lon = lons[left_idxs[0]] + lat = lats[left_idxs[-1]] + lon = lons[left_idxs[-1]] xy = _lltoxy(map_proj, truelat1, truelat2, stdlon, ref_lat_val, ref_lon_val, pole_lat, pole_lon, @@ -548,14 +547,10 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None, raise ValueError("'x' and 'y' must be the same length") if ref_lat.size == 1: - #outdim = [x_arr.size, 2] - #extra_dims = [outdim[0]] outdim = [2, x_arr.size] extra_dims = [outdim[1]] else: # Moving domain will have moving ref_lats/ref_lons - #outdim = [x_arr.size, ref_lat.size, 2] - #extra_dims = outdim[0:2] outdim = [2, ref_lat.size, x_arr.size] extra_dims = outdim[1:] @@ -570,11 +565,11 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None, ref_lat_val = ref_lat[0] ref_lon_val = ref_lon[0] else: - ref_lat_val = ref_lat[left_idxs[-1]] - ref_lon_val = ref_lon[left_idxs[-1]] + ref_lat_val = ref_lat[left_idxs[-2]] + ref_lon_val = ref_lon[left_idxs[-2]] - x_val = x_arr[left_idxs[0]] - y_val = y_arr[left_idxs[0]] + x_val = x_arr[left_idxs[-1]] + y_val = y_arr[left_idxs[-1]] ll = _xytoll(map_proj, truelat1, truelat2, stdlon, ref_lat_val, ref_lon_val, pole_lat, pole_lon, known_x, known_y, diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index 40e8e49..771164b 100644 --- a/src/wrf/metadecorators.py +++ b/src/wrf/metadecorators.py @@ -10,7 +10,7 @@ import numpy.ma as ma from .extension import _interpline from .util import (extract_vars, either, from_args, arg_location, is_coordvar, latlon_coordvars, to_np, - from_var, iter_left_indexes) + from_var, iter_left_indexes, is_mapping) from .coordpair import CoordPair from .py3compat import viewkeys, viewitems, py3range, ucode from .interputils import get_xy_z_params, get_xy, to_xy_coords @@ -586,8 +586,10 @@ def set_latlon_metadata(xy=False): @wrapt.decorator def func_wrapper(wrapped, instance, args, kwargs): - do_meta = from_args(wrapped, ("meta",), *args, **kwargs)["meta"] - + argvars = from_args(wrapped, ("wrfin", "meta"), *args, **kwargs) + # If it's a mapping, then this is handled as a special case in g_latlon + do_meta = (not is_mapping(argvars["wrfin"]) and argvars["meta"]) + if do_meta is None: do_meta = True diff --git a/src/wrf/util.py b/src/wrf/util.py index 7bf9927..2c0fd9e 100644 --- a/src/wrf/util.py +++ b/src/wrf/util.py @@ -3073,7 +3073,8 @@ def get_id(obj, prefix=''): # For sequences, the hashing string will be the list ID and the # path for the first file in the sequence if not is_mapping(obj): - _next = next(iter(obj)) + _obj = get_iterable(obj) + _next = next(iter(_obj)) return get_id(_next, prefix + str(id(obj))) # For each key in the mapping, recursively call get_id until diff --git a/test/comp_utest.py b/test/comp_utest.py index 864bb7a..de46a5f 100644 --- a/test/comp_utest.py +++ b/test/comp_utest.py @@ -642,15 +642,11 @@ def test_cape2d_1d(wrfnc): ter_follow = 1 result = cape_2d(p_hpa, tkel, qv, z, ter, psfc_hpa, ter_follow) - - print ("RESULT", result) ref = getvar(wrfnc, "cape_2d") ref = ref[(slice(None),) + col_idxs[1:]] - print ("REF", ref) - nt.assert_allclose(to_np(result), to_np(ref)) return func @@ -658,7 +654,7 @@ def test_cape2d_1d(wrfnc): if __name__ == "__main__": from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule, omp_set_dynamic, omp_get_num_procs, OMP_SCHED_STATIC) - omp_set_num_threads(omp_get_num_procs()-1) + omp_set_num_threads(omp_get_num_procs()//2) omp_set_schedule(OMP_SCHED_STATIC, 0) omp_set_dynamic(False) diff --git a/test/ipynb/Doc_Examples.ipynb b/test/ipynb/Doc_Examples.ipynb index 4b8d794..81ba5bd 100644 --- a/test/ipynb/Doc_Examples.ipynb +++ b/test/ipynb/Doc_Examples.ipynb @@ -1275,7 +1275,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", - "version": "2.7.13" + "version": "2.7.14" } }, "nbformat": 4, diff --git a/test/ipynb/Test_CTT.ipynb b/test/ipynb/Test_CTT.ipynb index 533ea13..c0c1a36 100644 --- a/test/ipynb/Test_CTT.ipynb +++ b/test/ipynb/Test_CTT.ipynb @@ -206,7 +206,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", - "version": "2.7.13" + "version": "2.7.14" } }, "nbformat": 4, diff --git a/test/ipynb/WRF_python_demo.ipynb b/test/ipynb/WRF_python_demo.ipynb index 977aaed..8f073d0 100644 --- a/test/ipynb/WRF_python_demo.ipynb +++ b/test/ipynb/WRF_python_demo.ipynb @@ -251,7 +251,7 @@ " return self\n", " \n", " def next(self):\n", - " if self._i > self._total:\n", + " if self._i >= self._total:\n", " raise StopIteration\n", " else:\n", " val = self.ncfile[self._i]\n", @@ -653,7 +653,7 @@ "metadata": {}, "outputs": [], "source": [ - "from wrf.latlon import xy_to_ll, ll_to_xy \n", + "from wrf import xy_to_ll, ll_to_xy \n", "\n", "a = xy_to_ll(ncfile, 400, 200)\n", "a1 = ll_to_xy(ncfile, a[0], a[1])\n", @@ -1134,7 +1134,7 @@ "metadata": {}, "outputs": [], "source": [ - "from wrf.latlon import xy_to_ll, ll_to_xy \n", + "from wrf import xy_to_ll, ll_to_xy \n", "\n", "a = xy_to_ll(ncfiles, 400, 200)\n", "a = xy_to_ll(ncfiles, [400,105], [200,205])\n", @@ -1196,6 +1196,27 @@ "print (\"\\n\")\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": null, @@ -1206,21 +1227,21 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 2", + "display_name": "Python 3", "language": "python", - "name": "python2" + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.13" + "pygments_lexer": "ipython3", + "version": "3.6.4" } }, "nbformat": 4, diff --git a/test/test_inputs.py b/test/test_inputs.py new file mode 100644 index 0000000..1fa986e --- /dev/null +++ b/test/test_inputs.py @@ -0,0 +1,181 @@ +import unittest as ut +import numpy.testing as nt +import numpy as np +import numpy.ma as ma +import os +import sys +import subprocess + +from wrf import (getvar, interplevel, interpline, vertcross, vinterp, + disable_xarray, xarray_enabled, to_np, + xy_to_ll, ll_to_xy, xy_to_ll_proj, ll_to_xy_proj, + extract_global_attrs, viewitems, CoordPair, ll_points) +from wrf.util import is_multi_file + +IN_DIR = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi" +TEST_FILES = [os.path.join(IN_DIR, "wrfout_d02_2005-08-28_00:00:00"), + os.path.join(IN_DIR, "wrfout_d02_2005-08-28_12:00:00"), + os.path.join(IN_DIR, "wrfout_d02_2005-08-29_00:00:00")] + +def wrfin_gen(wrf_in): + for x in wrf_in: + yield x + + +class wrf_in_iter_class(object): + def __init__(self, wrf_in): + self._wrf_in = wrf_in + self._total = len(wrf_in) + self._i = 0 + + def __iter__(self): + return self + + def next(self): + if self._i >= self._total: + raise StopIteration + else: + result = self._wrf_in[self._i] + self._i += 1 + return result + + # Python 3 + def __next__(self): + return self.next() + + +def make_test(varname, wrf_in): + def test(self): + x = getvar(wrf_in, varname, 0) + xa = getvar(wrf_in, varname, None) + + return test + +def make_interp_test(varname, wrf_in): + def test(self): + + # Only testing vinterp since other interpolation just use variables + if (varname == "vinterp"): + for timeidx in (0, None): + eth = getvar(wrf_in, "eth", timeidx=timeidx) + interp_levels = [850,500,5] + field = vinterp(wrf_in, + field=eth, + vert_coord="pressure", + interp_levels=interp_levels, + extrapolate=True, + field_type="theta-e", + timeidx=timeidx, + log_p=True) + else: + pass + + + return test + + +def make_latlon_test(testid, wrf_in): + def test(self): + if testid == "xy_out": + # Lats/Lons taken from NCL script, just hard-coding for now + lats = [-55, -60, -65] + lons = [25, 30, 35] + xy = ll_to_xy(wrf_in, lats, lons, timeidx=0) + xy = ll_to_xy(wrf_in, lats, lons, timeidx=None) + else: + i_s = np.asarray([10, 100, 150], int) + j_s = np.asarray([10, 100, 150], int) + ll = xy_to_ll(wrf_in, i_s, j_s, timeidx=0) + ll = xy_to_ll(wrf_in, i_s, j_s, timeidx=None) + + return test + + +class WRFVarsTest(ut.TestCase): + longMessage = True + +class WRFInterpTest(ut.TestCase): + longMessage = True + +class WRFLatLonTest(ut.TestCase): + longMessage = True + +if __name__ == "__main__": + from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule, + omp_set_dynamic, omp_get_num_procs, OMP_SCHED_STATIC) + omp_set_num_threads(omp_get_num_procs()-1) + omp_set_schedule(OMP_SCHED_STATIC, 0) + omp_set_dynamic(False) + + ignore_vars = [] # Not testable yet + wrf_vars = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz", + "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", + "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", + "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", + "wa", "uvmet10", "uvmet", "z", "cfrac", "zstag", "geopt_stag"] + interp_methods = ["interplevel", "vertcross", "interpline", "vinterp"] + latlon_tests = ["xy_out", "ll_out"] + + + for nc_lib in ("netcdf4", "pynio", "scipy"): + if nc_lib == "netcdf4": + try: + from netCDF4 import Dataset + except ImportError: + print ("netcdf4-python not installed") + continue + else: + test_in = [Dataset(x) for x in TEST_FILES] + elif nc_lib == "pynio": + try: + from Nio import open_file + except ImportError: + print ("PyNIO not installed") + continue + else: + test_in = [open_file(x +".nc", "r") for x in TEST_FILES] + elif nc_lib == "scipy": + try: + from scipy.io.netcdf import netcdf_file + except ImportError: + print ("scipy.io.netcdf not installed") + else: + test_in = [netcdf_file(x, mmap=False) for x in TEST_FILES] + + input0 = test_in[0] + input1 = test_in + input2 = tuple(input1) + input3 = wrfin_gen(test_in) + input4 = wrf_in_iter_class(test_in) + input5 = {"A" : input1, + "B" : input2} + input6 = {"A" : {"AA" : input1}, + "B" : {"BB" : input2}} + + for i,input in enumerate((input0, input1, input2, + input3, input4, input5, input6)): + for var in wrf_vars: + if var in ignore_vars: + continue + test_func1 = make_test(var, input) + + setattr(WRFVarsTest, "test_{0}_input{1}_{2}".format(nc_lib, + i, var), + test_func1) + + + for method in interp_methods: + test_interp_func1 = make_interp_test(method, input) + setattr(WRFInterpTest, + "test_{0}_input{1}_{2}".format(nc_lib, + i, method), + test_interp_func1) + + for testid in latlon_tests: + test_ll_func = make_latlon_test(testid, input) + test_name = "test_{0}_input{1}_{2}".format(nc_lib, i, testid) + setattr(WRFLatLonTest, test_name, test_ll_func) + + ut.main() + + \ No newline at end of file diff --git a/test/utests.py b/test/utests.py index 4962ef5..ed1d452 100644 --- a/test/utests.py +++ b/test/utests.py @@ -682,7 +682,7 @@ class WRFLatLonTest(ut.TestCase): if __name__ == "__main__": from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule, omp_set_dynamic, omp_get_num_procs, OMP_SCHED_STATIC) - omp_set_num_threads(omp_get_num_procs()-1) + omp_set_num_threads(omp_get_num_procs()//2) omp_set_schedule(OMP_SCHED_STATIC, 0) omp_set_dynamic(False) From 2c35edd7d7e42d5388615ae2b876345ae8b6502e Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Thu, 29 Mar 2018 16:00:28 -0600 Subject: [PATCH 9/9] Added ctt test --- test/ctt_test.py | 50 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 test/ctt_test.py diff --git a/test/ctt_test.py b/test/ctt_test.py new file mode 100644 index 0000000..6e97909 --- /dev/null +++ b/test/ctt_test.py @@ -0,0 +1,50 @@ +from netCDF4 import Dataset +import matplotlib.pyplot as plt +from matplotlib.cm import get_cmap +import cartopy.crs as crs +from cartopy.feature import NaturalEarthFeature + +from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords + +# Open the NetCDF file +ncfile = Dataset("/Users/ladwig/Documents/wrf_files/problem_files/cfrac_bug/wrfout_d02_1987-10-01_00:00:00") + +# Get the sea level pressure +ctt = getvar(ncfile, "ctt") + +# Get the latitude and longitude points +lats, lons = latlon_coords(ctt) + +# Get the cartopy mapping object +cart_proj = get_cartopy(ctt) + +# Create a figure +fig = plt.figure(figsize=(12,9)) +# Set the GeoAxes to the projection used by WRF +ax = plt.axes(projection=cart_proj) + +# Download and add the states and coastlines +states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none', + name='admin_1_states_provinces_shp') +ax.add_feature(states, linewidth=.5) +ax.coastlines('50m', linewidth=0.8) + +# Make the contour outlines and filled contours for the smoothed sea level pressure. +plt.contour(to_np(lons), to_np(lats), to_np(ctt), 10, colors="black", + transform=crs.PlateCarree()) +plt.contourf(to_np(lons), to_np(lats), to_np(ctt), 10, transform=crs.PlateCarree(), + cmap=get_cmap("jet")) + +# Add a color bar +plt.colorbar(ax=ax, shrink=.62) + +# Set the map limits. Not really necessary, but used for demonstration. +ax.set_xlim(cartopy_xlim(ctt)) +ax.set_ylim(cartopy_ylim(ctt)) + +# Add the gridlines +ax.gridlines(color="black", linestyle="dotted") + +plt.title("Cloud Top Temperature") + +plt.show()