diff --git a/doc/source/_static/images/cartopy_slp.png b/doc/source/_static/images/cartopy_slp.png new file mode 100644 index 0000000..d107bd0 Binary files /dev/null and b/doc/source/_static/images/cartopy_slp.png differ diff --git a/doc/source/_static/images/matthew.png b/doc/source/_static/images/matthew.png new file mode 100644 index 0000000..c2362e8 Binary files /dev/null and b/doc/source/_static/images/matthew.png differ diff --git a/doc/source/api.rst b/doc/source/api.rst index a205c49..544c28f 100644 --- a/doc/source/api.rst +++ b/doc/source/api.rst @@ -3,7 +3,7 @@ API Reference .. toctree:: - :maxdepth: 1 + :maxdepth: 2 user_api/index.rst internal_api/index.rst diff --git a/doc/source/basic_usage.rst b/doc/source/basic_usage.rst new file mode 100644 index 0000000..e89dc84 --- /dev/null +++ b/doc/source/basic_usage.rst @@ -0,0 +1,1176 @@ +How To Use +============ + +Basic Usage +---------------- + +.. _diagnostic-usage: + +Computing Diagnostic Variables +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The primary use for the :meth:`wrf.getvar` function is to return diagnostic +variables that require a calculation, since WRF does not produce these variables +natively. These diagnostics include CAPE, storm relative helicity, +omega, sea level pressure, etc. A table of all available diagnostics can be +found here: :ref:`diagnostic-table`. + +In the example below, sea level pressure is calculated and printed. + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + # Get the Sea Level Pressure + slp = getvar(ncfile, "slp") + + print(slp) + +Result: + +.. code-block:: none + + + array([[ 1012.22033691, 1012.29815674, 1012.24786377, ..., + 1010.13201904, 1009.93231201, 1010.06707764], + [ 1012.43286133, 1012.44476318, 1012.33666992, ..., + 1010.1072998 , 1010.10845947, 1010.04760742], + [ 1012.39544678, 1012.38085938, 1012.41705322, ..., + 1010.22937012, 1010.05596924, 1010.02679443], + ..., + [ 1009.0423584 , 1009.06921387, 1008.98779297, ..., + 1019.19281006, 1019.14434814, 1019.1105957 ], + [ 1009.22485352, 1009.07513428, 1008.98638916, ..., + 1019.07189941, 1019.04266357, 1019.0612793 ], + [ 1009.18896484, 1009.1071167 , 1008.97979736, ..., + 1018.91778564, 1018.95684814, 1019.04748535]], dtype=float32) + Coordinates: + XLONG (south_north, west_east) float32 -122.72 -122.693 -122.666 ... + XLAT (south_north, west_east) float32 21.1381 21.1451 21.1521 ... + Time datetime64[ns] 2016-10-07 + * south_north (south_north) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ... + * west_east (west_east) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + Attributes: + FieldType: 104 + MemoryOrder: XY + description: sea level pressure + units: hPa + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + +.. _extract_ncvars: + +Extracting WRF NetCDF Variables +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In addition to computing diagnostic variables (see :ref:`diagnostic-usage`), +the :meth:`wrf.getvar` function can be used to extract regular WRF-ARW output +NetCDF variables. + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + p = getvar(ncfile, "P") + + print(p) + +Result: + +.. code-block:: none + + + array([[[ 1.21753906e+03, 1.22532031e+03, 1.22030469e+03, ..., + 1.00760156e+03, 9.87640625e+02, 1.00111719e+03], + [ 1.23877344e+03, 1.24004688e+03, 1.22926562e+03, ..., + 1.00519531e+03, 1.00529688e+03, 9.99171875e+02], + [ 1.23503906e+03, 1.23367188e+03, 1.23731250e+03, ..., + 1.01739844e+03, 1.00005469e+03, 9.97093750e+02], + ..., + [ 1.77978516e+00, 1.77050781e+00, 1.79003906e+00, ..., + 4.22949219e+00, 4.25659180e+00, 4.13647461e+00], + [ 1.73291016e+00, 1.76879883e+00, 1.77978516e+00, ..., + 4.24047852e+00, 4.24707031e+00, 4.13549805e+00], + [ 1.71533203e+00, 1.65722656e+00, 1.67480469e+00, ..., + 4.06884766e+00, 4.03637695e+00, 4.04785156e+00]]], dtype=float32) + Coordinates: + XLONG (south_north, west_east) float32 -122.72 -122.693 -122.666 ... + XLAT (south_north, west_east) float32 21.1381 21.1451 21.1521 ... + Time datetime64[ns] 2016-10-07 + * bottom_top (bottom_top) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + * south_north (south_north) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ... + * west_east (west_east) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + Attributes: + FieldType: 104 + MemoryOrder: XYZ + description: perturbation pressure + units: Pa + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + + +Disabling xarray and metadata +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Sometimes you just want a regular numpy array and don't care about metadata. +This is often the case when you are working with compiled extensions. Metadata +can be disabled in one of two ways. + +#. disable xarray completely +#. set the *meta* function parameter to False. + +The example below illustrates both. + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, disable_xarray + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + # Disable xarray completely + disable_xarray() + p_no_meta = getvar(ncfile, "P") + print (type(p_no_meta)) + enable_xarray() + + # Disable by using the meta parameter + p_no_meta = getvar(ncfile, "P", meta=False) + print (type(p_no_meta)) + +Result: + +.. code-block:: none + + + + +Extracting a Numpy Array from a DataArray +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +If you need to convert an :class:`xarray.DataArray` to a :class:`numpy.ndarray`, +wrf-python provides the :meth:`wrf.npvalues` function for this purpose. Although +an :class:`xarray.DataArary` object already contains the +:attr:`xarray.DataArray.values` attribute to extract the Numpy array, there is a +problem when working with compiled extensions. The behavior for xarray (and pandas) +is to convert missing/fill values to NaN, which may cause crashes when working +with compiled extensions. Also, some existing code may be designed to work with +:class:`numpy.ma.MaskedArray`, and numpy arrays with NaN may not work with it. + +The :meth:`wrf.npvalues` function does the following: + +#. If no missing/fill values are used, :meth:`wrf.npvalues` simply returns the + :attr:`xarray.DataArray.values` attribute. + +#. If missing/fill values are used, then :meth:`wrf.npvalues` replaces the NaN + values with the _FillValue found in the :attr:`xarray.DataArray.attrs` + attribute (required) and a :class:`numpy.ma.MaskedArray` is returned. + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + # Get the Sea Level Pressure + cape_3d = getvar(ncfile, "cape_3d") + + cape_3d_ndarray = npvalues(cape_3d) + + print(type(cape_3d_ndarray)) + + +Result: + +.. code-block:: none + + + + +Sequences of Files +---------------------- + +Combining Multiple Files Using the 'cat' Method +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The 'cat' (concatenate) method aggregates all files in the sequence along the +'Time' dimension, which will be the leftmost dimension for the output array. +To include all of the times, in all of the files, in the output array, set the +*timeidx* parameter to :data:`wrf.ALL_TIMES` (an alias for None). If a single +value is specified for *timeidx*, then the time index is assumed to be taken from +the concatenation of all times for all files. + +It is import to note that no sorting is performed in the :meth:`wrf.getvar` +routine, so all files in the sequence must be sorted prior to calling this +function. + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, ALL_TIMES + + + # Creating a simple test list with three timesteps + wrflist = [Dataset("wrfout_d01_2016-10-07_00_00_00"), + Dataset("wrfout_d01_2016-10-07_01_00_00"), + Dataset("wrfout_d01_2016-10-07_02_00_00")] + + # Extract the 'P' variable for all times + p_cat = getvar(wrflist, "P", timeidx=ALL_TIMES, method="cat") + + print(p_cat) + +Result: + +.. code-block:: none + + + array([[[[ 1.21753906e+03, 1.22532031e+03, 1.22030469e+03, ..., + 1.00760156e+03, 9.87640625e+02, 1.00111719e+03], + [ 1.23877344e+03, 1.24004688e+03, 1.22926562e+03, ..., + 1.00519531e+03, 1.00529688e+03, 9.99171875e+02], + [ 1.23503906e+03, 1.23367188e+03, 1.23731250e+03, ..., + 1.01739844e+03, 1.00005469e+03, 9.97093750e+02], + ..., + [ 1.77978516e+00, 1.77050781e+00, 1.79003906e+00, ..., + 4.22949219e+00, 4.25659180e+00, 4.13647461e+00], + [ 1.73291016e+00, 1.76879883e+00, 1.77978516e+00, ..., + 4.24047852e+00, 4.24707031e+00, 4.13549805e+00], + [ 1.71533203e+00, 1.65722656e+00, 1.67480469e+00, ..., + 4.06884766e+00, 4.03637695e+00, 4.04785156e+00]]]], dtype=float32) + Coordinates: + XLONG (south_north, west_east) float32 -122.72 -122.693 -122.666 ... + XLAT (south_north, west_east) float32 21.1381 21.1451 21.1521 ... + * Time (Time) datetime64[ns] 2016-10-07 2016-10-07 2016-10-07 + * bottom_top (bottom_top) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + * south_north (south_north) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ... + * west_east (west_east) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + datetime (Time) datetime64[ns] 2016-10-07T00:00:00 ... + Attributes: + FieldType: 104 + MemoryOrder: XYZ + description: perturbation pressure + units: Pa + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + + +Combining Multiple Files Using the 'join' Method +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The 'join' method combines a sequence of files by adding a new leftmost +dimension for the file/sequence index. In situations where there are multiple +files with multiple times, and the last file contains less times than the +previous files, the remaining arrays will be arrays filled with missing values. +There are checks in place within the wrf-python algorithms to look for these missing +arrays, but be careful when calling compiled routines outside of wrf-python. + +In most cases, *timeidx* parameter should be set to :data:`wrf.ALL_TIMES`. If +a *timeidx* value is specified, then this time index is used when extracting the +variable from each file. In cases where there are multiple files with multiple +time steps, this is probably nonsensical, since the nth time index for each +file represents a different time. + +In general, join is rarely used, so the concatenate method should be used +for most cases. + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, ALL_TIMES + + + # Creating a simple test list with three timesteps + wrflist = [Dataset("wrfout_d01_2016-10-07_00_00_00"), + Dataset("wrfout_d01_2016-10-07_01_00_00"), + Dataset("wrfout_d01_2016-10-07_02_00_00")] + + # Extract the 'P' variable for all times + p_join = getvar(wrflist, "P", timeidx=ALL_TIMES, method="join") + + print(p_join) + +Result: + +.. code-block:: none + + + array([[[[ 1.21753906e+03, 1.22532031e+03, 1.22030469e+03, ..., + 1.00760156e+03, 9.87640625e+02, 1.00111719e+03], + [ 1.23877344e+03, 1.24004688e+03, 1.22926562e+03, ..., + 1.00519531e+03, 1.00529688e+03, 9.99171875e+02], + [ 1.23503906e+03, 1.23367188e+03, 1.23731250e+03, ..., + 1.01739844e+03, 1.00005469e+03, 9.97093750e+02], + ..., + [ 1.77978516e+00, 1.77050781e+00, 1.79003906e+00, ..., + 4.22949219e+00, 4.25659180e+00, 4.13647461e+00], + [ 1.73291016e+00, 1.76879883e+00, 1.77978516e+00, ..., + 4.24047852e+00, 4.24707031e+00, 4.13549805e+00], + [ 1.71533203e+00, 1.65722656e+00, 1.67480469e+00, ..., + 4.06884766e+00, 4.03637695e+00, 4.04785156e+00]]]], dtype=float32) + Coordinates: + XLONG (south_north, west_east) float32 -122.72 -122.693 -122.666 ... + XLAT (south_north, west_east) float32 21.1381 21.1451 21.1521 ... + * bottom_top (bottom_top) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + * south_north (south_north) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ... + * west_east (west_east) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + * file (file) int64 0 1 2 + datetime (file) datetime64[ns] 2016-10-07T00:00:00 ... + Time int64 0 + Attributes: + FieldType: 104 + MemoryOrder: XYZ + description: perturbation pressure + units: Pa + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + + +Note how the 'Time' dimension was replaced with the 'file' dimension, due to the +numpy's automatic squeezing of the single 'Time' dimension. To maintain the +'Time' dimension, set the *squeeze* parameter to False. + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, ALL_TIMES + + + # Creating a simple test list with three timesteps + wrflist = [Dataset("wrfout_d01_2016-10-07_00_00_00"), + Dataset("wrfout_d01_2016-10-07_01_00_00"), + Dataset("wrfout_d01_2016-10-07_02_00_00")] + + # Extract the 'P' variable for all times + p_join = getvar(wrflist, "P", timeidx=ALL_TIMES, method="join", squeeze=False) + + print(p_join) + +Result + +.. code-block:: none + + + array([[[[[ 1.21753906e+03, 1.22532031e+03, 1.22030469e+03, ..., + 1.00760156e+03, 9.87640625e+02, 1.00111719e+03], + [ 1.23877344e+03, 1.24004688e+03, 1.22926562e+03, ..., + 1.00519531e+03, 1.00529688e+03, 9.99171875e+02], + [ 1.23503906e+03, 1.23367188e+03, 1.23731250e+03, ..., + 1.01739844e+03, 1.00005469e+03, 9.97093750e+02], + ..., + [ 1.77978516e+00, 1.77050781e+00, 1.79003906e+00, ..., + 4.22949219e+00, 4.25659180e+00, 4.13647461e+00], + [ 1.73291016e+00, 1.76879883e+00, 1.77978516e+00, ..., + 4.24047852e+00, 4.24707031e+00, 4.13549805e+00], + [ 1.71533203e+00, 1.65722656e+00, 1.67480469e+00, ..., + 4.06884766e+00, 4.03637695e+00, 4.04785156e+00]]]]], dtype=float32) + Coordinates: + XLONG (south_north, west_east) float32 -122.72 -122.693 -122.666 ... + XLAT (south_north, west_east) float32 21.1381 21.1451 21.1521 ... + * bottom_top (bottom_top) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + * south_north (south_north) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ... + * west_east (west_east) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + * file (file) int64 0 1 2 + datetime (file, Time) datetime64[ns] 2016-10-07T00:00:00 ... + * Time (Time) int64 0 + Attributes: + FieldType: 104 + MemoryOrder: XYZ + description: perturbation pressure + units: Pa + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + + +Dictionaries of WRF File Sequences +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Dictionaries can also be used as input to the :meth:`wrf.getvar` functions. +This can be useful when working with ensembles. However, all WRF files in the +dictionary must have the same dimensions. The result is an array where the +leftmost dimension is the keys from the dictionary. Nested dictionaries +are allowed. + +The *method* argument is used to describe how each sequence in the dictionary +will be combined. + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, ALL_TIMES + + wrf_dict = {"ens1" : [Dataset("ens1/wrfout_d01_2016-10-07_00_00_00"), + Dataset("ens1/wrfout_d01_2016-10-07_01_00_00"), + Dataset("ens1/wrfout_d01_2016-10-07_02_00_00")], + "ens2" : [Dataset("ens2/wrfout_d01_2016-10-07_00_00_00"), + Dataset("ens2/wrfout_d01_2016-10-07_01_00_00"), + Dataset("ens2/wrfout_d01_2016-10-07_02_00_00")] + } + + p = getvar(wrf_dict, "P", timeidx=ALL_TIMES) + + print(p) + +Result: + +.. code-block:: none + + + array([[[[[ 1.21753906e+03, 1.22532031e+03, 1.22030469e+03, ..., + 1.00760156e+03, 9.87640625e+02, 1.00111719e+03], + [ 1.23877344e+03, 1.24004688e+03, 1.22926562e+03, ..., + 1.00519531e+03, 1.00529688e+03, 9.99171875e+02], + [ 1.23503906e+03, 1.23367188e+03, 1.23731250e+03, ..., + 1.01739844e+03, 1.00005469e+03, 9.97093750e+02], + ..., + [ 1.77978516e+00, 1.77050781e+00, 1.79003906e+00, ..., + 4.22949219e+00, 4.25659180e+00, 4.13647461e+00], + [ 1.73291016e+00, 1.76879883e+00, 1.77978516e+00, ..., + 4.24047852e+00, 4.24707031e+00, 4.13549805e+00], + [ 1.71533203e+00, 1.65722656e+00, 1.67480469e+00, ..., + 4.06884766e+00, 4.03637695e+00, 4.04785156e+00]]]]], dtype=float32) + Coordinates: + XLONG (south_north, west_east) float32 -122.72 -122.693 -122.666 ... + XLAT (south_north, west_east) float32 21.1381 21.1451 21.1521 ... + * Time (Time) datetime64[ns] 2016-10-07T00:00:00 ... + * bottom_top (bottom_top) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + * south_north (south_north) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ... + * west_east (west_east) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + datetime (Time) datetime64[ns] 2016-10-07T00:00:00 ... + * key_0 (key_0) + array([[ 5882.16992188, 5881.87939453, 5881.81005859, ..., + 5890.14501953, 5890.23583984, 5890.33349609], + [ 5882.71777344, 5882.17529297, 5882.1171875 , ..., + 5890.37695312, 5890.38525391, 5890.27978516], + [ 5883.32177734, 5882.47119141, 5882.34130859, ..., + 5890.48339844, 5890.42871094, 5890.17724609], + ..., + [ 5581.45800781, 5580.46826172, 5579.32617188, ..., + 5788.93554688, 5788.70507812, 5788.64453125], + [ 5580.32714844, 5579.51611328, 5578.34863281, ..., + 5788.15869141, 5787.87304688, 5787.65527344], + [ 5579.64404297, 5578.30957031, 5576.98632812, ..., + 5787.19384766, 5787.10888672, 5787.06933594]], dtype=float32) + Coordinates: + XLONG (south_north, west_east) float32 -122.72 -122.693 -122.666 ... + XLAT (south_north, west_east) float32 21.1381 21.1451 21.1521 ... + Time datetime64[ns] 2016-10-07 + * south_north (south_north) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ... + * west_east (west_east) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + Attributes: + FieldType: 104 + units: m + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + level: 500 hPa + missing_value: 9.96920996839e+36 + _FillValue: 9.96920996839e+36 + + +.. _vert_cross_interp: + +Vertical Cross Sections +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The :meth:`wrf.vertcross` function is used to create vertical cross sections. +To define a cross section, a start point and an end point needs to be specified. +Alternatively, a pivot point and an angle may be used. The start point, +end point, and pivot point are specified using a :class:`wrf.CoordPair` object, +and coordinates can either be in grid (x,y) coordinates or (latitude,longitude) +coordinates. When using (latitude,longitude) coordinates, a NetCDF file object or +a :class:`wrf.WrfProj` object must be provided. + +The vertical levels can also be specified using the *levels* parameter. If +not specified, then approximately 100 levels will be chosen in 1% increments. + +Example Using Start Point and End Point +***************************************** + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, vertcross, CoordPair + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + # Get the geopotential height (m) and pressure (hPa). + z = getvar(ncfile, "z") + p = getvar(ncfile, "pressure") + + # Define a start point and end point in grid coordinates + start_point = CoordPair(x=0, y=(z.shape[-2]-1)//2) + end_point = CoordPair(x=-1, y=(z.shape[-2]-1)//2) + + # Calculate the vertical cross section. By setting latlon to True, this + # also calculates the latitude and longitude coordinates along the cross + # section line and adds them to the 'xy_loc' metadata to help with plotting. + p_vert = vertcross(p, z, start_point=start_point, end_point=end_point, latlon=True) + + print(p_vert) + +Result: + +.. code-block:: none + + + array([[ nan, nan, nan, ..., nan, + nan, nan], + [ 989.66168213, 989.66802979, 989.66351318, ..., 988.05737305, + 987.99151611, 987.96917725], + [ 959.49450684, 959.50109863, 959.50030518, ..., 958.96948242, + 958.92980957, 958.89294434], + ..., + [ 24.28092003, 24.27359581, 24.27034378, ..., 24.24800491, + 24.2486496 , 24.24947357], + [ 23.2868309 , 23.27933884, 23.27607918, ..., 23.25231361, + 23.2530098 , 23.25384521], + [ nan, nan, nan, ..., nan, + nan, nan]], dtype=float32) + Coordinates: + Time datetime64[ns] 2016-10-07 + xy_loc (idx) object CoordPair(x=0.0, y=529.0, lat=34.5279502869, lon=-127.398925781) ... + * vertical (vertical) float32 0.0 261.828 523.656 785.484 1047.31 1309.14 ... + * idx (idx) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ... + Attributes: + FieldType: 104 + description: pressure + units: hPa + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + orientation: (0.0, 529.0) to (1797.0, 529.0) + missing_value: 9.96920996839e+36 + _FillValue: 9.96920996839e+36 + + +Example Using Pivot Point and Angle +************************************* + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, vertcross, CoordPair + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + # Get the geopotential height (m) and pressure (hPa). + z = getvar(ncfile, "z") + p = getvar(ncfile, "pressure") + + # Define a pivot point and angle in grid coordinates, with the + # pivot point being the center of the grid. + pivot_point = CoordPair(x=(z.shape[-1]-1)//2, y=(z.shape[-2]-1)//2) + angle = 90.0 + + # Calculate the vertical cross section. By setting latlon to True, this + # also calculates the latitude and longitude coordinates along the line + # and adds them to the metadata to help with plotting labels. + p_vert = vertcross(p, z, pivot_point=pivot_point, angle=angle, latlon=True) + + print (p_vert) + +Result: + +.. code-block:: none + + + array([[ nan, nan, nan, ..., nan, + nan, nan], + [ 989.66168213, 989.66802979, 989.66351318, ..., 988.05737305, + 987.99151611, 987.96917725], + [ 959.49450684, 959.50109863, 959.50030518, ..., 958.96948242, + 958.92980957, 958.89294434], + ..., + [ 24.28092003, 24.27359581, 24.27034378, ..., 24.24800491, + 24.2486496 , 24.24947357], + [ 23.2868309 , 23.27933884, 23.27607918, ..., 23.25231361, + 23.2530098 , 23.25384521], + [ nan, nan, nan, ..., nan, + nan, nan]], dtype=float32) + Coordinates: + Time datetime64[ns] 2016-10-07 + xy_loc (idx) object CoordPair(x=0.0, y=529.0, lat=34.5279502869, lon=-127.398925781) ... + * vertical (vertical) float32 0.0 261.828 523.656 785.484 1047.31 1309.14 ... + * idx (idx) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ... + Attributes: + FieldType: 104 + description: pressure + units: hPa + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + orientation: (0.0, 529.0) to (1797.0, 529.0) ; center=CoordPair(x=899.0, y=529.0) ; angle=90.0 + missing_value: 9.96920996839e+36 + _FillValue: 9.96920996839e+36 + + +Example Using Lat/Lon Coordinates +************************************* + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, vertcross, CoordPair + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + # Get the geopotential height (m) and pressure (hPa). + z = getvar(ncfile, "z") + p = getvar(ncfile, "pressure") + lats = getvar(ncfile, "lat") + lons = getvar(ncfile, "lon") + + # Making the same horizontal line, but with lats/lons + start_lat = lats[(lats.shape[-2]-1)//2, 0] + end_lat = lats[(lats.shape[-2]-1)//2, -1] + start_lon = lons[(lats.shape[-2]-1)//2, 0] + end_lon = lons[(lats.shape[-2]-1)//2, -1] + + # Cross section line using start_point and end_point. + start_point = CoordPair(lat=start_lat, lon=start_lon) + end_point = CoordPair(lat=end_lat, lon=end_lon) + + # When using lat/lon coordinates, you must supply a netcdf file object, or a + # projection object. + p_vert = vertcross(p, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True) + print(p_vert) + +Result: + +.. code-block:: none + + + array([[ nan, nan, nan, ..., nan, + nan, nan], + [ 989.66168213, 989.66802979, 989.66351318, ..., 988.05737305, + 987.99151611, 987.96917725], + [ 959.49450684, 959.50109863, 959.50030518, ..., 958.96948242, + 958.92980957, 958.89294434], + ..., + [ 24.28092003, 24.27359581, 24.27034378, ..., 24.24800491, + 24.2486496 , 24.24947357], + [ 23.2868309 , 23.27933884, 23.27607918, ..., 23.25231361, + 23.2530098 , 23.25384521], + [ nan, nan, nan, ..., nan, + nan, nan]], dtype=float32) + Coordinates: + Time datetime64[ns] 2016-10-07 + xy_loc (idx) object CoordPair(x=0.0, y=529.0, lat=34.5279502869, lon=-127.398925781) ... + * vertical (vertical) float32 0.0 261.828 523.656 785.484 1047.31 1309.14 ... + * idx (idx) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ... + Attributes: + FieldType: 104 + description: pressure + units: hPa + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + orientation: (0.0, 529.0) to (1797.0, 529.0) + missing_value: 9.96920996839e+36 + _FillValue: 9.96920996839e+36 + + +Example Using Specified Vertical Levels +***************************************** + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, vertcross, CoordPair + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + # Get the geopotential height (m) and pressure (hPa). + z = getvar(ncfile, "z") + p = getvar(ncfile, "pressure") + lats = getvar(ncfile, "lat") + lons = getvar(ncfile, "lon") + + # Making the same horizontal line, but with lats/lons + start_lat = lats[(lats.shape[-2]-1)//2, 0] + end_lat = lats[(lats.shape[-2]-1)//2, -1] + start_lon = lons[(lats.shape[-2]-1)//2, 0] + end_lon = lons[(lats.shape[-2]-1)//2, -1] + + # Pressure using start_point and end_point. These were obtained using + start_point = CoordPair(lat=start_lat, lon=start_lon) + end_point = CoordPair(lat=end_lat, lon=end_lon) + + # Specify vertical levels + levels = [1000., 2000., 3000.] + + # Calculate the cross section + p_vert = vertcross(p, z, wrfin=ncfile, levels=levels, start_point=start_point, end_point=end_point, latlon=True) + + print(p_vert) + +Result: + +.. code-block:: none + + + array([[ 906.375 , 906.38043213, 906.39367676, ..., 907.6661377 , + 907.63006592, 907.59191895], + [ 804.24737549, 804.26885986, 804.28076172, ..., 806.98632812, + 806.95556641, 806.92608643], + [ 713.24578857, 713.2722168 , 713.27886963, ..., 716.09594727, + 716.06610107, 716.03503418]], dtype=float32) + Coordinates: + Time datetime64[ns] 2016-10-07 + xy_loc (idx) object CoordPair(x=0.0, y=529.0, lat=34.5279502869, lon=-127.398925781) ... + * vertical (vertical) float32 1000.0 2000.0 3000.0 + * idx (idx) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ... + Attributes: + FieldType: 104 + description: pressure + units: hPa + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + orientation: (0.0, 529.0) to (1797.0, 529.0) + missing_value: 9.96920996839e+36 + _FillValue: 9.96920996839e+36 + + +Interpolating Two-Dimensional Fields to a Line +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Two-dimensional fields can be interpolated along a line, in a manner similar to +the vertical cross section (see :ref:`vert_cross_interp`), using the +:meth:`wrf.interpline` function. To define the line +to interpolate along, a start point and an end point needs to be specified. +Alternatively, a pivot point and an angle may be used. The start point, +end point, and pivot point are specified using a :class:`wrf.CoordPair` object, +and coordinates can either be in grid (x,y) coordinates or (latitude,longitude) +coordinates. When using (latitude,longitude) coordinates, a NetCDF file object or +a :class:`wrf.WrfProj` object must also be provided. + +Example Using Start Point and End Point +***************************************** + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, interpline, CoordPair + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + # Get the 2m temperature + t2 = getvar(ncfile, "T2") + + # Create a south-north line in the center of the domain using + # start point and end point + start_point = CoordPair(x=(t2.shape[-1]-1)//2, y=0) + end_point = CoordPair(x=(t2.shape[-1]-1)//2, y=-1) + + # Calculate the vertical cross section. By setting latlon to True, this + # also calculates the latitude and longitude coordinates along the line + # and adds them to the metadata to help with plotting labels. + t2_line = interpline(t2, start_point=start_point, end_point=end_point, latlon=True) + + print(t2_line, "\n") + +Result: + +.. code-block:: none + + + array([ 302.07214355, 302.08505249, 302.08688354, ..., 279.18557739, + 279.1998291 , 279.23132324], dtype=float32) + Coordinates: + Time datetime64[ns] 2016-10-07 + xy_loc (line_idx) object CoordPair(x=899.0, y=0.0, lat=24.3645858765, lon=-97.5) ... + * line_idx (line_idx) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ... + Attributes: + FieldType: 104 + description: TEMP at 2 M + units: K + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + orientation: (899.0, 0.0) to (899.0, 1057.0) + + +Example Using Pivot Point and Angle +***************************************** + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, interpline, CoordPair + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + # Get the 2m temperature + t2 = getvar(ncfile, "T2") + + # Create a south-north line using pivot point and angle + pivot_point = CoordPair((t2.shape[-1]-1)//2, (t2.shape[-2]-1)//2) + angle = 0.0 + + # Calculate the vertical cross section. By setting latlon to True, this + # also calculates the latitude and longitude coordinates along the line + # and adds them to the metadata to help with plotting labels. + t2_line = interpline(t2, start_point=start_point, end_point=end_point, latlon=True) + + print(t2_line, "\n") + +Result: + +.. code-block:: none + + + array([ 302.07214355, 302.08505249, 302.08688354, ..., 279.18557739, + 279.1998291 , 279.23132324], dtype=float32) + Coordinates: + Time datetime64[ns] 2016-10-07 + xy_loc (line_idx) object CoordPair(x=899.0, y=0.0, lat=24.3645858765, lon=-97.5) ... + * line_idx (line_idx) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ... + Attributes: + FieldType: 104 + description: TEMP at 2 M + units: K + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + orientation: (899.0, 0.0) to (899.0, 1057.0) ; center=CoordPair(x=899, y=529) ; angle=0.0 + + +Example Using Lat/Lon Coordinates +************************************* + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, interpline, CoordPair + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + t2 = getvar(ncfile, "T2") + lats = getvar(ncfile, "lat") + lons = getvar(ncfile, "lon") + + # Select the latitude,longitude points for a vertical line through + # the center of the domain. + start_lat = lats[0, (lats.shape[-1]-1)//2] + end_lat = lats[-1, (lats.shape[-1]-1)//2] + start_lon = lons[0, (lons.shape[-1]-1)//2] + end_lon = lons[-1, (lons.shape[-1]-1)//2] + + # Create the CoordPairs + start_point = CoordPair(lat=start_lat, lon=start_lon) + end_point = CoordPair(lat=end_lat, lon=end_lon) + + # Calculate the vertical cross section. By setting latlon to True, this + # also calculates the latitude and longitude coordinates along the line + # and adds them to the metadata to help with plotting labels. + t2_line = interpline(t2, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True) + + print (t2_line) + +Result: + +.. code-block:: none + + + array([ 302.07214355, 302.08505249, 302.08688354, ..., 279.18557739, + 279.1998291 , 279.23132324], dtype=float32) + Coordinates: + Time datetime64[ns] 2016-10-07 + xy_loc (line_idx) object CoordPair(x=899.0, y=0.0, lat=24.3645858765, lon=-97.5) ... + * line_idx (line_idx) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 ... + Attributes: + FieldType: 104 + description: TEMP at 2 M + units: K + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + orientation: (899.0, 0.0) to (899.0, 1057.0) + + +Interpolating a 3D Field to a Surface Type +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The :meth:`wrf.vinterp` is used to interpolate a field to a type of surface. +The available surfaces are pressure, geopotential height, theta, and theta-e. +The surface levels to interpolate also need to be specified. + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, interpline, CoordPair + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + # Interpolate tk to theta-e levels + interp_levels = [200, 300, 500, 1000] + + interp_field = vinterp(ncfile, + field=tk, + vert_coord="eth", + interp_levels=interp_levels, + extrapolate=True, + field_type="tk", + log_p=True) + + print(interp_field) + +Result: + +.. code-block:: none + + + array([[[ 296.12872314, 296.1166687 , 296.08905029, ..., 301.71026611, + 301.67956543, 301.67791748], + [ 296.11352539, 295.95581055, 295.91555786, ..., 301.63052368, + 301.62905884, 301.65887451], + [ 296.07556152, 295.91577148, 295.88214111, ..., 301.61499023, + 301.60287476, 301.63961792], + ..., + [ 219.11134338, 219.08581543, 219.08602905, ..., 218.29879761, + 218.30923462, 218.3787384 ], + [ 219.09260559, 219.07765198, 219.08340454, ..., 218.2855072 , + 218.30444336, 218.37931824], + [ 219.07936096, 219.08181763, 219.10089111, ..., 218.31173706, + 218.34288025, 218.3687439 ]]], dtype=float32) + Coordinates: + XLONG (south_north, west_east) float32 -122.72 -122.693 -122.666 ... + XLAT (south_north, west_east) float32 21.1381 21.1451 21.1521 ... + Time datetime64[ns] 2016-10-07 + * south_north (south_north) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ... + * west_east (west_east) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ... + * interp_level (interp_level) int64 200 300 500 1000 + Attributes: + FieldType: 104 + MemoryOrder: XYZ + description: temperature + units: K + stagger: + coordinates: XLONG XLAT + projection: LambertConformal(bottom_left=(21.138123, -122.71953), + top_right=(47.843636, -60.901367), stand_lon=-97.5, + moad_cen_lat=38.5000038147, truelat1=38.5, truelat2=38.5, + pole_lat=90.0, pole_lon=0.0) + vert_interp_type: eth + + +Lat/Lon <-> XY Routines +-------------------------- + +wrf-python includes a set of routines for converting back and forth between +latitude,longitude space and x,y space. The methods are :meth:`wrf.xy_to_ll`, +:meth:`wrf.xy_to_ll_proj`, :meth:`wrf.ll_to_xy`, :meth:`wrf.ll_to_xy_proj`. +The *latitude*, *longitude*, *x*, and *y* parameters to these methods +can contain sequences if multiple points are desired to be converted. + +Example With Single Coordinates +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + from wrf import getvar, interpline, CoordPair, xy_to_ll, ll_to_xy + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + lat_lon = xy_to_ll(ncfile, 400, 200) + + print(lat_lon) + + x_y = ll_to_xy(ncfile, lat_lon[0], lat_lon[1]) + + print (x_y) + +Result: + +.. code-block:: none + + + array([ 28.55816408, -112.67827617]) + Coordinates: + * lat_lon (lat_lon) + array([400, 200]) + Coordinates: + latlon_coord object CoordPair(lat=28.5581640822, lon=-112.678276173) + * x_y (x_y) + array([[ 28.55816408, 27.03835783], + [-112.67827617, -121.36392174]]) + Coordinates: + * lat_lon (lat_lon) + array([[400, 105], + [200, 205]]) + Coordinates: + latlon_coord (idx) object CoordPair(lat=28.5581640822, lon=-112.678276173) ... + * x_y (x_y) `_:: + + $ conda install -c bladwig wrf-python + + +Installing via Source Code +-------------------------- + +Installation via source code will require a Fortran and C compiler in order +to run f2py. You can get them +`here `_. + +The source code is available via github: + +https://github.com/NCAR/wrf-python + +To install, change to the wrf-python directory and run:: + + $ pip install . + + diff --git a/doc/source/internal_api/index.rst b/doc/source/internal_api/index.rst index 9bbba0f..11fa312 100644 --- a/doc/source/internal_api/index.rst +++ b/doc/source/internal_api/index.rst @@ -1,8 +1,6 @@ Internal API ============= -------------- - Routines ------------- diff --git a/doc/source/new.rst b/doc/source/new.rst new file mode 100644 index 0000000..9ae2363 --- /dev/null +++ b/doc/source/new.rst @@ -0,0 +1,10 @@ +What's New +=========== + +v1.0a3 +----------- + +- Alpha release 3 +- Added docstrings +- Now uses CoordPair for cross sections so that lat/lon can be used +- Renamed some functions and or arguments diff --git a/doc/source/plot.rst b/doc/source/plot.rst new file mode 100644 index 0000000..50452d4 --- /dev/null +++ b/doc/source/plot.rst @@ -0,0 +1,91 @@ +Plotting Examples +================= + +The examples below show how wrf-python can be used to make plots with +matplotlib (with basemap and cartopy) and PyNGL. None of these examples +make use of xarray's builtin plotting functions, since additional work is most +likely needed to extend xarray in order to work correctly. This is planned +for a future release. + +Matplotlib With Cartopy +------------------------- + +Cartopy is becoming the main tool for base mapping with matplotlib, but you should +be aware of a few shortcomings when working with WRF data. + +- The builtin tranformations of coordinates when calling the contouring functions + do not work correctly with the rotated pole projection. The + transform_points method needs to be called manually on the latitude and + longitude arrays. + +- The rotated pole projection requires the x and y limits to be set manually + using set_xlim and set_ylim. + +- You can't place latitude and longitude labels on the axes when using + any projection other than Mercator or LatLon. + + +Plotting a Two-dimensional Field +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. image:: _static/images/cartopy_slp.png + :scale: 100% + :align: center + +.. code-block:: python + + from __future__ import (absolute_import, division, print_function, unicode_literals) + + from netCDF4 import Dataset + import matplotlib.pyplot as plt + from matplotlib.cm import get_cmap + import cartopy.crs as crs + from cartopy.feature import NaturalEarthFeature + + from wrf import npvalues, getvar, smooth2d + + ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") + + # Get the sea level pressure + slp = getvar(ncfile, "slp") + + # Smooth the sea level pressure since it tends to be noisy near the mountains + smooth_slp = smooth2d(slp, 3) + + # Get the numpy array from the XLAT and XLONG coordinates + lats = npvalues(slp.coords["XLAT"]) + lons = npvalues(slp.coords["XLONG"]) + + # Get the wrf.WrfProj object + wrf_proj = slp.attrs["projection"] + + # The WrfProj.cartopy() method returns a cartopy.crs projection object + cart_proj = wrf_proj.cartopy() + + # Create a figure that's 10x10 + fig = plt.figure(figsize=(10,10)) + + # Get the GeoAxes set to the projection used by WRF + ax = plt.axes(projection=cart_proj) + + # Download and add the states and coastlines + states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none', + name='admin_1_states_provinces_shp') + ax.add_feature(states, linewidth=.5) + ax.coastlines('50m', linewidth=0.8) + + # Make the contour outlines and filled contours for the smoothed sea level pressure. + # The transform keyword indicates that the lats and lons arrays are lat/lon coordinates and tells + # cartopy to transform the points in to the WRF projection set for the GeoAxes. + plt.contour(lons, lats, npvalues(smooth_slp), 10, colors="black", transform=crs.PlateCarree()) + plt.contourf(lons, lats, npvalues(smooth_slp), 10, transform=crs.PlateCarree()) + + # Add a color bar + plt.colorbar(ax=ax, shrink=.47) + + # Set the map limits + ax.set_xlim(wrf_proj.cartopy_xlim()) + ax.set_ylim(wrf_proj.cartopy_ylim()) + + # Add the gridlines + ax.gridlines() diff --git a/doc/source/user_api/index.rst b/doc/source/user_api/index.rst index d5a1ea3..7e30cbe 100644 --- a/doc/source/user_api/index.rst +++ b/doc/source/user_api/index.rst @@ -1,8 +1,6 @@ User API ============= ------------------ - Routines ------------------ diff --git a/src/wrf/api.py b/src/wrf/api.py index 3ab32b4..6ef5ec9 100644 --- a/src/wrf/api.py +++ b/src/wrf/api.py @@ -29,6 +29,7 @@ from .projection import (WrfProj, NullProjection, LambertConformal, Mercator, PolarStereographic, LatLon, RotatedLatLon, getproj) from .coordpair import CoordPair +from .interputils import to_xy_coords from cache import cache_item, get_cached_item from .version import __version__ @@ -63,6 +64,7 @@ __all__ += ["npvalues", "extract_global_attrs", "is_standard_wrf_var", __all__ += ["WrfProj", "NullProjection", "LambertConformal", "Mercator", "PolarStereographic", "LatLon", "RotatedLatLon", "getproj"] __all__ += ["CoordPair"] +__all__ += ["to_xy_coords"] __all__ += ["cache_item", "get_cached_item"] __all__ += ["__version__"] diff --git a/src/wrf/interputils.py b/src/wrf/interputils.py index 0faa989..0da2c66 100644 --- a/src/wrf/interputils.py +++ b/src/wrf/interputils.py @@ -408,24 +408,26 @@ def to_xy_coords(pairs, wrfin=None, timeidx=0, stagger=None, projection=None): latinc = 0.0 loninc = 0.0 - xy_vals = _ll_to_xy(lat, lon, meta=False, squeeze=True, - as_int=True, - map_proj=projection.map_proj, - truelat1=projection.truelat1, - truelat2=projection.truelat2, - stand_lon=projection.stand_lon, - ref_lat=projection.ll_lat, - ref_lon=projection.ll_lon, - pole_lat=pole_lat, - pole_lon=pole_lon, - known_x=0, - known_y=0, - dx=projection.dx, - dy=projection.dy, - latinc=latinc, - loninc=loninc) - + as_int=True, + map_proj=projection.map_proj, + truelat1=projection.truelat1, + truelat2=projection.truelat2, + stand_lon=projection.stand_lon, + ref_lat=projection.ll_lat, + ref_lon=projection.ll_lon, + pole_lat=pole_lat, + pole_lon=pole_lon, + known_x=0, + known_y=0, + dx=projection.dx, + dy=projection.dy, + latinc=latinc, + loninc=loninc) + + xy_vals = xy_vals.squeeze() + + if xy_vals.ndim == 1: return CoordPair(x=xy_vals[0], y=xy_vals[1]) else: diff --git a/src/wrf/latlonutils.py b/src/wrf/latlonutils.py index 00aee6b..b21bbd6 100644 --- a/src/wrf/latlonutils.py +++ b/src/wrf/latlonutils.py @@ -9,7 +9,8 @@ from .constants import Constants, ProjectionTypes from .extension import _lltoxy, _xytoll from .util import (extract_vars, extract_global_attrs, either, is_moving_domain, is_multi_time_req, - iter_left_indexes, is_mapping, is_multi_file) + iter_left_indexes, is_mapping, is_multi_file, + npvalues) from .py3compat import viewkeys, viewitems from .projutils import dict_keys_to_upper @@ -372,11 +373,10 @@ def _ll_to_xy(latitude, longitude, wrfin=None, timeidx=0, lats = np.asarray(latitude) lons = np.asarray(longitude) - if lats.ndim > 1: - lats = lats.ravel() + # Note: For scalars, this will make a single element array + lats = lats.ravel() - if lons.ndim > 1: - lons = lons.ravel() + lons = lons.ravel() if (lats.size != lons.size): raise ValueError("'latitude' and 'longitude' " @@ -541,11 +541,9 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None, x_arr = x_arr + 1 y_arr = y_arr + 1 - if x_arr.ndim > 1: - x_arr = x_arr.ravel() + x_arr = x_arr.ravel() - if y_arr.ndim > 1: - y_arr = y_arr.ravel() + y_arr = y_arr.ravel() if (x_arr.size != y_arr.size): raise ValueError("'x' and 'y' must be the same length") diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index 0eb8e3a..3e8cce1 100644 --- a/src/wrf/metadecorators.py +++ b/src/wrf/metadecorators.py @@ -783,7 +783,7 @@ def _set_horiz_meta(wrapped, instance, args, kwargs): outname = "field3d_{0}".format(name_levelstr) outattrs = OrderedDict() - outattrs["PlotLevelID"] = levelstr + outattrs["level"] = levelstr outattrs["missing_value"] = missingval outattrs["_FillValue"] = missingval diff --git a/src/wrf/projection.py b/src/wrf/projection.py index d706181..99b6e02 100644 --- a/src/wrf/projection.py +++ b/src/wrf/projection.py @@ -47,7 +47,7 @@ if cartopy_enabled(): extent of the projection. Default is -80.0. max_latitude (:obj:`float`, optional): The maximum northerly - extent of the projection. Default is 84.0.. + extent of the projection. Default is 84.0. globe (:class:`cartopy.crs.Globe`, optional): A globe object. If omitted, a default globe is created. @@ -280,10 +280,13 @@ class WrfProj(object): def __repr__(self): args = ("bottom_left={}, top_right={}, " "stand_lon={}, moad_cen_lat={}, " + "truelat1={}, truelat2={}, " "pole_lat={}, pole_lon={}".format((self.ll_lat, self.ll_lon), (self.ur_lat, self.ur_lon), self.stand_lon, self.moad_cen_lat, + self.truelat1, + self.truelat2, self.pole_lat, self.pole_lon)) return "{}({})".format(self.__class__.__name__, args) diff --git a/src/wrf/util.py b/src/wrf/util.py index b4a1333..061a76a 100644 --- a/src/wrf/util.py +++ b/src/wrf/util.py @@ -920,14 +920,11 @@ def _combine_dict(wrfdict, varname, timeidx, method, meta, _key): is_moving = is_moving_domain(wrfdict, varname, _key=_key) - # Not quite sure how to handle coord caching with dictionaries, so - # disabling it for now by setting _key to 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]) - # Create the output data numpy array based on the first array outdims = [numkeys] outdims += first_array.shape @@ -1587,7 +1584,6 @@ def _cat_files(wrfseq, varname, timeidx, is_moving, squeeze, meta, _key): outdata[startidx:endidx, :] = vardata[:] if xarray_enabled() and meta: - # XTIME new in 3.7 if timename is not None and not timecached: xtimedata = wrfnc.variables[timename][:] outxtimes[startidx:endidx] = xtimedata[:] @@ -2414,7 +2410,7 @@ def extract_times(wrfin, timeidx, method="cat", squeeze=True, cache=None, dims=outdimnames, attrs=outattrs) else: - outarr = np.asarray(time_list) + outarr = np.asarray(time_list, dtype="datetime64[ns]") if not multitime: return outarr[timeidx] diff --git a/src/wrf/uvmet.py b/src/wrf/uvmet.py index 3783bc4..5ac3001 100755 --- a/src/wrf/uvmet.py +++ b/src/wrf/uvmet.py @@ -18,7 +18,7 @@ from .util import extract_vars, extract_global_attrs, either @convert_units("wind", "m s-1") def _get_uvmet(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, meta=True, _key=None, - ten_m=False, units ="m s-1"): + ten_m=False, units="m s-1"): """Return the u,v wind components rotated to earth coordinates. The leftmost dimension of the returned array represents two different @@ -426,7 +426,7 @@ def get_uvmet_wspd_wdir(wrfin, timeidx=0, method="cat", squeeze=True, """ uvmet = _get_uvmet(wrfin, timeidx, method, squeeze, - cache, meta, _key, False, units) + cache, meta, _key, False, units="m s-1") return _calc_wspd_wdir(uvmet[0,...,:,:,:], uvmet[1,...,:,:,:], False, units) @@ -509,7 +509,7 @@ def get_uvmet10_wspd_wdir(wrfin, timeidx=0, method="cat", squeeze=True, """ uvmet10 = _get_uvmet(wrfin, timeidx, method, squeeze, cache, meta, _key, - True, units) + True, units="m s-1") return _calc_wspd_wdir(uvmet10[0,...,:,:], uvmet10[1,...,:,:], True, units) diff --git a/test/ipynb/Doc_Examples.ipynb b/test/ipynb/Doc_Examples.ipynb new file mode 100644 index 0000000..1db43ca --- /dev/null +++ b/test/ipynb/Doc_Examples.ipynb @@ -0,0 +1,415 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "from __future__ import (absolute_import, division, print_function, unicode_literals)\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.cm import get_cmap\n", + "import cartopy.crs as crs\n", + "import cartopy.feature as cfeature\n", + "from netCDF4 import Dataset as nc\n", + "\n", + "from wrf import npvalues, getvar, smooth2d\n", + "\n", + "# Open the output NetCDF with netcdf-python\n", + "filename = \"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00\"\n", + "ncfile = nc(filename)\n", + "\n", + "# Get sea level pressure and cloud top temperature\n", + "slp = getvar(ncfile, \"slp\")\n", + "ctt = getvar(ncfile, \"ctt\")\n", + "\n", + "# Smooth the SLP\n", + "smooth_slp = smooth2d(slp, 3)\n", + "\n", + "# Extract the latitude and longitude coordinate arrays as regular numpy array instead of xarray.DataArray\n", + "lons = npvalues(slp.coords[\"XLONG\"])\n", + "lats = npvalues(slp.coords[\"XLAT\"])\n", + "\n", + "# Get the cartopy projection class\n", + "wrf_proj = slp.attrs[\"projection\"]\n", + "cart_proj = wrf_proj.cartopy()\n", + "\n", + "# Create the figure\n", + "fig = plt.figure(figsize=(4,4))\n", + "ax = plt.axes([0.1,0.1,0.8,0.8], projection=cart_proj)\n", + "\n", + "# Download and create the states, land, and oceans using cartopy features\n", + "states = cfeature.NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n", + " name='admin_1_states_provinces_shp')\n", + "land = cfeature.NaturalEarthFeature(category='physical', name='land', scale='50m', \n", + " facecolor=cfeature.COLORS['land'])\n", + "ocean = cfeature.NaturalEarthFeature(category='physical', name='ocean', scale='50m', \n", + " facecolor=cfeature.COLORS['water'])\n", + "\n", + "# Make the pressure contours.\n", + "contour_levels = [960, 965, 970, 975, 980, 990]\n", + "c1 = plt.contour(lons, lats, npvalues(smooth_slp), levels=contour_levels, colors=\"white\", \n", + " transform=crs.PlateCarree(), zorder=3, linewidths=1.0)\n", + "\n", + "# Add pressure contour labels\n", + "#plt.clabel(c1, contour_levels, inline=True, fmt='%.0f', fontsize=7)\n", + "\n", + "# Create the filled cloud top temperature contours\n", + "contour_levels = [-80, -70, -60, -50, -40, -30, -20, -10, 0, 10]\n", + "plt.contourf(lons, lats, npvalues(ctt), contour_levels, cmap=get_cmap(\"Greys\"),\n", + " transform=crs.PlateCarree(), zorder=2)\n", + "\n", + "plt.plot([-80,-77.8], [26.75,26.75], color=\"yellow\", marker=\"o\", transform=crs.PlateCarree(), zorder=3)\n", + "\n", + "# Create the label bar for cloud top temperature\n", + "#cb2 = plt.colorbar(ax=ax, fraction=0.046, pad=0.04)\n", + "cb2 = plt.colorbar(ax=ax)\n", + "\n", + "# Draw the oceans, land, and states\n", + "ax.add_feature(ocean)\n", + "ax.add_feature(land)\n", + "ax.add_feature(states, linewidth=.5, edgecolor=\"black\")\n", + "\n", + "# Crop the domain to the region around the hurricane\n", + "ax.set_extent([-85., -75.0, np.amin(lats), 30.0], crs=crs.PlateCarree())\n", + "ax.gridlines(crs=crs.PlateCarree(), draw_labels=False)\n", + "\n", + "# Add the title and show the image\n", + "#plt.title(\"Hurricane Matthew Cloud Top Temperature (degC) \")\n", + "plt.savefig(\"/Users/ladwig/Documents/workspace/wrf_python/doc/source/_static/images/matthew.png\",\n", + " transparent=True, bbox_inches=\"tight\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "from __future__ import (absolute_import, division, print_function, unicode_literals)\n", + "\n", + "import numpy as np\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", + "from Nio import open_file\n", + "\n", + "from wrf import npvalues, getvar, smooth2d, ll_to_xy, CoordPair, vertcross, getproj, get_proj_params, to_xy_coords\n", + "\n", + "# Open the output NetCDF file with PyNIO\n", + "filename = \"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00\"\n", + "pynio_filename = filename + \".nc\"\n", + "ncfile = open_file(pynio_filename)\n", + "\n", + "# Extract pressure and model height\n", + "z = getvar(ncfile, \"z\", timeidx=0)\n", + "dbz = getvar(ncfile, \"dbz\", timeidx=0)\n", + "\n", + "wspd = getvar(ncfile, \"uvmet_wspd_wdir\", units=\"kt\")[0,:]\n", + "Z = 10**(dbz/10.)\n", + "\n", + "start_point = CoordPair(lat=26.75, lon=-80.0)\n", + "end_point = CoordPair(lat=26.75, lon=-77.8)\n", + "\n", + "# Compute the vertical cross-section interpolation. Also, include the lat/lon points along the cross-section.\n", + "z_cross = vertcross(Z, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)\n", + "wspd_cross = vertcross(wspd, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)\n", + "dbz_cross = 10.0 * np.log10(z_cross)\n", + "\n", + "# Create the figure\n", + "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(4,4))\n", + "#ax = plt.axes([0.1,0.1,0.8,0.8])\n", + "\n", + "# Define the contour levels [0, 50, 100, 150, ....]\n", + "levels = [5 + 5*n for n in range(15)]\n", + "\n", + "# Make the contour plot\n", + "a = axes[0].contourf(npvalues(wspd_cross))\n", + "# Add the color bar\n", + "fig.colorbar(a, ax=axes[0])\n", + "\n", + "b = axes[1].contourf(npvalues(dbz_cross), levels=levels)\n", + "fig.colorbar(b, ax=axes[1])\n", + "\n", + "# Set the x-ticks to use latitude and longitude labels.\n", + "coord_pairs = npvalues(dbz_cross.coords[\"xy_loc\"])\n", + "x_ticks = np.arange(coord_pairs.shape[0])\n", + "x_labels = [pair.latlon_str() for pair in npvalues(coord_pairs)]\n", + "axes[0].set_xticks(x_ticks[::20])\n", + "axes[0].set_xticklabels([], rotation=45)\n", + "axes[1].set_xticks(x_ticks[::20])\n", + "axes[1].set_xticklabels(x_labels[::20], rotation=45, fontsize=6) \n", + "\n", + "\n", + "# Set the y-ticks to be height.\n", + "vert_vals = npvalues(dbz_cross.coords[\"vertical\"])\n", + "v_ticks = np.arange(vert_vals.shape[0])\n", + "axes[0].set_yticks(v_ticks[::20])\n", + "axes[0].set_yticklabels(vert_vals[::20], fontsize=6) \n", + "axes[1].set_yticks(v_ticks[::20])\n", + "axes[1].set_yticklabels(vert_vals[::20], fontsize=6) \n", + "\n", + "# Set the x-axis and y-axis labels\n", + "axes[1].set_xlabel(\"Latitude, Longitude\", fontsize=7)\n", + "axes[0].set_ylabel(\"Height (m)\", fontsize=7)\n", + "axes[1].set_ylabel(\"Height (m)\", fontsize=7)\n", + "\n", + "# Add a title\n", + "axes[0].set_title(\"Cross-Section of Wind Speed (kt)\", {\"fontsize\" : 10})\n", + "axes[1].set_title(\"Cross-Section of Reflectivity (dBZ)\", {\"fontsize\" : 10})\n", + "\n", + "plt.savefig(\"/Users/ladwig/Documents/workspace/wrf_python/doc/source/_static/images/matthew_cross.png\",\n", + " transparent=True, bbox_inches=\"tight\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "from __future__ import (absolute_import, division, print_function, unicode_literals)\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.cm import get_cmap\n", + "import cartopy.crs as crs\n", + "import cartopy.feature as cfeature\n", + "from Nio import open_file\n", + "\n", + "from wrf import getvar, npvalues, vertcross, smooth2d, CoordPair\n", + "\n", + "# Open the output NetCDF file with PyNIO\n", + "filename = \"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00\"\n", + "pynio_filename = filename + \".nc\"\n", + "ncfile = open_file(pynio_filename)\n", + "\n", + "# Get the WRF variables\n", + "slp = getvar(ncfile, \"slp\")\n", + "smooth_slp = smooth2d(slp, 3)\n", + "ctt = getvar(ncfile, \"ctt\")\n", + "z = getvar(ncfile, \"z\", timeidx=0)\n", + "dbz = getvar(ncfile, \"dbz\", timeidx=0)\n", + "Z = 10**(dbz/10.)\n", + "wspd = getvar(ncfile, \"uvmet_wspd_wdir\", units=\"kt\")[0,:]\n", + "\n", + "# Set the start point and end point for the cross section\n", + "start_point = CoordPair(lat=26.75, lon=-80.0)\n", + "end_point = CoordPair(lat=26.75, lon=-77.8)\n", + "\n", + "# Compute the vertical cross-section interpolation. Also, include the lat/lon points along the cross-section.\n", + "z_cross = vertcross(Z, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)\n", + "wspd_cross = vertcross(wspd, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True, meta=True)\n", + "dbz_cross = 10.0 * np.log10(z_cross)\n", + "\n", + "# Extract the latitude and longitude coordinate arrays as regular numpy array instead of xarray.DataArray\n", + "lons = npvalues(slp.coords[\"XLONG\"])\n", + "lats = npvalues(slp.coords[\"XLAT\"])\n", + "\n", + "# Get the cartopy projection class\n", + "wrf_proj = slp.attrs[\"projection\"]\n", + "cart_proj = wrf_proj.cartopy()\n", + "\n", + "# Create the figure which will have 3 subplots\n", + "fig = plt.figure(figsize=(7,5))\n", + "ax_ctt = fig.add_subplot(1,2,1,projection=cart_proj)\n", + "ax_wspd = fig.add_subplot(2,2,2)\n", + "ax_dbz = fig.add_subplot(2,2,4)\n", + "\n", + "## Plot the cloud top temperature\n", + "\n", + "# Download and create the states, land, and oceans using cartopy features\n", + "states = cfeature.NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n", + " name='admin_1_states_provinces_shp')\n", + "land = cfeature.NaturalEarthFeature(category='physical', name='land', scale='50m', \n", + " facecolor=cfeature.COLORS['land'])\n", + "ocean = cfeature.NaturalEarthFeature(category='physical', name='ocean', scale='50m', \n", + " facecolor=cfeature.COLORS['water'])\n", + "\n", + "# Make the pressure contours.\n", + "contour_levels = [960, 965, 970, 975, 980, 990]\n", + "c1 = ax_ctt.contour(lons, lats, npvalues(smooth_slp), levels=contour_levels, colors=\"white\", \n", + " transform=crs.PlateCarree(), zorder=3, linewidths=1.0)\n", + "\n", + "# Create the filled cloud top temperature contours\n", + "contour_levels = [-80.0, -70.0, -60, -50, -40, -30, -20, -10, 0, 10]\n", + "ctt_contours = ax_ctt.contourf(lons, lats, npvalues(ctt), contour_levels, cmap=get_cmap(\"Greys\"),\n", + " transform=crs.PlateCarree(), zorder=2)\n", + "\n", + "ax_ctt.plot([start_point.lon, end_point.lon], [start_point.lat, end_point.lat], color=\"yellow\", \n", + " marker=\"o\", transform=crs.PlateCarree(), zorder=3)\n", + "\n", + "# Create the label bar for cloud top temperature\n", + "cb_ctt = fig.colorbar(ctt_contours, ax=ax_ctt, shrink=.5)\n", + "cb_ctt.ax.tick_params(labelsize=5)\n", + "\n", + "# Draw the oceans, land, and states\n", + "ax_ctt.add_feature(ocean)\n", + "ax_ctt.add_feature(land)\n", + "ax_ctt.add_feature(states, linewidth=.5, edgecolor=\"black\")\n", + "\n", + "# Crop the domain to the region around the hurricane\n", + "ax_ctt.set_extent([-85., -75.0, np.amin(lats), 30.0], crs=crs.PlateCarree())\n", + "ax_ctt.gridlines(crs=crs.PlateCarree(), draw_labels=False)\n", + "\n", + "## Plot the cross sections\n", + "\n", + "# Make the contour plot for wspd\n", + "wspd_contours = ax_wspd.contourf(npvalues(wspd_cross))\n", + "# Add the color bar\n", + "cb_wspd = fig.colorbar(wspd_contours, ax=ax_wspd)\n", + "cb_wspd.ax.tick_params(labelsize=5)\n", + "\n", + "# Make the contour plot for dbz\n", + "levels = [5 + 5*n for n in range(15)]\n", + "dbz_contours = ax_dbz.contourf(npvalues(dbz_cross), levels=levels)\n", + "cb_dbz = fig.colorbar(dbz_contours, ax=ax_dbz)\n", + "cb_dbz.ax.tick_params(labelsize=5)\n", + "\n", + "# Set the x-ticks to use latitude and longitude labels.\n", + "coord_pairs = npvalues(dbz_cross.coords[\"xy_loc\"])\n", + "x_ticks = np.arange(coord_pairs.shape[0])\n", + "x_labels = [pair.latlon_str() for pair in npvalues(coord_pairs)]\n", + "ax_wspd.set_xticks(x_ticks[::20])\n", + "ax_wspd.set_xticklabels([], rotation=45)\n", + "ax_dbz.set_xticks(x_ticks[::20])\n", + "ax_dbz.set_xticklabels(x_labels[::20], rotation=45, fontsize=4) \n", + "\n", + "\n", + "# Set the y-ticks to be height.\n", + "vert_vals = npvalues(dbz_cross.coords[\"vertical\"])\n", + "v_ticks = np.arange(vert_vals.shape[0])\n", + "ax_wspd.set_yticks(v_ticks[::20])\n", + "ax_wspd.set_yticklabels(vert_vals[::20], fontsize=4) \n", + "ax_dbz.set_yticks(v_ticks[::20])\n", + "ax_dbz.set_yticklabels(vert_vals[::20], fontsize=4) \n", + "\n", + "# Set the x-axis and y-axis labels\n", + "ax_dbz.set_xlabel(\"Latitude, Longitude\", fontsize=5)\n", + "ax_wspd.set_ylabel(\"Height (m)\", fontsize=5)\n", + "ax_dbz.set_ylabel(\"Height (m)\", fontsize=5)\n", + "\n", + "# Add a title\n", + "ax_ctt.set_title(\"Cloud Top Temperature (degC)\", {\"fontsize\" : 7})\n", + "ax_wspd.set_title(\"Cross-Section of Wind Speed (kt)\", {\"fontsize\" : 7})\n", + "ax_dbz.set_title(\"Cross-Section of Reflectivity (dBZ)\", {\"fontsize\" : 7})\n", + "\n", + "plt.savefig(\"/Users/ladwig/Documents/workspace/wrf_python/doc/source/_static/images/matthew.png\",\n", + " transparent=True, bbox_inches=\"tight\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# SLP\n", + "from __future__ import (absolute_import, division, print_function, unicode_literals)\n", + " \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 npvalues, getvar, smooth2d\n", + "\n", + "ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00\")\n", + "\n", + "# Get the sea level pressure\n", + "slp = getvar(ncfile, \"slp\")\n", + "\n", + "# Smooth the sea level pressure since it tends to be noisey near the mountains\n", + "smooth_slp = smooth2d(slp, 3)\n", + "\n", + "# Get the numpy array from the XLAT and XLONG coordinates\n", + "lats = npvalues(slp.coords[\"XLAT\"])\n", + "lons = npvalues(slp.coords[\"XLONG\"])\n", + "\n", + "# Get the wrf.WrfProj object\n", + "wrf_proj = slp.attrs[\"projection\"]\n", + "\n", + "# The cartopy() method returns a cartopy.crs projection object\n", + "cart_proj = wrf_proj.cartopy()\n", + "\n", + "# Create a figure that's 10x10\n", + "fig = plt.figure(figsize=(10,10))\n", + "# Get the GeoAxes set 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", + "# The transform keyword indicates that the lats and lons arrays are lat/lon coordinates and tells \n", + "# cartopy to transform the points in to grid space.\n", + "plt.contour(lons, lats, npvalues(smooth_slp), 10, colors=\"black\", transform=crs.PlateCarree())\n", + "plt.contourf(lons, lats, npvalues(smooth_slp), 10, transform=crs.PlateCarree())\n", + "\n", + "# Add a color bar\n", + "plt.colorbar(ax=ax, shrink=.47)\n", + "\n", + "# Set the map limits\n", + "ax.set_xlim(wrf_proj.cartopy_xlim())\n", + "ax.set_ylim(wrf_proj.cartopy_ylim())\n", + "\n", + "# Add the gridlines\n", + "ax.gridlines()" + ] + } + ], + "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.12" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/test/ipynb/WRF_Workshop_Demo.ipynb b/test/ipynb/WRF_Workshop_Demo.ipynb index a2f0a3e..601b1c6 100644 --- a/test/ipynb/WRF_Workshop_Demo.ipynb +++ b/test/ipynb/WRF_Workshop_Demo.ipynb @@ -69,7 +69,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": { "collapsed": true }, @@ -81,11 +81,24 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": { "collapsed": false }, - "outputs": [], + "outputs": [ + { + "ename": "NIOError", + "evalue": "Unable to open file", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mNIOError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0mfilename\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"wrfout_d01_2010-06-13_21-00-00\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0mpynio_filename\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfilename\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m\".nc\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0mncfile\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen_file\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpynio_filename\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 11\u001b[0m \u001b[0;31m# Alternative using netCDF4-python (for reference)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/Users/ladwig/miniconda2/lib/python2.7/site-packages/PyNIO/Nio.pyc\u001b[0m in \u001b[0;36mopen_file\u001b[0;34m(filename, mode, options, history, format)\u001b[0m\n\u001b[1;32m 733\u001b[0m \u001b[0mmask_above_value\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_get_option_value\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0m_Nio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moption_defaults\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'MaskAboveValue'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 734\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 735\u001b[0;31m \u001b[0mfile\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_Nio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen_file\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mhistory\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 736\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 737\u001b[0m \u001b[0mfile_proxy\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_proxy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfile\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'str'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m__del__\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0m__del__\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcreate_variable\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcreate_variable\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcreate_group\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcreate_group\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mclose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mclose\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mNIOError\u001b[0m: Unable to open file" + ] + } + ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", @@ -424,21 +437,21 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 2", "language": "python", - "name": "python3" + "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 3 + "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.5.2" + "pygments_lexer": "ipython2", + "version": "2.7.12" } }, "nbformat": 4, diff --git a/test/ipynb/WRF_python_demo.ipynb b/test/ipynb/WRF_python_demo.ipynb index b647315..520ea13 100644 --- a/test/ipynb/WRF_python_demo.ipynb +++ b/test/ipynb/WRF_python_demo.ipynb @@ -19,14 +19,16 @@ "\n", "from wrf import getvar\n", "from netCDF4 import Dataset as nc\n", - "ncfile = nc(\"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-02-25_18_00_00\")\n" + "#ncfile = nc(\"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-02-25_18_00_00\")\n", + "ncfile = nc(\"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00\")\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false + "collapsed": false, + "scrolled": false }, "outputs": [], "source": [ @@ -148,8 +150,9 @@ }, "outputs": [], "source": [ + "from wrf import ALL_TIMES\n", "wrflist = [ncfile, ncfile, ncfile]\n", - "p_cat = getvar(wrflist, \"P\", method=\"cat\")\n", + "p_cat = getvar(wrflist, \"P\", timeidx=ALL_TIMES, method=\"cat\")\n", "print(p_cat)\n", "del p_cat" ] @@ -169,7 +172,7 @@ }, "outputs": [], "source": [ - "p_join = getvar(wrflist, \"P\", method=\"join\")\n", + "p_join = getvar(wrflist, \"P\", timeidx=ALL_TIMES, method=\"join\")\n", "print(p_join)" ] }, @@ -191,7 +194,8 @@ }, "outputs": [], "source": [ - "p_join = getvar(wrflist, \"P\", timeidx=0, method=\"join\", squeeze=False)\n", + "from wrf import ALL_TIMES\n", + "p_join = getvar(wrflist, \"P\", timeidx=ALL_TIMES, method=\"join\", squeeze=False)\n", "print(p_join)\n", "del p_join" ] @@ -213,7 +217,7 @@ "source": [ "wrf_dict = {\"label1\" : [ncfile, ncfile],\n", " \"label2\" : [ncfile, ncfile]}\n", - "p_dict = getvar(wrf_dict, \"P\")\n", + "p_dict = getvar(wrf_dict, \"P\", timeidx=ALL_TIMES)\n", "print(p_dict)\n", "del p_dict" ] @@ -308,8 +312,9 @@ " \"pvo\", \"pw\", \"rh2\", \"rh\", \"slp\", \"ter\", \"td2\", \"td\", \"tc\",\n", " \"theta\", \"tk\", \"tv\", \"twb\", \"updraft_helicity\", \"ua\", \"va\", \n", " \"wa\", \"uvmet10\", \"uvmet\", \"z\", \"ctt\"]\n", + "wrf_vars = [\"slp\"]\n", "\n", - "vard = {varname: getvar(ncfile, varname, method=\"join\", squeeze=False) for varname in wrf_vars}\n", + "vard = {varname: getvar(ncfile, varname, method=\"cat\", squeeze=True) for varname in wrf_vars}\n", "for varname in wrf_vars:\n", " print(vard[varname])\n", " print (\"\\n\")\n" @@ -331,8 +336,8 @@ "outputs": [], "source": [ "from wrf import npvalues\n", - "masked_ndarray = npvalues(vard[\"cape_2d\"])\n", - "print(masked_ndarray)\n", + "masked_ndarray = npvalues(vard[\"slp\"])\n", + "print(type(masked_ndarray))\n", "del masked_ndarray" ] }, @@ -398,11 +403,11 @@ "outputs": [], "source": [ "# Pressure using pivot and angle\n", - "from wrf import getvar, vertcross\n", + "from wrf import getvar, vertcross, CoordPair\n", "\n", "z = getvar(ncfile, \"z\")\n", "p = getvar(ncfile, \"pressure\")\n", - "pivot_point = (z.shape[-1] / 2, z.shape[-2] / 2) \n", + "pivot_point = CoordPair((z.shape[-1]-1) // 2, (z.shape[-2] - 1) // 2) \n", "angle = 90.0\n", "\n", "p_vert = vertcross(p, z, pivot_point=pivot_point, angle=angle, latlon=True)\n", @@ -411,14 +416,122 @@ "del p_vert\n", "\n", "# Pressure using start_point and end_point\n", - "start_point = (0, z.shape[-2]/2)\n", - "end_point = (-1, z.shape[-2]/2)\n", + "start_point = CoordPair(0, (z.shape[-2]-1) // 2)\n", + "end_point = CoordPair(-1, (z.shape[-2]-1) // 2)\n", "\n", "p_vert = vertcross(p, z, start_point=start_point, end_point=end_point, latlon=True)\n", "print(p_vert)\n", "del p_vert, p, z" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Pressure using pivot and angle\n", + "from wrf import getvar, vertcross, CoordPair, xy_to_ll\n", + "\n", + "z = getvar(ncfile, \"z\")\n", + "p = getvar(ncfile, \"pressure\")\n", + "lats = getvar(ncfile, \"lat\")\n", + "lons = getvar(ncfile, \"lon\")\n", + "\n", + "#print ((lats.shape[-2]-1) / 2)\n", + "#print ((lats.shape[-1]-1) / 2)\n", + "\n", + "#print (npvalues(lats[529, 899]))\n", + "#print (npvalues(lons[529, 899]))\n", + "\n", + "#print (npvalues(lats[529, 0]))\n", + "#print (npvalues(lons[529, 0]))\n", + "\n", + "#print (npvalues(lats[529, -1]))\n", + "#print (npvalues(lons[529, -1]))\n", + "\n", + "pivot_point = CoordPair(lat=38.5, lon=-97.5) \n", + "angle = 90.0\n", + "\n", + "p_vert = vertcross(p, z, wrfin=ncfile, pivot_point=pivot_point, angle=angle, latlon=True)\n", + "print (p_vert)\n", + "print (\"\\n\")\n", + "\n", + "start_lat = lats[(lats.shape[-2]-1)//2, 0]\n", + "end_lat = lats[(lats.shape[-2]-1)//2, -1]\n", + "start_lon = lons[(lats.shape[-2]-1)//2, 0]\n", + "end_lon = lons[(lats.shape[-2]-1)//2, -1]\n", + "\n", + "print (start_lat)\n", + "print (end_lat)\n", + "print (start_lon)\n", + "print (end_lon)\n", + "\n", + "# Pressure using start_point and end_point\n", + "start_point = CoordPair(lat=start_lat, lon=start_lon)\n", + "end_point = CoordPair(lat=end_lat, lon=end_lon)\n", + "\n", + "p_vert = vertcross(p, z, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True)\n", + "print(p_vert)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Pressure using pivot and angle\n", + "from wrf import getvar, vertcross, CoordPair, xy_to_ll\n", + "\n", + "z = getvar(ncfile, \"z\")\n", + "p = getvar(ncfile, \"pressure\")\n", + "lats = getvar(ncfile, \"lat\")\n", + "lons = getvar(ncfile, \"lon\")\n", + "\n", + "#print ((lats.shape[-2]-1) / 2)\n", + "#print ((lats.shape[-1]-1) / 2)\n", + "\n", + "#print (npvalues(lats[529, 899]))\n", + "#print (npvalues(lons[529, 899]))\n", + "\n", + "#print (npvalues(lats[529, 0]))\n", + "#print (npvalues(lons[529, 0]))\n", + "\n", + "#print (npvalues(lats[529, -1]))\n", + "#print (npvalues(lons[529, -1]))\n", + "\n", + "pivot_point = CoordPair(lat=38.5, lon=-97.5) \n", + "angle = 90.0\n", + "\n", + "p_vert = vertcross(p, z, wrfin=ncfile, pivot_point=pivot_point, angle=angle, latlon=True)\n", + "print (p_vert)\n", + "print (\"\\n\")\n", + "\n", + "start_lat = lats[(lats.shape[-2]-1)//2, 0]\n", + "end_lat = lats[(lats.shape[-2]-1)//2, -1]\n", + "start_lon = lons[(lats.shape[-2]-1)//2, 0]\n", + "end_lon = lons[(lats.shape[-2]-1)//2, -1]\n", + "\n", + "print (start_lat)\n", + "print (end_lat)\n", + "print (start_lon)\n", + "print (end_lon)\n", + "\n", + "# Pressure using start_point and end_point\n", + "start_point = CoordPair(lat=start_lat, lon=start_lon)\n", + "end_point = CoordPair(lat=end_lat, lon=end_lon)\n", + "\n", + "levels = [1000., 2000., 3000.]\n", + "\n", + "p_vert = vertcross(p, z, wrfin=ncfile, levels=levels, start_point=start_point, end_point=end_point, latlon=True)\n", + "print(p_vert)" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -435,11 +548,11 @@ "outputs": [], "source": [ "# T2 using pivot and angle\n", - "from wrf import interpline, getvar\n", + "from wrf import interpline, getvar, CoordPair\n", "\n", "t2 = getvar(ncfile, \"T2\")\n", - "pivot_point = (t2.shape[-2] / 2, t2.shape[-1] / 2) \n", - "angle = 90.0\n", + "pivot_point = CoordPair((t2.shape[-1]-1)//2, (t2.shape[-2]-1)//2) \n", + "angle = 0.0\n", "\n", "t2_line = interpline(t2, pivot_point=pivot_point, angle=angle, latlon=True)\n", "print(t2_line, \"\\n\")\n", @@ -447,13 +560,30 @@ "del t2_line\n", "\n", "# T2 using start_point and end_point\n", - "start_point = (t2.shape[-2]/2, 0)\n", - "end_point = (t2.shape[-2]/2, -1)\n", + "start_point = CoordPair((t2.shape[-1]-1)//2, 0)\n", + "end_point = CoordPair((t2.shape[-1]-1)//2, -1)\n", "\n", "t2_line = interpline(t2, start_point=start_point, end_point=end_point, latlon=True)\n", "print(t2_line, \"\\n\")\n", "\n", - "del t2_line, t2" + "del t2_line\n", + "\n", + "t2 = getvar(ncfile, \"T2\")\n", + "lats = getvar(ncfile, \"lat\")\n", + "lons = getvar(ncfile, \"lon\")\n", + "\n", + "start_lat = lats[0, (lats.shape[-1]-1)//2]\n", + "end_lat = lats[-1, (lats.shape[-1]-1)//2]\n", + "start_lon = lons[0, (lons.shape[-1]-1)//2]\n", + "end_lon = lons[-1, (lons.shape[-1]-1)//2]\n", + "\n", + "start_point = CoordPair(lat=start_lat, lon=start_lon)\n", + "end_point = CoordPair(lat=end_lat, lon=end_lon)\n", + "\n", + "t2_line = interpline(t2, wrfin=ncfile, start_point=start_point, end_point=end_point, latlon=True)\n", + "print (t2_line)\n", + "\n", + "del t2_line, t2\n" ] }, { @@ -566,7 +696,16 @@ "from wrf.latlon import xy_to_ll, ll_to_xy \n", "\n", "a = xy_to_ll(ncfile, 400, 200)\n", + "a1 = ll_to_xy(ncfile, a[0], a[1])\n", + "\n", + "#print(a)\n", + "#print(\"\\n\")\n", + "#print(a1)\n", + "#print(\"\\n\")\n", + "\n", "a = xy_to_ll(ncfile, [400,105], [200,205])\n", + "a1 = ll_to_xy(ncfile, a[0,:], a[1,:])\n", + "\n", "b = ll_to_xy(ncfile, 45.5, -110.8, as_int=True)\n", "\n", "# Note: Lists/Dictionaries of files will add a new dimension ('domain') only if the domain is moving\n", @@ -577,6 +716,8 @@ "\n", "print(a)\n", "print(\"\\n\")\n", + "print(a1)\n", + "print(\"\\n\")\n", "print(b)\n", "print(\"\\n\")\n", "print(c)\n", @@ -618,17 +759,18 @@ "import cartopy.crs as crs\n", "from cartopy.feature import NaturalEarthFeature\n", "\n", - "from wrf import npvalues, getvar\n", + "from wrf import npvalues, getvar, smooth2d\n", "\n", "slp = getvar(ncfile, \"slp\")\n", + "smooth_slp = smooth2d(slp, 3)\n", "lons = npvalues(slp.coords[\"XLONG\"])\n", "lats = npvalues(slp.coords[\"XLAT\"])\n", "\n", "wrf_proj = slp.attrs[\"projection\"]\n", "cart_proj = wrf_proj.cartopy()\n", "\n", - "fig = plt.figure(figsize=(20,20))\n", - "ax = plt.axes([0.1,0.1,0.8,0.8], projection=cart_proj)\n", + "fig = plt.figure(figsize=(10,10))\n", + "ax = plt.axes(projection=cart_proj)\n", "\n", "states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n", " name='admin_1_states_provinces_shp')\n", @@ -640,14 +782,80 @@ "x = xform_coords[:,:,0]\n", "y = xform_coords[:,:,1]\n", "\n", - "plt.contour(x, y, npvalues(slp), 20, cmap=get_cmap(\"gist_ncar\"))\n", - "plt.colorbar(ax=ax, shrink=.7)\n", + "plt.contour(x, y, npvalues(smooth_slp), 10, colors=\"black\")\n", + "plt.contourf(x, y, npvalues(smooth_slp), 10)\n", + "plt.colorbar(ax=ax, shrink=.47)\n", "\n", "ax.set_xlim(wrf_proj.cartopy_xlim())\n", "ax.set_ylim(wrf_proj.cartopy_ylim())\n", "ax.gridlines()\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# SLP\n", + "from __future__ import (absolute_import, division, print_function, unicode_literals)\n", + " \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 npvalues, getvar, smooth2d\n", + "\n", + "ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrfout_d01_2016-10-07_00_00_00\")\n", + "\n", + "# Get the sea level pressure\n", + "slp = getvar(ncfile, \"slp\")\n", + "\n", + "# Smooth the sea level pressure since it tends to be noisey near the mountains\n", + "smooth_slp = smooth2d(slp, 3)\n", + "\n", + "# Get the numpy array from the XLAT and XLONG coordinates\n", + "lats = npvalues(slp.coords[\"XLAT\"])\n", + "lons = npvalues(slp.coords[\"XLONG\"])\n", + "\n", + "# Get the wrf.WrfProj object\n", + "wrf_proj = slp.attrs[\"projection\"]\n", + "\n", + "# The cartopy() method returns a cartopy.crs projection object\n", + "cart_proj = wrf_proj.cartopy()\n", + "\n", + "# Create a figure that's 10x10\n", + "fig = plt.figure(figsize=(10,10))\n", + "# Get the GeoAxes set 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", + "# The transform keyword indicates that the lats and lons arrays are lat/lon coordinates and tells \n", + "# cartopy to transform the points in to grid space.\n", + "plt.contour(lons, lats, npvalues(smooth_slp), 10, colors=\"black\", transform=crs.PlateCarree())\n", + "plt.contourf(lons, lats, npvalues(smooth_slp), 10, transform=crs.PlateCarree())\n", + "\n", + "# Add a color bar\n", + "plt.colorbar(ax=ax, shrink=.47)\n", + "\n", + "# Set the map limits\n", + "ax.set_xlim(wrf_proj.cartopy_xlim())\n", + "ax.set_ylim(wrf_proj.cartopy_ylim())\n", + "\n", + "# Add the gridlines\n", + "ax.gridlines()" + ] + }, { "cell_type": "code", "execution_count": null, @@ -708,7 +916,8 @@ "cell_type": "code", "execution_count": null, "metadata": { - "collapsed": false + "collapsed": false, + "scrolled": false }, "outputs": [], "source": [ @@ -717,15 +926,15 @@ "import matplotlib.pyplot as plt\n", "from matplotlib.cm import get_cmap\n", "\n", - "from wrf import getvar, vertcross, npvalues\n", + "from wrf import getvar, vertcross, npvalues, CoordPair\n", "\n", "p = getvar(ncfile, \"pressure\")\n", "z = getvar(ncfile, \"z\", units=\"dm\")\n", "\n", - "pivot_point = (z.shape[-1] / 2, z.shape[-2] / 2) \n", + "pivot_point = CoordPair(z.shape[-1] / 2, z.shape[-2] / 2) \n", "angle = 90.0\n", "\n", - "p_vert = vertcross(p, z, pivot_point=pivot_point, angle=angle)\n", + "p_vert = vertcross(p, z, pivot_point=pivot_point, angle=angle, levels=[1000,850,500])\n", "\n", "fig = plt.figure(figsize=(20,8))\n", "ax = plt.axes([0.1,0.1,0.8,0.8])\n", @@ -1038,21 +1247,21 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 2", "language": "python", - "name": "python3" + "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 3 + "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.5.2" + "pygments_lexer": "ipython2", + "version": "2.7.12" } }, "nbformat": 4, diff --git a/test/utests.py b/test/utests.py index a78940e..ce502bf 100644 --- a/test/utests.py +++ b/test/utests.py @@ -8,7 +8,7 @@ import subprocess from wrf import (getvar, interplevel, interpline, vertcross, vinterp, disable_xarray, xarray_enabled, npvalues, xy_to_ll, ll_to_xy, xy_to_ll_proj, ll_to_xy_proj, - extract_global_attrs, viewitems) + extract_global_attrs, viewitems, CoordPair) from wrf.util import is_multi_file NCL_EXE = "/Users/ladwig/nclbuild/6.3.0/bin/ncl" @@ -258,7 +258,7 @@ def make_interp_test(varname, wrf_in, referent, multi=False, hts = getvar(in_wrfnc, "z", timeidx=timeidx) p = getvar(in_wrfnc, "pressure", timeidx=timeidx) - pivot_point = (hts.shape[-1] / 2, hts.shape[-2] / 2) + pivot_point = CoordPair(hts.shape[-1] / 2, hts.shape[-2] / 2) ht_cross = vertcross(hts, p, pivot_point=pivot_point, angle=90.) nt.assert_allclose(npvalues(ht_cross), ref_ht_cross, rtol=.01) @@ -270,8 +270,8 @@ def make_interp_test(varname, wrf_in, referent, multi=False, ref_p_cross, rtol=.01) # Test point to point - start_point = (0,hts.shape[-2]/2) - end_point = (-1,hts.shape[-2]/2) + start_point = CoordPair(0, hts.shape[-2]/2) + end_point = CoordPair(-1,hts.shape[-2]/2) p_cross2 = vertcross(p,hts,start_point=start_point, end_point=end_point) @@ -284,15 +284,15 @@ def make_interp_test(varname, wrf_in, referent, multi=False, ref_t2_line = _get_refvals(referent, "t2_line", repeat, multi) t2 = getvar(in_wrfnc, "T2", timeidx=timeidx) - pivot_point = (t2.shape[-1] / 2, t2.shape[-2] / 2) + pivot_point = CoordPair(t2.shape[-1] / 2, t2.shape[-2] / 2) t2_line1 = interpline(t2, pivot_point=pivot_point, angle=90.0) nt.assert_allclose(npvalues(t2_line1), ref_t2_line) # Test point to point - start_point = (0, t2.shape[-2]/2) - end_point = (-1, t2.shape[-2]/2) + start_point = CoordPair(0, t2.shape[-2]/2) + end_point = CoordPair(-1, t2.shape[-2]/2) t2_line2 = interpline(t2, start_point=start_point, end_point=end_point)