Browse Source

Merge branch 'release/1.1.3'

main 1.1.3
Bill Ladwig 7 years ago
parent
commit
31091c374e
  1. 6
      doc/source/_templates/product_table.txt
  2. 22
      doc/source/new.rst
  3. 3
      doc/source/tutorial.rst
  4. 4
      doc/source/tutorials/wrf_workshop_2017.rst
  5. 410
      doc/source/tutorials/wrf_workshop_2018.rst
  6. 28
      fortran/rip_cape.f90
  7. 56
      fortran/wrf_fctt.f90
  8. 6
      fortran/wrf_vinterp.f90
  9. 23
      fortran/wrffortran.pyf
  10. 73
      src/wrf/computation.py
  11. 55
      src/wrf/extension.py
  12. 6
      src/wrf/g_cape.py
  13. 38
      src/wrf/g_ctt.py
  14. 8
      src/wrf/g_dewpoint.py
  15. 2
      src/wrf/g_helicity.py
  16. 131
      src/wrf/g_latlon.py
  17. 8
      src/wrf/g_rh.py
  18. 4
      src/wrf/g_slp.py
  19. 22
      src/wrf/interp.py
  20. 27
      src/wrf/latlonutils.py
  21. 40
      src/wrf/metadecorators.py
  22. 5
      src/wrf/routines.py
  23. 45
      src/wrf/specialdec.py
  24. 10
      src/wrf/util.py
  25. 2
      src/wrf/version.py
  26. BIN
      test/ci_tests/ci_result_file.nc
  27. 2
      test/ci_tests/make_test_file.py
  28. 65
      test/comp_utest.py
  29. 50
      test/ctt_test.py
  30. 2
      test/ipynb/Doc_Examples.ipynb
  31. 214
      test/ipynb/Test_CTT.ipynb
  32. 214
      test/ipynb/Test_CTT_Prod.ipynb
  33. 37
      test/ipynb/WRF_python_demo.ipynb
  34. 181
      test/test_inputs.py
  35. 2
      test/utests.py

6
doc/source/_templates/product_table.txt

@ -13,11 +13,13 @@ @@ -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'. |
| | | | |

22
doc/source/new.rst

@ -4,6 +4,28 @@ What's New @@ -4,6 +4,28 @@ 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.
- Dictionary inputs now work correctly with xy_to_ll and ll_to_xy.
- The cape_2d diagnostic can now work with a single column of data, just like
cape_3d.
v1.1.2
^^^^^^^^^^^^^^

3
doc/source/tutorial.rst

@ -12,7 +12,7 @@ Upcoming Tutorials @@ -12,7 +12,7 @@ Upcoming Tutorials
.. toctree::
:maxdepth: 1
tutorials/tutorial_03_2018.rst
tutorials/wrf_workshop_2018.rst
Past Tutorials
@ -22,4 +22,5 @@ Past Tutorials @@ -22,4 +22,5 @@ Past Tutorials
:maxdepth: 1
tutorials/wrf_workshop_2017.rst
tutorials/tutorial_03_2018.rst

4
doc/source/tutorials/wrf_workshop_2017.rst

@ -1,5 +1,5 @@ @@ -1,5 +1,5 @@
WRF Workshop 2017
=====================
WRF Users' Workshop 2017
==========================
Welcome wrf-python tutorial attendees!

410
doc/source/tutorials/wrf_workshop_2018.rst

@ -0,0 +1,410 @@ @@ -0,0 +1,410 @@
WRF Users' Workshop 2018
=========================
Welcome WRF-Python tutorial attendees!
The instructions below should be completed prior to arriving at the tutorial.
There will not be enough time to do this during the tutorial.
Prerequisites
---------------
This tutorial assumes that you have basic knowledge of how to type commands
in to a command terminal using your preferred operating system. You
should know some basic directory commands like *cd*, *mkdir*, *cp*, *mv*.
Regarding Python, to understand the examples in this tutorial, you
should have some experience with Python basics. Below is a list of some
Python concepts that you will see in the examples, but don't worry if you aren't
familiar with everything.
- Opening a Python interpreter and entering commands.
- Importing packages via the import statement.
- Familiarity with some of the basic Python types: str, list, tuple, dict, bool, float, int, None.
- Creating a list, tuple, or dict with "[ ]", "( )", "{ }" syntax (e.g. my_list = [1,2,3,4,5]).
- Accessing dict/list/tuple items with the "x[ ]" syntax (e.g. my_list_item = my_list[0]).
- Slicing str/list/tuple with the ":" syntax (e.g. my_slice = my_list[1:3]).
- Using object methods and attributes with the "x.y" syntax (e.g. my_list.append(6)).
- Calling functions (e.g. result = some_function(x, y))
- Familiarity with numpy would be helpful, as only a very brief introduction
is provided.
- Familiarity with matplotlib would be helpful, as only a very brief
introduction is provided.
If you are completely new to Python, that shouldn't be a problem, since
most of the examples consist of basic container types and function calls. It
would be helpful to look at some introductory material before arriving at the
tutorial. If you've programmed before, picking up Python isn't too difficult.
Here are some links:
https://www.learnpython.org/
https://developers.google.com/edu/python/
Step 1: Open a Command Terminal
--------------------------------
To begin, you will first need to know how to open a command line terminal for
your operating system.
For Windows:
.. code-block:: none
WINDOWS + r
type cmd in the run window
For Mac:
.. code-block:: none
Finder -> Applications -> Utilities -> Terminal
For Linux:
.. code-block:: none
Try one of the following:
CTRL + ALT + T
CTRL + ALT + F2
Step 2: Download Miniconda
----------------------------
For this tutorial, you will need to download and install Miniconda. We are
going to use Python 2.7, but it will also work with Python 3.5+.
Please use the appropriate link below to download Miniconda for your operating
system.
.. note::
64-bit OS recommended
`Win64 <https://repo.continuum.io/miniconda/Miniconda2-latest-Windows-x86_64.exe>`_
`Mac <https://repo.continuum.io/miniconda/Miniconda2-latest-MacOSX-x86_64.sh>`_
`Linux <https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh>`_
For more information, see: https://conda.io/miniconda.html
.. note::
**What is Miniconda?**
If you have used the Anaconda distribution for Python before, then you will be
familiar with Miniconda. The Anaconda Python distribution includes numerous
scientific packages out of the box, which can be difficult for users to build and
install. More importantly, Anaconda includes the conda package manager.
The conda package manager is a utility (similar to yum or apt-get) that installs
packages from a repository of pre-compiled Python packages. These repositories
are called channels. Conda makes it easy for Python users to install and
uninstall packages, and also can be used to create isolated Python environments
(more on that later).
Miniconda is a bare bones implementation of Anaconda and only includes the
conda package manager. Since we are going to use the conda-forge channel to
install our scientific packages, Miniconda avoids any complications between
packages provided by Anaconda and conda-forge.
Step 3: Install Miniconda
----------------------------
Windows:
1. Browse to the directory where you downloaded Miniconda2-latest-Windows-x86_64.exe.
2. Double click on Miniconda2-latest-Windows-x86_64.exe.
3. Follow the instructions.
4. Restart your command terminal.
Mac and Linux:
For Mac and Linux, the installer is a bash script.
1. Using a terminal, you need to execute the bash shell script that you downloaded by
doing::
bash /path/to/Miniconda2-latest-MacOSX-x86_64.sh [Mac]
bash /path/to/Miniconda2-latest-Linux-x86_64.sh [Linux]
2. Follow the instructions.
3. At the end of the installation, it will ask if you want to add the
miniconda2 path to your bash environment. If you are unsure what to do,
you should say "yes". If you say "no", we're going to assume you know
what you are doing.
If you said "yes", then once you restart your shell, the miniconda2 Python
will be found instead of the system Python when you type the "python"
command. If you want to undo this later, then you can edit
either ~/.bash_profile or ~/.bashrc (depending on OS used) and
comment out the line that looks similar to::
# added by Miniconda2 4.1.11 installer
export PATH="/path/to/miniconda2/bin:$PATH"
4. Restart your command terminal.
5. [Linux and Mac Users Only] Miniconda only works with bash. If bash is
not your default shell, then you need to activate the bash shell by typing
the following in to your command terminal::
bash
6. Verify that your system is using the correct Python interpreter by typing
the following in to your command terminal::
which python
You should see the path to your miniconda installation. If not, see the
note below.
.. note::
If you have already installed another Python distribution, like Enthought
Canopy, you will need to comment out any PATH entries for that distribution
in your .bashrc or .bash_profile. Otherwise, your shell environment may
pick to wrong Python installation.
If bash is not your default shell type, and the PATH variable has been
set in .bash_profile by the miniconda installer, try executing
"bash -l" instead of the "bash" command in step 5.
Step 4: Set Up the Conda Environment
--------------------------------------
If you are new to the conda package manager, one of the nice features of conda
is that you can create isolated Python environments that prevent package
incompatibilities. This is similar to the *virtualenv* package that some
Python users may be familiar with. However, conda is not compatible with
virtualenv, so only use conda environments when working with conda.
The name of our conda environment for this tutorial is: **tutorial_2018**.
Follow the instructions below to create the tutorial_2018 environment.
1. Open a command terminal if you haven't done so.
2. [Linux and Mac Users Only] The conda package manager only works with bash,
so if bash is not your current shell, type::
bash
3. Add the conda-forge channel to your conda package manager.
Type or copy the command below in to your command terminal. You should
run this command even if you have already done it in the past.
This will ensure that conda-forge is set as the highest priority channel.
::
conda config --add channels conda-forge
.. note::
Conda-forge is a community driven collection of packages that are
continually tested to ensure compatibility. We highly recommend using
conda-forge when working with conda. See https://conda-forge.github.io/
for more details on this excellent project.
4. Create the conda environment for the tutorial.
Type or copy this command in to your command terminal::
conda create -n tutorial_2018 python=2.7 matplotlib cartopy netcdf4 jupyter git ffmpeg wrf-python
Type "y" when prompted. It will take several minutes to install everything.
This command creates an isolated Python environment named *tutorial_2018*, and installs
the python interpreter, matplotlib, cartopy, netcdf4, jupyter, git, ffmpeg, and wrf-python
packages.
.. note::
When the installation completes, your command terminal might post a message similar to:
.. code-block:: none
If this is your first install of dbus, automatically load on login with:
mkdir -p ~/Library/LaunchAgents
cp /path/to/miniconda2/envs/tutorial_test/org.freedesktop.dbus-session.plist ~/Library/LaunchAgents/
launchctl load -w ~/Library/LaunchAgents/org.freedesktop.dbus-session.plist
This is indicating that the dbus package can be set up to automatically load on login. You
can either ignore this message or type in the commands as indicated on your command terminal.
The tutorial should work fine in either case.
5. Activate the conda environment.
To activate the tutorial_2018 Python environment, type the following
in to the command terminal:
For Linux and Mac (using bash)::
source activate tutorial_2018
For Windows::
activate tutorial_2018
You should see (tutorial_2018) on your command prompt.
To deactivate your conda environment, type the following in to the
command terminal:
For Linux and Mac::
source deactivate
For Windows::
deactivate tutorial_2018
Step 5: Download the Student Workbook
---------------------------------------
The student workbook for the tutorial is available on GitHub. The tutorial_2018
conda environment includes the git application needed to download the repository.
These instructions download the tutorial in to your home directory. If you want
to place the tutorial in to another directory, we're going to assume you know
how to do this yourself.
To download the student workbook, follow these instructions:
1. Activate the tutorial_2018 conda environment following the instructions
in the previous step (*source activate tutorial_2018* or
*activate tutorial_2018*).
2. Change your working directory to the home directory by typing the
following command in to the command terminal:
For Linux and Mac::
cd ~
For Windows::
cd %HOMEPATH%
3. Download the git repository for the tutorial by typing the following
in to the command terminal::
git clone https://github.com/NCAR/wrf_python_tutorial.git
4. There may be additional changes to the tutorial after you have downloaded
it. To pull down the latest changes, type the following in to the
command terminal:
For Linux and Mac::
source activate tutorial_2018
cd ~/wrf_python_tutorial/wrf_workshop_2018
git pull
For Windows::
activate tutorial_2018
cd %HOMEPATH%\wrf_python_tutorial\wrf_workshop_2018
git pull
.. note::
If you try the "git pull" command and it returns an error indicating
that you have made changes to the workbook, this is probably because
you ran the workbook and it contains the cell output. To fix this,
first do a checkout of the workbook, then do the pull.
.. code-block:: none
git checkout -- .
git pull
Step 6: Verify Your Environment
----------------------------------
Verifying that your environment is correct involves importing a few
packages and checking for errors (you may see some warnings for matplotlib
or xarray, but you can safely ignore these).
1. Activate the tutorial_2018 conda environment if it isn't already active
(see instructions above).
2. Open a python terminal by typing the following in to the command
terminal::
python
3. Now type the following in to the Python interpreter::
>>> import netCDF4
>>> import matplotlib
>>> import xarray
>>> import wrf
4. You can exit the Python interpreter using **CTRL + D**
Step 7: Obtain WRF Output Files
----------------------------------
For this tutorial, we strongly recommend that you use your own WRF output files.
The tutorial includes an easy way to point to your own data files. The WRF
output files should all be from the same WRF run and use the same domain.
If your files are located on another system (e.g. yellowstone), then copy 2 or
3 of these files to your local computer prior to the tutorial.
If you do not have any of your own WRF output files, then you can download the
instructor data files from a link that should have been provided to you in an
email prior to the tutorial.
If you are using the link provided in the email for your data, you can follow
the instructions below to place your data in the default location for your
workbook.
1. The link in the email should take you to a location on an Amazon cloud
drive.
2. If you hover your mouse over the wrf_tutorial_data.zip file, you'll see
an empty check box appear next to the file name. Click this check
box.
3. At the bottom of the screen, you'll see a Download button next to a
cloud icon. Click this button to start the download.
4. The download was most likely placed in to your ~/Downloads folder
[%HOMEPATH%\\Downloads for Windows]. Using your preferred method of choice
for unzipping files, unzip this file in to your home directory. Your data
should now be in ~/wrf_tutorial_data
[%HOMEPATH%\\wrf_tutorial_data for Windows].
5. Verify that you have three WRF output files in that directory.
Getting Help
----------------
If you experience problems during this installation, please send a question
to the :ref:`google-group` support mailing list.
We look forward to seeing you at the tutorial!

28
fortran/rip_cape.f90

@ -272,6 +272,7 @@ END SUBROUTINE DPFCALC @@ -272,6 +272,7 @@ END SUBROUTINE DPFCALC
! NCLFORTSTART
SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
prsf, prs_new, tmk_new, qvp_new, ght_new,&
cmsg,mix,mjy,mkzh,ter_follow,&
psafile, errstat, errmsg)
USE wrf_constants, ONLY : CELKEL, G, EZERO, ESLCON1, ESLCON2, &
@ -293,6 +294,14 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -293,6 +294,14 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) ::sfp
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cape
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cin
REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prsf
REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prs_new
REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: tmk_new
REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: qvp_new
REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: ght_new
REAL(KIND=8), INTENT(IN) :: cmsg
CHARACTER(LEN=*), INTENT(IN) :: psafile
INTEGER, INTENT(INOUT) :: errstat
@ -308,15 +317,10 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -308,15 +317,10 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat
REAL(KIND=8) :: benamin, dz
REAL(KIND=8), DIMENSION(150) :: buoy, zrel, benaccum
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prsf
REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs
REAL(KIND=8), DIMENSION(150,150) :: psaditmk
LOGICAL :: elfound
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prs_new
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: tmk_new
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: qvp_new
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: ght_new
! To remove compiler warnings
tmkpari = 0
@ -597,6 +601,7 @@ END SUBROUTINE DCAPECALC3D @@ -597,6 +601,7 @@ END SUBROUTINE DCAPECALC3D
! NCLFORTSTART
SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
prsf, prs_new, tmk_new, qvp_new, ght_new,&
cmsg,mix,mjy,mkzh,ter_follow,&
psafile, errstat, errmsg)
USE wrf_constants, ONLY : CELKEL, G, EZERO, ESLCON1, ESLCON2, &
@ -618,6 +623,13 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -618,6 +623,13 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) ::sfp
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cape
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cin
REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prsf
REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prs_new
REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: tmk_new
REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: qvp_new
REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: ght_new
REAL(KIND=8), INTENT(IN) :: cmsg
CHARACTER(LEN=*), INTENT(IN) :: psafile
INTEGER, INTENT(INOUT) :: errstat
@ -635,15 +647,11 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& @@ -635,15 +647,11 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat
REAL(KIND=8) :: benamin, dz, pup, pdn
REAL(KIND=8), DIMENSION(150) :: buoy, zrel, benaccum
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prsf
REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs
REAL(KIND=8), DIMENSION(150,150) :: psaditmk
LOGICAL :: elfound
REAL(KIND=8), DIMENSION(mkzh) :: eth_temp
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prs_new
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: tmk_new
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: qvp_new
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: ght_new
! To remove compiler warnings
errstat = 0

56
fortran/wrf_fctt.f90

@ -1,5 +1,6 @@ @@ -1,5 +1,6 @@
!NCLFORTSTART
SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew)
SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci,&
fill_nocloud, missing, opt_thresh, nz, ns, ew)
USE wrf_constants, ONLY : EPS, USSALR, RD, G, ABSCOEFI, ABSCOEF, CELKEL
IMPLICIT NONE
@ -7,10 +8,14 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew @@ -7,10 +8,14 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew
!f2py threadsafe
!f2py intent(in,out) :: ctt
INTEGER, INTENT(IN) :: nz, ns, ew, haveqci
INTEGER, INTENT(IN) :: nz, ns, ew, haveqci, fill_nocloud
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: ght, prs, tk, qci, qcw, qvp
REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: ter
REAL(KIND=8), DIMENSION(ew,ns), INTENT(OUT) :: ctt
REAL(KIND=8), INTENT(IN) :: missing
REAL(KIND=8), INTENT(IN) :: opt_thresh
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(INOUT) :: pf
!NCLEND
@ -20,7 +25,7 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew @@ -20,7 +25,7 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew
INTEGER i,j,k,ripk
REAL(KIND=8) :: opdepthu, opdepthd, dp, arg1, fac, prsctt, ratmix
REAL(KIND=8) :: arg2, agl_hgt, vt
REAL(KIND=8), DIMENSION(ew,ns,nz) :: pf
REAL(KIND=8) :: p1, p2
!$OMP PARALLEL
@ -58,12 +63,12 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew @@ -58,12 +63,12 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew
DO i=1,ew
opdepthd = 0.D0
k = 0
prsctt = 0
prsctt = -1
! Integrate downward from model top, calculating path at full
! model vertical levels.
DO k=1, nz
DO k=2,nz
opdepthu = opdepthd
ripk = nz - k + 1
@ -74,21 +79,23 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew @@ -74,21 +79,23 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew
END IF
IF (haveqci .EQ. 0) then
IF (tk(i,j,k) .LT. CELKEL) then
IF (tk(i,j,ripk) .LT. CELKEL) then
! Note: abscoefi is m**2/g, qcw is g/kg, so no convrsion needed
opdepthd = opdepthu + ABSCOEFI*qcw(i,j,k) * dp/G
opdepthd = opdepthu + ABSCOEFI*qcw(i,j,ripk) * dp/G
ELSE
opdepthd = opdepthu + ABSCOEF*qcw(i,j,k) * dp/G
opdepthd = opdepthu + ABSCOEF*qcw(i,j,ripk) * dp/G
END IF
ELSE
opdepthd = opdepthd + (ABSCOEF*qcw(i,j,ripk) + ABSCOEFI*qci(i,j,ripk))*dp/G
END IF
IF (opdepthd .LT. 1. .AND. k .LT. nz) THEN
IF (opdepthd .LT. opt_thresh .AND. k .LT. nz) THEN
CYCLE
ELSE IF (opdepthd .LT. 1. .AND. k .EQ. nz) THEN
prsctt = prs(i,j,1)
ELSE IF (opdepthd .LT. opt_thresh .AND. k .EQ. nz) THEN
IF (fill_nocloud .EQ. 0) THEN
prsctt = prs(i,j,1)
ENDIF
EXIT
ELSE
fac = (1. - opdepthu)/(opdepthd - opdepthu)
@ -98,17 +105,22 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew @@ -98,17 +105,22 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew
END IF
END DO
DO k=2,nz
ripk = nz - k + 1
p1 = prs(i,j,ripk+1)
p2 = prs(i,j,ripk)
IF (prsctt .GE. p1 .AND. prsctt .LE. p2) THEN
fac = (prsctt - p1)/(p2 - p1)
arg1 = fac*(tk(i,j,ripk) - tk(i,j,ripk+1)) - CELKEL
ctt(i,j) = tk(i,j,ripk+1) + arg1
EXIT
END IF
END DO
! prsctt should only be 0 if fill values are used
IF (prsctt .GT. -1) THEN
DO k=2,nz
ripk = nz - k + 1
p1 = prs(i,j,ripk+1)
p2 = prs(i,j,ripk)
IF (prsctt .GE. p1 .AND. prsctt .LE. p2) THEN
fac = (prsctt - p1)/(p2 - p1)
arg1 = fac*(tk(i,j,ripk) - tk(i,j,ripk+1)) - CELKEL
ctt(i,j) = tk(i,j,ripk+1) + arg1
EXIT
END IF
END DO
ELSE
ctt(i,j) = missing
END IF
END DO
END DO
!$OMP END DO

6
fortran/wrf_vinterp.f90

@ -127,7 +127,7 @@ END FUNCTION wrf_intrp_value @@ -127,7 +127,7 @@ END FUNCTION wrf_intrp_value
!NCLFORTSTART
SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
sfp, smsfp, vcarray, interp_levels, numlevels,&
icase, ew, ns, nz, extrap, vcor, logp, rmsg,&
icase, ew, ns, nz, extrap, vcor, logp, tempout, rmsg,&
errstat, errmsg)
USE wrf_constants, ONLY : ALGERR, SCLHT, EXPON, EXPONI, GAMMA, GAMMAMD, TLCLC1, &
TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3, &
@ -146,6 +146,9 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& @@ -146,6 +146,9 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
REAL(KIND=8), DIMENSION(ew,ns,numlevels), INTENT(OUT) :: dataout
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: vcarray
REAL(KIND=8), DIMENSION(numlevels), INTENT(IN) :: interp_levels
REAL(KIND=8), DIMENSION(ew,ns), INTENT(INOUT) :: tempout
REAL(KIND=8), INTENT(IN) :: rmsg
INTEGER, INTENT(INOUT) :: errstat
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg
@ -156,7 +159,6 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& @@ -156,7 +159,6 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
INTEGER :: i, j, k, kupper !itriv
INTEGER :: ifound, isign !miy,mjx
INTEGER :: log_errcnt, interp_errcnt, interp_errstat
REAL(KIND=8), DIMENSION(ew,ns) :: tempout
REAL(KIND=8) :: rlevel, vlev, diff
REAL(KIND=8) :: tmpvlev
REAL(KIND=8) :: vcp1, vcp0, valp0, valp1

23
fortran/wrffortran.pyf

@ -231,7 +231,7 @@ python module _wrffortran ! in @@ -231,7 +231,7 @@ python module _wrffortran ! in
integer, optional,intent(in),check(shape(prs,0)==mkzh),depend(prs) :: mkzh=shape(prs,0)
integer intent(in) :: ter_follow
end subroutine dpfcalc
subroutine dcapecalc3d(prs,tmk,qvp,ght,ter,sfp,cape,cin,cmsg,mix,mjy,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90
subroutine dcapecalc3d(prs,tmk,qvp,ght,ter,sfp,cape,cin,prsf,prs_new,tmk_new,qvp_new,ght_new,cmsg,mix,mjy,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90
threadsafe
use wrf_constants, only: tlclc2,gamma,tlclc1,rgasmd,tlclc4,g,tlclc3,thtecon3,eps,rd,cpmd,celkel,gammamd,eslcon2,eslcon1,cp,thtecon1,ezero,thtecon2
real(kind=8) dimension(mix,mjy,mkzh),intent(in) :: prs
@ -242,6 +242,11 @@ python module _wrffortran ! in @@ -242,6 +242,11 @@ python module _wrffortran ! in
real(kind=8) dimension(mix,mjy),intent(in),depend(mix,mjy) :: sfp
real(kind=8) dimension(mix,mjy,mkzh),intent(out,in),depend(mix,mjy,mkzh) :: cape
real(kind=8) dimension(mix,mjy,mkzh),intent(out,in),depend(mix,mjy,mkzh) :: cin
real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: prsf
real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: prs_new
real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: tmk_new
real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: qvp_new
real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: ght_new
real(kind=8) intent(in) :: cmsg
integer, optional,intent(in),check(shape(prs,0)==mix),depend(prs) :: mix=shape(prs,0)
integer, optional,intent(in),check(shape(prs,1)==mjy),depend(prs) :: mjy=shape(prs,1)
@ -251,7 +256,7 @@ python module _wrffortran ! in @@ -251,7 +256,7 @@ python module _wrffortran ! in
integer intent(inout) :: errstat
character*(*) intent(inout) :: errmsg
end subroutine dcapecalc3d
subroutine dcapecalc2d(prs,tmk,qvp,ght,ter,sfp,cape,cin,cmsg,mix,mjy,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90
subroutine dcapecalc2d(prs,tmk,qvp,ght,ter,sfp,cape,cin,prsf,prs_new,tmk_new,qvp_new,ght_new,cmsg,mix,mjy,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90
threadsafe
use wrf_constants, only: tlclc2,gamma,tlclc1,rgasmd,tlclc4,g,tlclc3,thtecon3,eps,rd,cpmd,celkel,gammamd,eslcon2,eslcon1,cp,thtecon1,ezero,thtecon2
real(kind=8) dimension(mix,mjy,mkzh),intent(in) :: prs
@ -262,6 +267,11 @@ python module _wrffortran ! in @@ -262,6 +267,11 @@ python module _wrffortran ! in
real(kind=8) dimension(mix,mjy),intent(in),depend(mix,mjy) :: sfp
real(kind=8) dimension(mix,mjy,mkzh),intent(out,in),depend(mix,mjy,mkzh) :: cape
real(kind=8) dimension(mix,mjy,mkzh),intent(out,in),depend(mix,mjy,mkzh) :: cin
real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: prsf
real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: prs_new
real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: tmk_new
real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: qvp_new
real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: ght_new
real(kind=8) intent(in) :: cmsg
integer, optional,intent(in),check(shape(prs,0)==mix),depend(prs) :: mix=shape(prs,0)
integer, optional,intent(in),check(shape(prs,1)==mjy),depend(prs) :: mjy=shape(prs,1)
@ -349,7 +359,7 @@ python module _wrffortran ! in @@ -349,7 +359,7 @@ python module _wrffortran ! in
real(kind=8), parameter,optional,depend(cp,rd) :: gamma=0.285714285714
real(kind=8), parameter,optional,depend(expon,rd,ussalr,g) :: exponi=5.25864379523
end module wrf_constants
subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt,haveqci,nz,ns,ew) ! in :_wrffortran:wrf_fctt.f90
subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt,pf,haveqci,fill_nocloud,missing,opt_thresh,nz,ns,ew) ! in :_wrffortran:wrf_fctt.f90
threadsafe
use wrf_constants, only: g,eps,rd,ussalr,abscoefi,abscoef,celkel
real(kind=8) dimension(ew,ns,nz),intent(in) :: prs
@ -360,7 +370,11 @@ python module _wrffortran ! in @@ -360,7 +370,11 @@ python module _wrffortran ! in
real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: ght
real(kind=8) dimension(ew,ns),intent(in),depend(ew,ns) :: ter
real(kind=8) dimension(ew,ns),intent(out,in),depend(ew,ns) :: ctt
real(kind=8) dimension(ew,ns,nz),intent(inout),depend(ew,ns,nz) :: pf
integer intent(in) :: haveqci
integer intent(in) :: fill_nocloud
real(kind=8) intent(in) :: missing
real(kind=8) intent(in) :: opt_thresh
integer, optional,intent(in),check(shape(prs,2)==nz),depend(prs) :: nz=shape(prs,2)
integer, optional,intent(in),check(shape(prs,1)==ns),depend(prs) :: ns=shape(prs,1)
integer, optional,intent(in),check(shape(prs,0)==ew),depend(prs) :: ew=shape(prs,0)
@ -730,7 +744,7 @@ python module _wrffortran ! in @@ -730,7 +744,7 @@ python module _wrffortran ! in
integer intent(inout) :: errstat
real(kind=8) :: wrf_intrp_value
end function wrf_intrp_value
subroutine wrf_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,sfp,smsfp,vcarray,interp_levels,numlevels,icase,ew,ns,nz,extrap,vcor,logp,rmsg,errstat,errmsg) ! in :_wrffortran:wrf_vinterp.f90
subroutine wrf_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,sfp,smsfp,vcarray,interp_levels,numlevels,icase,ew,ns,nz,extrap,vcor,logp,tempout,rmsg,errstat,errmsg) ! in :_wrffortran:wrf_vinterp.f90
threadsafe
use wrf_constants, only: tlclc2,tlclc3,sclht,tlclc4,thtecon3,eps,tlclc1,celkel,exponi,gammamd,ussalr,expon,thtecon1,algerr,gamma,thtecon2
real(kind=8) dimension(ew,ns,nz),intent(in) :: datain
@ -752,6 +766,7 @@ python module _wrffortran ! in @@ -752,6 +766,7 @@ python module _wrffortran ! in
integer intent(in) :: extrap
integer intent(in) :: vcor
integer intent(in) :: logp
real(kind=8) dimension(ew,ns),intent(inout),depend(ew,ns) :: tempout
real(kind=8) intent(in) :: rmsg
integer intent(inout) :: errstat
character*(*) intent(inout) :: errmsg

73
src/wrf/computation.py

@ -714,10 +714,10 @@ def smooth2d(field, passes, meta=True): @@ -714,10 +714,10 @@ def smooth2d(field, passes, meta=True):
@set_cape_alg_metadata(is2d=True, copyarg="pres_hpa")
def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow,
missing=default_fill(np.float64), meta=True):
"""Return the two-dimensional CAPE, CIN, LCL, and LFC.
"""Return the two-dimensional MCAPE, MCIN, LCL, and LFC.
This function calculates the maximum convective available potential
energy (CAPE), maximum convective inhibition (CIN),
energy (MCAPE), maximum convective inhibition (MCIN),
lifted condensation level (LCL), and level of free convection (LFC). This
function uses the RIP [Read/Interpolate/plot] code to calculate
potential energy (CAPE) and convective inhibition
@ -725,18 +725,28 @@ def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow, @@ -725,18 +725,28 @@ def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow,
in the column (i.e. something akin to Colman's MCAPE). CAPE is defined as
the accumulated buoyant energy from the level of free convection (LFC) to
the equilibrium level (EL). CIN is defined as the accumulated negative
buoyant energy from the parcel starting point to the LFC. The word 'parcel'
here refers to a 500 meter deep parcel, with actual temperature and
moisture averaged over that depth.
buoyant energy from the parcel starting point to the LFC.
The cape_2d algorithm works by first finding the maximum theta-e height
level in the lowest 3000 m. A parcel with a depth of 500 m is then
calculated and centered over this maximum theta-e height level. The
parcel's moisture and temperature characteristics are calculated by
averaging over the depth of this 500 m parcel. This 'maximum' parcel
is then used to compute MCAPE, MCIN, LCL and LFC.
The leftmost dimension of the returned array represents four different
quantities:
- return_val[0,...] will contain CAPE [J kg-1]
- return_val[1,...] will contain CIN [J kg-1]
- return_val[0,...] will contain MCAPE [J kg-1]
- return_val[1,...] will contain MCIN [J kg-1]
- return_val[2,...] will contain LCL [m]
- return_val[3,...] will contain LFC [m]
This function also supports computing MCAPE along a single vertical
column. In this mode, the *pres_hpa*, *tkel*, *qv* and *height* arguments
must be one-dimensional vertical columns, and the *terrain* and
*psfc_hpa* arguments must be scalar values
(:obj:`float`, :class:`numpy.float32` or :class:`numpy.float64`).
This is the raw computational algorithm and does not extract any variables
from WRF output files. Use :meth:`wrf.getvar` to both extract and compute
@ -749,6 +759,9 @@ def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow, @@ -749,6 +759,9 @@ def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow,
least three dimensions. The rightmost dimensions can be
top_bottom x south_north x west_east or bottom_top x south_north x
west_east.
When operating on only a single column of values, the vertical
column can be bottom_top or top_bottom. In this case, *terrain*
and *psfc_hpa* must be scalars.
Note:
@ -774,12 +787,16 @@ def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow, @@ -774,12 +787,16 @@ def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow,
terrain (:class:`xarray.DataArray` or :class:`numpy.ndarray`):
Terrain height in [m]. This is at least a two-dimensional array
with the same dimensionality as *pres_hpa*, excluding the vertical
(bottom_top/top_bottom) dimension.
(bottom_top/top_bottom) dimension. When operating on a single
vertical column, this argument must be a scalar (:obj:`float`,
:class:`numpy.float32`, or :class:`numpy.float64`).
psfc_hpa (:class:`xarray.DataArray` or :class:`numpy.ndarray`):
The surface pressure in [hPa]. This is at least a two-dimensional
array with the same dimensionality as *pres_hpa*, excluding the
vertical (bottom_top/top_bottom) dimension.
vertical (bottom_top/top_bottom) dimension. When operating on a
singlevertical column, this argument must be a scalar
(:obj:`float`, :class:`numpy.float32`, or :class:`numpy.float64`).
Note:
@ -1014,6 +1031,11 @@ def cloudfrac(vert, relh, vert_inc_w_height, low_thresh, mid_thresh, @@ -1014,6 +1031,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 +1069,9 @@ def cloudfrac(vert, relh, vert_inc_w_height, low_thresh, mid_thresh, @@ -1047,8 +1069,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 +1117,27 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True, @@ -1094,6 +1117,27 @@ 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 amount of optical
depth (integrated from top down) required to trigger a cloud top
temperature calculation. The cloud top temperature is calculated at
the vertical level where this threshold is met. Vertical columns
with less than this threshold will be treated as cloud free areas.
In general, the larger the value is for this
threshold, the lower the altitude will be for the cloud top
temperature calculation, and therefore higher cloud top
temperature values. Default is 1.0, which should be sufficient for
most users.
meta (:obj:`bool`): Set to False to disable metadata and return
:class:`numpy.ndarray` instead of
@ -1128,8 +1172,13 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True, @@ -1128,8 +1172,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)

55
src/wrf/extension.py

@ -618,20 +618,33 @@ def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow, @@ -618,20 +618,33 @@ def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow,
else:
cape_routine = dcapecalc2d
# Work arrays
k_left_shape = (p_hpa.shape[2], p_hpa.shape[0], p_hpa.shape[1])
prsf = np.empty(k_left_shape, np.float64, order="F")
prs_new = np.empty(k_left_shape, np.float64, order="F")
tmk_new = np.empty(k_left_shape, np.float64, order="F")
qvp_new = np.empty(k_left_shape, np.float64, order="F")
ght_new = np.empty(k_left_shape, np.float64, order="F")
# note that p_hpa, tk, qv, and ht have the vertical flipped
result = cape_routine(p_hpa,
tk,
qv,
ht,
ter,
sfp,
capeview,
cinview,
missing,
ter_follow,
psafile,
errstat,
errmsg)
tk,
qv,
ht,
ter,
sfp,
capeview,
cinview,
prsf,
prs_new,
tmk_new,
qvp_new,
ght_new,
missing,
ter_follow,
psafile,
errstat,
errmsg)
if int(errstat) != 0:
raise DiagnosticError("".join(npbytes_to_str(errmsg)).strip())
@ -754,10 +767,11 @@ def _xytoll(map_proj, truelat1, truelat2, stdlon, lat1, lon1, @@ -754,10 +767,11 @@ def _xytoll(map_proj, truelat1, truelat2, stdlon, lat1, lon1,
@check_args(0, 3, (3,3,3,3,3,3,2))
@left_iteration(3, 2, ref_var_idx=0, ignore_args=(7,))
@left_iteration(3, 2, ref_var_idx=0, ignore_args=(7,8,9,10))
@cast_type(arg_idxs=(0,1,2,3,4,5,6))
@extract_and_transpose()
def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None):
def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, fill_nocloud,
missing, opt_thresh, outview=None):
"""Wrapper for wrfcttcalc.
Located in wrf_fctt.f90.
@ -766,6 +780,8 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): @@ -766,6 +780,8 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None):
if outview is None:
outview = np.empty_like(ter)
pf = np.empty(p_hpa.shape[0:3], np.float64, order="F")
result = wrfcttcalc(p_hpa,
tk,
qice,
@ -774,7 +790,11 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): @@ -774,7 +790,11 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None):
ght,
ter,
outview,
haveqci)
pf,
haveqci,
fill_nocloud,
missing,
opt_thresh)
return result
@ -858,7 +878,9 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, @@ -858,7 +878,9 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp,
if outview is None:
outdims = field.shape[0:2] + interp_levels.shape
outview = np.empty(outdims, field.dtype, order="F")
tempout = np.zeros(field.shape[0:2], np.float64, order="F")
errstat = np.array(0)
errmsg = np.zeros(Constants.ERRLEN, "c")
@ -877,6 +899,7 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, @@ -877,6 +899,7 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp,
extrap,
vcor,
logp,
tempout,
missing,
errstat,
errmsg)

6
src/wrf/g_cape.py

@ -13,13 +13,13 @@ from .metadecorators import set_cape_metadata @@ -13,13 +13,13 @@ from .metadecorators import set_cape_metadata
@set_cape_metadata(is2d=True)
def get_2dcape(wrfin, timeidx=0, method="cat", squeeze=True, cache=None,
meta=True, _key=None, missing=default_fill(np.float64)):
"""Return the 2d fields of CAPE, CIN, LCL, and LFC.
"""Return the 2d fields of MCAPE, MCIN, LCL, and LFC.
The leftmost dimension of the returned array represents four different
quantities:
- return_val[0,...] will contain CAPE [J kg-1]
- return_val[1,...] will contain CIN [J kg-1]
- return_val[0,...] will contain MCAPE [J kg-1]
- return_val[1,...] will contain MCIN [J kg-1]
- return_val[2,...] will contain LCL [m]
- return_val[3,...] will contain LFC [m]

38
src/wrf/g_ctt.py

@ -1,11 +1,12 @@ @@ -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 @@ -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", @@ -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", @@ -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", @@ -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)

8
src/wrf/g_dewpoint.py

@ -76,7 +76,9 @@ def get_dp(wrfin, timeidx=0, method="cat", squeeze=True, @@ -76,7 +76,9 @@ def get_dp(wrfin, timeidx=0, method="cat", squeeze=True,
p = ncvars["P"]
pb = ncvars["PB"]
qvapor = ncvars["QVAPOR"]
# Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to
# break with every release
qvapor = ncvars["QVAPOR"].copy()
# Algorithm requires hPa
full_p = .01*(p + pb)
@ -152,7 +154,9 @@ def get_dp_2m(wrfin, timeidx=0, method="cat", squeeze=True, @@ -152,7 +154,9 @@ def get_dp_2m(wrfin, timeidx=0, method="cat", squeeze=True,
# Algorithm requires hPa
psfc = .01*(ncvars["PSFC"])
q2 = ncvars["Q2"]
# Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to
# break with every release
q2 = ncvars["Q2"].copy()
q2[q2 < 0] = 0
td = _td(psfc, q2)

2
src/wrf/g_helicity.py

@ -176,7 +176,7 @@ def get_uh(wrfin, timeidx=0, method="cat", squeeze=True, @@ -176,7 +176,7 @@ def get_uh(wrfin, timeidx=0, method="cat", squeeze=True,
"""
ncvars = extract_vars(wrfin, timeidx, ("W", "PH", "PHB", "MAPFAC_M"),
method, squeeze, cache, meta=False, _key=None)
method, squeeze, cache, meta=False, _key=_key)
wstag = ncvars["W"]
ph = ncvars["PH"]

131
src/wrf/g_latlon.py

@ -1,9 +1,16 @@ @@ -1,9 +1,16 @@
from __future__ import (absolute_import, division, print_function,
unicode_literals)
from .util import extract_vars, get_id
from .latlonutils import (_lat_varname, _lon_varname, _ll_to_xy, _xy_to_ll)
from collections import OrderedDict
import numpy as np
from xarray import DataArray
from .util import extract_vars, get_id, get_iterable, is_mapping, to_np
from .py3compat import viewkeys
from .latlonutils import _lat_varname, _lon_varname, _ll_to_xy, _xy_to_ll
from .metadecorators import set_latlon_metadata
from .config import xarray_enabled
def get_lat(wrfin, timeidx=0, method="cat", squeeze=True,
@ -151,7 +158,101 @@ def get_lon(wrfin, timeidx=0, method="cat", squeeze=True, @@ -151,7 +158,101 @@ def get_lon(wrfin, timeidx=0, method="cat", squeeze=True,
return lon_var[varname]
# TODO: Do we need the user to know about method, squeeze, cache for this?
def _llxy_mapping(wrfin, x_or_lat, y_or_lon, func, timeidx, stagger,
squeeze, meta, as_int=None):
keynames = []
# This might not work once mapping iterators are implemented
numkeys = len(wrfin)
key_iter = iter(viewkeys(wrfin))
first_key = next(key_iter)
keynames.append(first_key)
first_args = [wrfin[first_key], x_or_lat, y_or_lon, timeidx, squeeze,
meta, stagger]
if as_int is not None:
first_args.append(as_int)
first_array = func(*first_args)
# Create the output data numpy array based on the first array
outdims = [numkeys]
outdims += first_array.shape
outdata = np.empty(outdims, first_array.dtype)
outdata[0,:] = first_array[:]
idx = 1
while True:
try:
key = next(key_iter)
except StopIteration:
break
else:
keynames.append(key)
args = [wrfin[first_key], x_or_lat, y_or_lon, timeidx, squeeze,
meta, stagger]
if as_int is not None:
args.append(as_int)
result_array = func(*args)
if outdata.shape[1:] != result_array.shape:
raise ValueError("data sequences must have the "
"same size for all dictionary keys")
outdata[idx,:] = to_np(result_array)[:]
idx += 1
if xarray_enabled() and meta:
outname = str(first_array.name)
# Note: assumes that all entries in dict have same coords
outcoords = OrderedDict(first_array.coords)
# First find and store all the existing key coord names/values
# This is applicable only if there are nested dictionaries.
key_coordnames = []
coord_vals = []
existing_cnt = 0
while True:
key_coord_name = "key_{}".format(existing_cnt)
if key_coord_name not in first_array.dims:
break
key_coordnames.append(key_coord_name)
coord_vals.append(to_np(first_array.coords[key_coord_name]))
existing_cnt += 1
# Now add the key coord name and values for THIS dictionary.
# Put the new key_n name at the bottom, but the new values will
# be at the top to be associated with key_0 (left most). This
# effectively shifts the existing 'key_n' coordinate values to the
# right one dimension so *this* dicionary's key coordinate values
# are at 'key_0'.
key_coordnames.append(key_coord_name)
coord_vals.insert(0, keynames)
# make it so that key_0 is leftmost
outdims = key_coordnames + list(first_array.dims[existing_cnt:])
# Create the new 'key_n', value pairs
for coordname, coordval in zip(key_coordnames, coord_vals):
outcoords[coordname] = coordval
outattrs = OrderedDict(first_array.attrs)
outarr = DataArray(outdata, name=outname, coords=outcoords,
dims=outdims, attrs=outattrs)
else:
outarr = outdata
return outarr
@set_latlon_metadata(xy=True)
def ll_to_xy(wrfin, latitude, longitude, timeidx=0,
@ -215,9 +316,15 @@ def ll_to_xy(wrfin, latitude, longitude, timeidx=0, @@ -215,9 +316,15 @@ def ll_to_xy(wrfin, latitude, longitude, timeidx=0,
be a :class:`numpy.ndarray` object with no metadata.
"""
if is_mapping(wrfin):
return _llxy_mapping(wrfin, latitude, longitude, ll_to_xy,
timeidx, stagger, squeeze, meta, as_int)
_key = get_id(wrfin)
return _ll_to_xy(latitude, longitude, wrfin, timeidx, stagger, "cat",
_wrfin = get_iterable(wrfin)
return _ll_to_xy(latitude, longitude, _wrfin, timeidx, stagger, "cat",
squeeze, None, _key, as_int, **{})
@set_latlon_metadata(xy=True)
@ -319,10 +426,10 @@ def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True, @@ -319,10 +426,10 @@ def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, as_int=True,
return _ll_to_xy(latitude, longitude, None, 0, True, "cat", squeeze, None,
None, as_int, **projparams)
@set_latlon_metadata(xy=False)
def xy_to_ll(wrfin, x, y, timeidx=0, stagger=None, squeeze=True, meta=True):
def xy_to_ll(wrfin, x, y, timeidx=0, squeeze=True, meta=True, stagger=None):
"""Return the latitude and longitude for specified x,y coordinates.
The *x* and *y* arguments can be a single value or a sequence of values.
@ -370,9 +477,6 @@ def xy_to_ll(wrfin, x, y, timeidx=0, stagger=None, squeeze=True, meta=True): @@ -370,9 +477,6 @@ def xy_to_ll(wrfin, x, y, timeidx=0, stagger=None, squeeze=True, meta=True):
- 'v': Use the same staggered grid as the v wind component,
which has a staggered south_north (y) dimension.
as_int (:obj:`bool`): Set to True to return the x,y values as
:obj:`int`, otherwise they will be returned as :obj:`float`.
Returns:
:class:`xarray.DataArray` or :class:`numpy.ndarray`: The
latitude and longitude values whose leftmost dimension is 2
@ -382,8 +486,13 @@ def xy_to_ll(wrfin, x, y, timeidx=0, stagger=None, squeeze=True, meta=True): @@ -382,8 +486,13 @@ def xy_to_ll(wrfin, x, y, timeidx=0, stagger=None, squeeze=True, meta=True):
be a :class:`numpy.ndarray` object with no metadata.
"""
if is_mapping(wrfin):
return _llxy_mapping(wrfin, x, y, xy_to_ll,
timeidx, stagger, squeeze, meta)
_key = get_id(wrfin)
return _xy_to_ll(x, y, wrfin, timeidx, stagger, "cat", True, None,
_wrfin = get_iterable(wrfin)
return _xy_to_ll(x, y, _wrfin, timeidx, stagger, "cat", True, None,
_key, **{})

8
src/wrf/g_rh.py

@ -71,7 +71,9 @@ def get_rh(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, @@ -71,7 +71,9 @@ def get_rh(wrfin, timeidx=0, method="cat", squeeze=True, cache=None,
t = ncvars["T"]
p = ncvars["P"]
pb = ncvars["PB"]
qvapor = ncvars["QVAPOR"]
# Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to
# break with every release
qvapor = ncvars["QVAPOR"].copy()
full_t = t + Constants.T_BASE
full_p = p + pb
@ -144,7 +146,9 @@ def get_rh_2m(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, @@ -144,7 +146,9 @@ def get_rh_2m(wrfin, timeidx=0, method="cat", squeeze=True, cache=None,
meta=False, _key=_key)
t2 = ncvars["T2"]
psfc = ncvars["PSFC"]
q2 = ncvars["Q2"]
# Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to
# break with every release
q2 = ncvars["Q2"].copy()
q2[q2 < 0] = 0
rh = _rh(q2, psfc, t2)

4
src/wrf/g_slp.py

@ -83,7 +83,9 @@ def get_slp(wrfin, timeidx=0, method="cat", squeeze=True, @@ -83,7 +83,9 @@ def get_slp(wrfin, timeidx=0, method="cat", squeeze=True,
t = ncvars["T"]
p = ncvars["P"]
pb = ncvars["PB"]
qvapor = ncvars["QVAPOR"]
# Copy needed for the mmap nonsense of scipy.io.netcdf, which seems to
# break with every release
qvapor = ncvars["QVAPOR"].copy()
ph = ncvars["PH"]
phb = ncvars["PHB"]

22
src/wrf/interp.py

@ -8,7 +8,7 @@ from .extension import (_interpz3d, _vertcross, _interpline, _smooth2d, @@ -8,7 +8,7 @@ from .extension import (_interpz3d, _vertcross, _interpline, _smooth2d,
_monotonic, _vintrp)
from .metadecorators import set_interp_metadata
from .util import extract_vars, is_staggered, get_id, to_np
from .util import extract_vars, is_staggered, get_id, to_np, get_iterable
from .py3compat import py3range
from .interputils import get_xy, get_xy_z_params, to_xy_coords
from .constants import Constants, default_fill, ConversionFactors
@ -548,6 +548,8 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, @@ -548,6 +548,8 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
"""
_key = get_id(wrfin)
_wrfin = get_iterable(wrfin)
# Remove case sensitivity
field_type = field_type.lower() if field_type is not None else "none"
vert_coord = vert_coord.lower() if vert_coord is not None else "none"
@ -583,7 +585,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, @@ -583,7 +585,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
interp_levels = np.asarray(interp_levels, np.float64)
# TODO: Check if field is staggered
if is_staggered(wrfin, field):
if is_staggered(_wrfin, field):
raise RuntimeError("Please unstagger field in the vertical")
# Check for valid coord
@ -606,23 +608,23 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, @@ -606,23 +608,23 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
# Extract vriables
#timeidx = -1 # Should this be an argument?
ncvars = extract_vars(wrfin, timeidx, ("PSFC", "QVAPOR", "F"),
ncvars = extract_vars(_wrfin, timeidx, ("PSFC", "QVAPOR", "F"),
method, squeeze, cache, meta=False, _key=_key)
sfp = ncvars["PSFC"] * ConversionFactors.PA_TO_HPA
qv = ncvars["QVAPOR"]
coriolis = ncvars["F"]
terht = get_terrain(wrfin, timeidx, units="m",
terht = get_terrain(_wrfin, timeidx, units="m",
method=method, squeeze=squeeze, cache=cache,
meta=False, _key=_key)
tk = get_temp(wrfin, timeidx, units="k",
tk = get_temp(_wrfin, timeidx, units="k",
method=method, squeeze=squeeze, cache=cache,
meta=False, _key=_key)
p = get_pressure(wrfin, timeidx, units="pa",
p = get_pressure(_wrfin, timeidx, units="pa",
method=method, squeeze=squeeze, cache=cache,
meta=False, _key=_key)
ght = get_height(wrfin, timeidx, msl=True, units="m",
ght = get_height(_wrfin, timeidx, msl=True, units="m",
method=method, squeeze=squeeze, cache=cache,
meta=False, _key=_key)
@ -639,7 +641,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, @@ -639,7 +641,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
vcord_array = np.exp(-ght/sclht)
elif vert_coord == "ght_agl":
ht_agl = get_height(wrfin, timeidx, msl=False, units="m",
ht_agl = get_height(_wrfin, timeidx, msl=False, units="m",
method=method, squeeze=squeeze, cache=cache,
meta=False, _key=_key)
@ -647,7 +649,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, @@ -647,7 +649,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
vcord_array = np.exp(-ht_agl/sclht)
elif vert_coord in ("theta", "th"):
t = get_theta(wrfin, timeidx, units="k",
t = get_theta(_wrfin, timeidx, units="k",
method=method, squeeze=squeeze, cache=cache,
meta=False, _key=_key)
@ -671,7 +673,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False, @@ -671,7 +673,7 @@ def vinterp(wrfin, field, vert_coord, interp_levels, extrapolate=False,
idir = 1
delta = 0.01
eth = get_eth(wrfin, timeidx, method=method, squeeze=squeeze,
eth = get_eth(_wrfin, timeidx, method=method, squeeze=squeeze,
cache=cache, meta=False, _key=_key)
p_hpa = p * ConversionFactors.PA_TO_HPA

27
src/wrf/latlonutils.py

@ -131,8 +131,9 @@ def _get_proj_params(wrfin, timeidx, stagger, method, squeeze, cache, _key): @@ -131,8 +131,9 @@ def _get_proj_params(wrfin, timeidx, stagger, method, squeeze, cache, _key):
"""
if timeidx < 0:
raise ValueError("'timeidx' must be greater than 0")
if timeidx is not None:
if timeidx < 0:
raise ValueError("'timeidx' must be greater than 0")
attrs = extract_global_attrs(wrfin, attrs=("MAP_PROJ", "TRUELAT1",
"TRUELAT2", "STAND_LON",
@ -383,12 +384,10 @@ def _ll_to_xy(latitude, longitude, wrfin=None, timeidx=0, @@ -383,12 +384,10 @@ def _ll_to_xy(latitude, longitude, wrfin=None, timeidx=0,
if ref_lat.size == 1:
outdim = [2, lats.size]
#outdim = [lats.size, 2]
extra_dims = [outdim[1]]
else:
# Moving domain will have moving ref_lats/ref_lons
outdim = [2, ref_lat.size, lats.size]
#outdim = [lats.size, ref_lat.size, 2]
extra_dims = outdim[1:]
result = np.empty(outdim, np.float64)
@ -402,11 +401,11 @@ def _ll_to_xy(latitude, longitude, wrfin=None, timeidx=0, @@ -402,11 +401,11 @@ def _ll_to_xy(latitude, longitude, wrfin=None, timeidx=0,
ref_lat_val = ref_lat[0]
ref_lon_val = ref_lon[0]
else:
ref_lat_val = ref_lat[left_idxs[-1]]
ref_lon_val = ref_lon[left_idxs[-1]]
ref_lat_val = ref_lat[left_idxs[-2]]
ref_lon_val = ref_lon[left_idxs[-2]]
lat = lats[left_idxs[0]]
lon = lons[left_idxs[0]]
lat = lats[left_idxs[-1]]
lon = lons[left_idxs[-1]]
xy = _lltoxy(map_proj, truelat1, truelat2, stdlon,
ref_lat_val, ref_lon_val, pole_lat, pole_lon,
@ -548,14 +547,10 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None, @@ -548,14 +547,10 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None,
raise ValueError("'x' and 'y' must be the same length")
if ref_lat.size == 1:
#outdim = [x_arr.size, 2]
#extra_dims = [outdim[0]]
outdim = [2, x_arr.size]
extra_dims = [outdim[1]]
else:
# Moving domain will have moving ref_lats/ref_lons
#outdim = [x_arr.size, ref_lat.size, 2]
#extra_dims = outdim[0:2]
outdim = [2, ref_lat.size, x_arr.size]
extra_dims = outdim[1:]
@ -570,11 +565,11 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None, @@ -570,11 +565,11 @@ def _xy_to_ll(x, y, wrfin=None, timeidx=0, stagger=None,
ref_lat_val = ref_lat[0]
ref_lon_val = ref_lon[0]
else:
ref_lat_val = ref_lat[left_idxs[-1]]
ref_lon_val = ref_lon[left_idxs[-1]]
ref_lat_val = ref_lat[left_idxs[-2]]
ref_lon_val = ref_lon[left_idxs[-2]]
x_val = x_arr[left_idxs[0]]
y_val = y_arr[left_idxs[0]]
x_val = x_arr[left_idxs[-1]]
y_val = y_arr[left_idxs[-1]]
ll = _xytoll(map_proj, truelat1, truelat2, stdlon, ref_lat_val,
ref_lon_val, pole_lat, pole_lon, known_x, known_y,

40
src/wrf/metadecorators.py

@ -10,7 +10,7 @@ import numpy.ma as ma @@ -10,7 +10,7 @@ import numpy.ma as ma
from .extension import _interpline
from .util import (extract_vars, either, from_args, arg_location,
is_coordvar, latlon_coordvars, to_np,
from_var, iter_left_indexes)
from_var, iter_left_indexes, is_mapping)
from .coordpair import CoordPair
from .py3compat import viewkeys, viewitems, py3range, ucode
from .interputils import get_xy_z_params, get_xy, to_xy_coords
@ -586,8 +586,10 @@ def set_latlon_metadata(xy=False): @@ -586,8 +586,10 @@ def set_latlon_metadata(xy=False):
@wrapt.decorator
def func_wrapper(wrapped, instance, args, kwargs):
do_meta = from_args(wrapped, ("meta",), *args, **kwargs)["meta"]
argvars = from_args(wrapped, ("wrfin", "meta"), *args, **kwargs)
# If it's a mapping, then this is handled as a special case in g_latlon
do_meta = (not is_mapping(argvars["wrfin"]) and argvars["meta"])
if do_meta is None:
do_meta = True
@ -1985,7 +1987,7 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): @@ -1985,7 +1987,7 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"):
p = argvals[copyarg]
missing = argvals["missing"]
# Note: cape_3d supports using only a single column of data
# Note: 2D/3D cape supports using only a single column of data
is1d = p.ndim == 1
# Need to squeeze the right dimensions for 1D cape
@ -1997,31 +1999,34 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): @@ -1997,31 +1999,34 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"):
outattrs = OrderedDict()
if is2d:
outname = "cape_2d"
if is1d:
outname = "cape_2d"
else:
outname = "cape_2d_column"
outattrs["description"] = "mcape ; mcin ; lcl ; lfc"
outattrs["units"] = "J kg-1 ; J kg-1 ; m ; m"
outattrs["MemoryOrder"] = "XY"
outattrs["MemoryOrder"] = ""
else:
if not is1d:
outname = "cape_3d"
outattrs["description"] = "cape; cin"
outattrs["units"] = "J kg-1 ; J kg-1"
outattrs["MemoryOrder"] = "XYZ"
else:
outname = "cape_column"
outattrs["description"] = "cape; cin"
outattrs["units"] = "J kg-1 ; J kg-1"
outname = "cape_3d_column"
outattrs["MemoryOrder"] = "Z"
outattrs["description"] = "cape; cin"
outattrs["units"] = "J kg-1 ; J kg-1"
if isinstance(p, DataArray):
if is2d:
# Right dims
outdims[-2:] = p.dims[-2:]
# Left dims
outdims[1:-2] = p.dims[0:-3]
if not is1d:
# Right dims
outdims[-2:] = p.dims[-2:]
# Left dims
outdims[1:-2] = p.dims[0:-3]
else:
if not is1d:
# Right dims
@ -2030,12 +2035,13 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): @@ -2030,12 +2035,13 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"):
outdims[1:-3] = p.dims[0:-3]
else:
outdims[1] = p.dims[0]
outcoords = {}
# Left-most is always cape_cin or cape_cin_lcl_lfc
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"]

5
src/wrf/routines.py

@ -112,7 +112,7 @@ _VALID_KARGS = {"cape2d" : ["missing"], @@ -112,7 +112,7 @@ _VALID_KARGS = {"cape2d" : ["missing"],
"wspd_wdir10" : ["units"],
"uvmet_wspd_wdir" : ["units"],
"uvmet10_wspd_wdir" : ["units"],
"ctt" : [],
"ctt" : ["fill_nocloud", "missing", "opt_thresh", "units"],
"cloudfrac" : ["vert_type", "low_thresh",
"mid_thresh", "high_thresh"],
"geopt_stag" : [],
@ -138,7 +138,8 @@ _ALIASES = {"cape_2d" : "cape2d", @@ -138,7 +138,8 @@ _ALIASES = {"cape_2d" : "cape2d",
"td2" : "dp2m",
"cfrac" : "cloudfrac",
"wspd_wdir_uvmet" : "uvmet_wspd_wdir",
"wspd_wdir_uvmet10" : "uvmet10_wspd_wdir"
"wspd_wdir_uvmet10" : "uvmet10_wspd_wdir",
"th" : "theta"
}
class ArgumentError(Exception):

45
src/wrf/specialdec.py

@ -229,7 +229,14 @@ def cape_left_iter(alg_dtype=np.float64): @@ -229,7 +229,14 @@ def cape_left_iter(alg_dtype=np.float64):
ter_follow = args[8]
is2d = i3dflag == 0
is1d = np.isscalar(sfp)
# Note: This should still work with DataArrays
is1d = np.isscalar(sfp) or np.size(sfp) == 1
# Make sure sfp and terrain are regular floats for 1D case
# This should also work with DataArrays
if is1d:
ter = float(ter)
sfp = float(sfp)
orig_dtype = p_hpa.dtype
@ -542,39 +549,27 @@ def check_cape_args(): @@ -542,39 +549,27 @@ def check_cape_args():
ter_follow = args[8]
is2d = False if i3dflag != 0 else True
is1d = ((np.isscalar(sfp) or np.size(sfp) == 1) or
(np.isscalar(ter) or np.size(ter) == 1))
if not (p_hpa.shape == tk.shape == qv.shape == ht.shape):
raise ValueError("arguments 0, 1, 2, 3 must be the same shape")
# 2D CAPE does not allow for scalars
if is2d:
if np.isscalar(ter) or np.isscalar(sfp):
raise ValueError("arguments 4 and 5 cannot be scalars with "
"cape_2d routine")
if not is1d:
if ter.ndim != p_hpa.ndim-1 or sfp.ndim != p_hpa.ndim-1:
raise ValueError("arguments 4 and 5 must have "
"{} dimensions".format(p_hpa.ndim-1))
# 3D cape is allowed to be just a vertical column with scalars
# for terrain and psfc_hpa.
else:
if np.isscalar(ter) and np.isscalar(sfp):
if p_hpa.ndim != 1:
raise ValueError("arguments 0-3 "
"must be 1-dimensional when "
"arguments 4 and 5 are scalars")
if is2d:
raise ValueError("argument 7 must be 0 when using 1D data")
else:
if ((np.isscalar(ter) and not np.isscalar(sfp)) or
(not np.isscalar(ter) and np.isscalar(sfp))):
raise ValueError("arguments 4 and 5 must both be scalars")
else:
if ter.ndim != p_hpa.ndim-1 or sfp.ndim != p_hpa.ndim-1:
raise ValueError("arguments 4 and 5 "
"must have {} dimensions".format(
p_hpa.ndim-1))
if np.size(ter) != np.size(sfp):
raise ValueError("arguments 4 and 5 must both be scalars or "
"both be arrays")
# Only need to test p_hpa since we assured args 0-3 have same ndim
if p_hpa.ndim != 1:
raise ValueError("arguments 0-3 "
"must be 1-dimensional when "
"arguments 4 and 5 are scalars")
return wrapped(*args, **kwargs)

10
src/wrf/util.py

@ -973,10 +973,12 @@ def _combine_dict(wrfdict, varname, timeidx, method, meta, _key): @@ -973,10 +973,12 @@ def _combine_dict(wrfdict, varname, timeidx, method, meta, _key):
is_moving = is_moving_domain(wrfdict, varname, _key=_key)
_cache_key = _key[first_key] if _key is not None else None
first_array = _extract_var(wrfdict[first_key], varname,
timeidx, is_moving=is_moving, method=method,
squeeze=False, cache=None, meta=meta,
_key=_key[first_key])
_key=_cache_key)
# Create the output data numpy array based on the first array
outdims = [numkeys]
@ -992,10 +994,11 @@ def _combine_dict(wrfdict, varname, timeidx, method, meta, _key): @@ -992,10 +994,11 @@ def _combine_dict(wrfdict, varname, timeidx, method, meta, _key):
break
else:
keynames.append(key)
_cache_key = _key[key] if _key is not None else None
vardata = _extract_var(wrfdict[key], varname, timeidx,
is_moving=is_moving, method=method,
squeeze=False, cache=None, meta=meta,
_key=_key[key])
_key=_cache_key)
if outdata.shape[1:] != vardata.shape:
raise ValueError("data sequences must have the "
@ -3070,7 +3073,8 @@ def get_id(obj, prefix=''): @@ -3070,7 +3073,8 @@ def get_id(obj, prefix=''):
# For sequences, the hashing string will be the list ID and the
# path for the first file in the sequence
if not is_mapping(obj):
_next = next(iter(obj))
_obj = get_iterable(obj)
_next = next(iter(_obj))
return get_id(_next, prefix + str(id(obj)))
# For each key in the mapping, recursively call get_id until

2
src/wrf/version.py

@ -1,2 +1,2 @@ @@ -1,2 +1,2 @@
__version__ = "1.1.2"
__version__ = "1.1.3"

BIN
test/ci_tests/ci_result_file.nc

Binary file not shown.

2
test/ci_tests/make_test_file.py

@ -12,7 +12,7 @@ from wrf import (getvar, interplevel, interpline, vertcross, vinterp, py2round, @@ -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",

65
test/comp_utest.py

@ -597,9 +597,67 @@ def test_cape3d_1d(wrfnc): @@ -597,9 +597,67 @@ def test_cape3d_1d(wrfnc):
nt.assert_allclose(to_np(result), to_np(ref))
return func
def test_cape2d_1d(wrfnc):
def func(self):
varnames = ("T", "P", "PB", "QVAPOR", "PH", "PHB", "HGT", "PSFC")
ncvars = extract_vars(wrfnc, 0, varnames, method="cat", squeeze=True,
cache=None, meta=True)
t = ncvars["T"]
p = ncvars["P"]
pb = ncvars["PB"]
qv = ncvars["QVAPOR"]
ph = ncvars["PH"]
phb = ncvars["PHB"]
ter = ncvars["HGT"]
psfc = ncvars["PSFC"]
col_idxs = (slice(None), t.shape[-2]//2, t.shape[-1]//2)
t = t[col_idxs]
p = p[col_idxs]
pb = pb[col_idxs]
qv = qv[col_idxs]
ph = ph[col_idxs]
phb = phb[col_idxs]
ter = float(ter[col_idxs[1:]])
psfc = float(psfc[col_idxs[1:]])
full_t = t + Constants.T_BASE
full_p = p + pb
tkel = tk(full_p, full_t, meta=False)
geopt = ph + phb
geopt_unstag = destagger(geopt, -1)
z = geopt_unstag/Constants.G
# Convert pressure to hPa
p_hpa = ConversionFactors.PA_TO_HPA * full_p
psfc_hpa = ConversionFactors.PA_TO_HPA * psfc
i3dflag = 0
ter_follow = 1
result = cape_2d(p_hpa, tkel, qv, z, ter, psfc_hpa, ter_follow)
ref = getvar(wrfnc, "cape_2d")
ref = ref[(slice(None),) + col_idxs[1:]]
nt.assert_allclose(to_np(result), to_np(ref))
return func
if __name__ == "__main__":
from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule,
omp_set_dynamic, omp_get_num_procs, OMP_SCHED_STATIC)
omp_set_num_threads(omp_get_num_procs()//2)
omp_set_schedule(OMP_SCHED_STATIC, 0)
omp_set_dynamic(False)
varnames = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz",
"geopt", "helicity", "lat", "lon", "omg", "p", "pressure",
"pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc",
@ -609,14 +667,14 @@ if __name__ == "__main__": @@ -609,14 +667,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)):
@ -639,6 +697,9 @@ if __name__ == "__main__": @@ -639,6 +697,9 @@ if __name__ == "__main__":
func = test_cape3d_1d(wrfnc)
setattr(WRFVarsTest, "test_cape3d_1d", func)
func = test_cape2d_1d(wrfnc)
setattr(WRFVarsTest, "test_cape2d_1d", func)
ut.main()

50
test/ctt_test.py

@ -0,0 +1,50 @@ @@ -0,0 +1,50 @@
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature
from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords
# Open the NetCDF file
ncfile = Dataset("/Users/ladwig/Documents/wrf_files/problem_files/cfrac_bug/wrfout_d02_1987-10-01_00:00:00")
# Get the sea level pressure
ctt = getvar(ncfile, "ctt")
# Get the latitude and longitude points
lats, lons = latlon_coords(ctt)
# Get the cartopy mapping object
cart_proj = get_cartopy(ctt)
# Create a figure
fig = plt.figure(figsize=(12,9))
# Set the GeoAxes to the projection used by WRF
ax = plt.axes(projection=cart_proj)
# Download and add the states and coastlines
states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
name='admin_1_states_provinces_shp')
ax.add_feature(states, linewidth=.5)
ax.coastlines('50m', linewidth=0.8)
# Make the contour outlines and filled contours for the smoothed sea level pressure.
plt.contour(to_np(lons), to_np(lats), to_np(ctt), 10, colors="black",
transform=crs.PlateCarree())
plt.contourf(to_np(lons), to_np(lats), to_np(ctt), 10, transform=crs.PlateCarree(),
cmap=get_cmap("jet"))
# Add a color bar
plt.colorbar(ax=ax, shrink=.62)
# Set the map limits. Not really necessary, but used for demonstration.
ax.set_xlim(cartopy_xlim(ctt))
ax.set_ylim(cartopy_ylim(ctt))
# Add the gridlines
ax.gridlines(color="black", linestyle="dotted")
plt.title("Cloud Top Temperature")
plt.show()

2
test/ipynb/Doc_Examples.ipynb

@ -1275,7 +1275,7 @@ @@ -1275,7 +1275,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.13"
"version": "2.7.14"
}
},
"nbformat": 4,

214
test/ipynb/Test_CTT.ipynb

@ -0,0 +1,214 @@ @@ -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.14"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

214
test/ipynb/Test_CTT_Prod.ipynb

@ -0,0 +1,214 @@ @@ -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
}

37
test/ipynb/WRF_python_demo.ipynb

@ -251,7 +251,7 @@ @@ -251,7 +251,7 @@
" return self\n",
" \n",
" def next(self):\n",
" if self._i > self._total:\n",
" if self._i >= self._total:\n",
" raise StopIteration\n",
" else:\n",
" val = self.ncfile[self._i]\n",
@ -653,7 +653,7 @@ @@ -653,7 +653,7 @@
"metadata": {},
"outputs": [],
"source": [
"from wrf.latlon import xy_to_ll, ll_to_xy \n",
"from wrf import xy_to_ll, ll_to_xy \n",
"\n",
"a = xy_to_ll(ncfile, 400, 200)\n",
"a1 = ll_to_xy(ncfile, a[0], a[1])\n",
@ -1134,7 +1134,7 @@ @@ -1134,7 +1134,7 @@
"metadata": {},
"outputs": [],
"source": [
"from wrf.latlon import xy_to_ll, ll_to_xy \n",
"from wrf import xy_to_ll, ll_to_xy \n",
"\n",
"a = xy_to_ll(ncfiles, 400, 200)\n",
"a = xy_to_ll(ncfiles, [400,105], [200,205])\n",
@ -1196,6 +1196,27 @@ @@ -1196,6 +1196,27 @@
"print (\"\\n\")\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
@ -1206,21 +1227,21 @@ @@ -1206,21 +1227,21 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"display_name": "Python 3",
"language": "python",
"name": "python2"
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.13"
"pygments_lexer": "ipython3",
"version": "3.6.4"
}
},
"nbformat": 4,

181
test/test_inputs.py

@ -0,0 +1,181 @@ @@ -0,0 +1,181 @@
import unittest as ut
import numpy.testing as nt
import numpy as np
import numpy.ma as ma
import os
import sys
import subprocess
from wrf import (getvar, interplevel, interpline, vertcross, vinterp,
disable_xarray, xarray_enabled, to_np,
xy_to_ll, ll_to_xy, xy_to_ll_proj, ll_to_xy_proj,
extract_global_attrs, viewitems, CoordPair, ll_points)
from wrf.util import is_multi_file
IN_DIR = "/Users/ladwig/Documents/wrf_files/wrf_vortex_multi"
TEST_FILES = [os.path.join(IN_DIR, "wrfout_d02_2005-08-28_00:00:00"),
os.path.join(IN_DIR, "wrfout_d02_2005-08-28_12:00:00"),
os.path.join(IN_DIR, "wrfout_d02_2005-08-29_00:00:00")]
def wrfin_gen(wrf_in):
for x in wrf_in:
yield x
class wrf_in_iter_class(object):
def __init__(self, wrf_in):
self._wrf_in = wrf_in
self._total = len(wrf_in)
self._i = 0
def __iter__(self):
return self
def next(self):
if self._i >= self._total:
raise StopIteration
else:
result = self._wrf_in[self._i]
self._i += 1
return result
# Python 3
def __next__(self):
return self.next()
def make_test(varname, wrf_in):
def test(self):
x = getvar(wrf_in, varname, 0)
xa = getvar(wrf_in, varname, None)
return test
def make_interp_test(varname, wrf_in):
def test(self):
# Only testing vinterp since other interpolation just use variables
if (varname == "vinterp"):
for timeidx in (0, None):
eth = getvar(wrf_in, "eth", timeidx=timeidx)
interp_levels = [850,500,5]
field = vinterp(wrf_in,
field=eth,
vert_coord="pressure",
interp_levels=interp_levels,
extrapolate=True,
field_type="theta-e",
timeidx=timeidx,
log_p=True)
else:
pass
return test
def make_latlon_test(testid, wrf_in):
def test(self):
if testid == "xy_out":
# Lats/Lons taken from NCL script, just hard-coding for now
lats = [-55, -60, -65]
lons = [25, 30, 35]
xy = ll_to_xy(wrf_in, lats, lons, timeidx=0)
xy = ll_to_xy(wrf_in, lats, lons, timeidx=None)
else:
i_s = np.asarray([10, 100, 150], int)
j_s = np.asarray([10, 100, 150], int)
ll = xy_to_ll(wrf_in, i_s, j_s, timeidx=0)
ll = xy_to_ll(wrf_in, i_s, j_s, timeidx=None)
return test
class WRFVarsTest(ut.TestCase):
longMessage = True
class WRFInterpTest(ut.TestCase):
longMessage = True
class WRFLatLonTest(ut.TestCase):
longMessage = True
if __name__ == "__main__":
from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule,
omp_set_dynamic, omp_get_num_procs, OMP_SCHED_STATIC)
omp_set_num_threads(omp_get_num_procs()-1)
omp_set_schedule(OMP_SCHED_STATIC, 0)
omp_set_dynamic(False)
ignore_vars = [] # Not testable yet
wrf_vars = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz",
"geopt", "helicity", "lat", "lon", "omg", "p", "pressure",
"pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc",
"theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va",
"wa", "uvmet10", "uvmet", "z", "cfrac", "zstag", "geopt_stag"]
interp_methods = ["interplevel", "vertcross", "interpline", "vinterp"]
latlon_tests = ["xy_out", "ll_out"]
for nc_lib in ("netcdf4", "pynio", "scipy"):
if nc_lib == "netcdf4":
try:
from netCDF4 import Dataset
except ImportError:
print ("netcdf4-python not installed")
continue
else:
test_in = [Dataset(x) for x in TEST_FILES]
elif nc_lib == "pynio":
try:
from Nio import open_file
except ImportError:
print ("PyNIO not installed")
continue
else:
test_in = [open_file(x +".nc", "r") for x in TEST_FILES]
elif nc_lib == "scipy":
try:
from scipy.io.netcdf import netcdf_file
except ImportError:
print ("scipy.io.netcdf not installed")
else:
test_in = [netcdf_file(x, mmap=False) for x in TEST_FILES]
input0 = test_in[0]
input1 = test_in
input2 = tuple(input1)
input3 = wrfin_gen(test_in)
input4 = wrf_in_iter_class(test_in)
input5 = {"A" : input1,
"B" : input2}
input6 = {"A" : {"AA" : input1},
"B" : {"BB" : input2}}
for i,input in enumerate((input0, input1, input2,
input3, input4, input5, input6)):
for var in wrf_vars:
if var in ignore_vars:
continue
test_func1 = make_test(var, input)
setattr(WRFVarsTest, "test_{0}_input{1}_{2}".format(nc_lib,
i, var),
test_func1)
for method in interp_methods:
test_interp_func1 = make_interp_test(method, input)
setattr(WRFInterpTest,
"test_{0}_input{1}_{2}".format(nc_lib,
i, method),
test_interp_func1)
for testid in latlon_tests:
test_ll_func = make_latlon_test(testid, input)
test_name = "test_{0}_input{1}_{2}".format(nc_lib, i, testid)
setattr(WRFLatLonTest, test_name, test_ll_func)
ut.main()

2
test/utests.py

@ -682,7 +682,7 @@ class WRFLatLonTest(ut.TestCase): @@ -682,7 +682,7 @@ class WRFLatLonTest(ut.TestCase):
if __name__ == "__main__":
from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule,
omp_set_dynamic, omp_get_num_procs, OMP_SCHED_STATIC)
omp_set_num_threads(omp_get_num_procs()-1)
omp_set_num_threads(omp_get_num_procs()//2)
omp_set_schedule(OMP_SCHED_STATIC, 0)
omp_set_dynamic(False)

Loading…
Cancel
Save