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/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..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, 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,10 +8,14 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew !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 @@ -20,7 +25,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 @@ -58,12 +63,12 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew 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 @@ -74,21 +79,23 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew 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) @@ -98,17 +105,22 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew 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/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/fortran/wrffortran.pyf b/fortran/wrffortran.pyf index 7211ff6..be27452 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,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 @@ -360,7 +370,11 @@ 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 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) @@ -730,7 +744,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 +766,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 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 167d27c..42c4734 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()) @@ -754,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. @@ -766,6 +780,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,7 +790,11 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): ght, ter, outview, - haveqci) + pf, + haveqci, + fill_nocloud, + missing, + opt_thresh) return result @@ -858,7 +878,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 +899,7 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, extrap, vcor, logp, + tempout, missing, errstat, errmsg) 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/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/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 5442956..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" : [], + "ctt" : ["fill_nocloud", "missing", "opt_thresh", "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): 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 " 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 +}