From cc86d68a35efb5de440426e595c4286c875d1108 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 16 Mar 2018 16:34:44 -0600 Subject: [PATCH] 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 +}