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 +}