Merge branch 'bugfix/proj_globe' into develop

src/wrf/projection.py

@@ -309,7 +309,8 @@ class WrfProj(object):
         return (None if not cartopy_enabled() 
                 else crs.Globe(ellipse=None,
                                semimajor_axis=Constants.WRF_EARTH_RADIUS,
-                               semiminor_axis=Constants.WRF_EARTH_RADIUS))
+                               semiminor_axis=Constants.WRF_EARTH_RADIUS,
+                               nadgrids="@null"))
      
     def cartopy_xlim(self, geobounds):
         """Return the x extents in projected coordinates for cartopy.
@@ -610,7 +611,7 @@ class LambertConformal(WrfProj):
                     else self.truelat2)
         
         _proj4 = ("+proj=lcc +units=meters +a={} +b={} +lat_1={} "
-                       "+lat_2={} +lat_0={} +lon_0={}".format(
+                  "+lat_2={} +lat_0={} +lon_0={} +nadgrids=@null".format(
                                             Constants.WRF_EARTH_RADIUS,
                                             Constants.WRF_EARTH_RADIUS,
                                             self.truelat1,
@@ -731,7 +732,7 @@ class Mercator(WrfProj):
     def _proj4(self):
         
         _proj4 = ("+proj=merc +units=meters +a={} +b={} "
-                       "+lon_0={} +lat_ts={}".format(
+                  "+lon_0={} +lat_ts={} +nadgrids=@null".format(
                                             Constants.WRF_EARTH_RADIUS,
                                             Constants.WRF_EARTH_RADIUS,
                                             self._stand_lon,
@@ -845,7 +846,7 @@ class PolarStereographic(WrfProj):
     
     def _proj4(self):
         _proj4 = ("+proj=stere +units=meters +a={} +b={} "
-                       "+lat0={} +lon_0={} +lat_ts={}".format(
+                  "+lat0={} +lon_0={} +lat_ts={} +nadgrids=@null".format(
                                             Constants.WRF_EARTH_RADIUS,
                                             Constants.WRF_EARTH_RADIUS,
                                             self._hemi,
@@ -951,7 +952,7 @@ class LatLon(WrfProj):
     
     def _proj4(self):
         _proj4 = ("+proj=eqc +units=meters +a={} +b={} "
-                       "+lon_0={}".format(Constants.WRF_EARTH_RADIUS,
+                  "+lon_0={} +nadgrids=@null".format(Constants.WRF_EARTH_RADIUS,
                                           Constants.WRF_EARTH_RADIUS,
                                           self.stand_lon))
         return _proj4

test/plot_tests/plot_all.ipynb

@@ -0,0 +1,325 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from __future__ import print_function\n",
+    "\n",
+    "# This jupyter notebook command inserts matplotlib graphics in \n",
+    "# to the workbook\n",
+    "%matplotlib inline\n",
+    "\n",
+    "# Modify these to point to your own files\n",
+    "WRF_DIRECTORY = \"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/moving_nest\"\n",
+    "\n",
+    "WRF_FILES = [\"wrfout_d02_2005-08-28_00:00:00\",\n",
+    "             \"wrfout_d02_2005-08-28_12:00:00\",\n",
+    "             \"wrfout_d02_2005-08-29_00:00:00\"]\n",
+    "\n",
+    "\n",
+    "# Do not modify the code below this line\n",
+    "#------------------------------------------------------\n",
+    "# Turn off annoying warnings\n",
+    "import warnings\n",
+    "warnings.filterwarnings('ignore')\n",
+    "\n",
+    "# Make sure the environment is good\n",
+    "import numpy\n",
+    "import cartopy\n",
+    "import matplotlib\n",
+    "from netCDF4 import Dataset\n",
+    "from xarray import DataArray\n",
+    "from wrf import (getvar, interplevel, vertcross, \n",
+    "                 vinterp, ALL_TIMES)\n",
+    "import os\n",
+    "\n",
+    "_WRF_FILES = [os.path.abspath(os.path.expanduser(\n",
+    "    os.path.join(WRF_DIRECTORY, f))) for f in WRF_FILES]\n",
+    "\n",
+    "# Check that the WRF files exist\n",
+    "try:\n",
+    "    for f in _WRF_FILES:\n",
+    "        if not os.path.exists(f):\n",
+    "            raise ValueError(\"{} does not exist. \"\n",
+    "                \"Check for typos or incorrect directory.\".format(f))\n",
+    "except ValueError:\n",
+    "    # If the directory ended up in the zip file, then \n",
+    "    # another 'wrf_tutorial_data' directory might be \n",
+    "    # used.\n",
+    "    WRF_DIRECTORY = os.path.join(WRF_DIRECTORY, \"wrf_tutorial_data\")\n",
+    "    _WRF_FILES = [os.path.abspath(os.path.expanduser(\n",
+    "        os.path.join(WRF_DIRECTORY, f))) for f in WRF_FILES]\n",
+    "    for f in _WRF_FILES:\n",
+    "        if not os.path.exists(f):\n",
+    "            raise\n",
+    "\n",
+    "            \n",
+    "# Create functions so that the WRF files only need\n",
+    "# to be specified using the WRF_FILES global above\n",
+    "def single_wrf_file():\n",
+    "    global _WRF_FILES\n",
+    "    return _WRF_FILES[0]\n",
+    "\n",
+    "def multiple_wrf_files():\n",
+    "    global _WRF_FILES\n",
+    "    return _WRF_FILES\n",
+    "\n",
+    "def save_fig(diagname):\n",
+    "    f = single_wrf_file()\n",
+    "    if f.find(\"_d01_\") > 0:\n",
+    "        fout = \"{}.png\".format(os.path.join(os.path.abspath(\".\"), \"d01\", diagname))\n",
+    "    else:\n",
+    "        fout = \"{}.png\".format(os.path.join(os.path.abspath(\".\"), \"d02\", diagname))\n",
+    "    \n",
+    "    matplotlib.pyplot.savefig(fout)\n",
+    "        \n",
+    "\n",
+    "print(\"All tests passed!\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def get_kwargs(diagname):\n",
+    "    kwargs = {}\n",
+    "    if diagname == \"ctt\":\n",
+    "        kwargs = {\"fill_nocloud\" : True}\n",
+    "    \n",
+    "    return kwargs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def contour_levels(diagname, diag, numlevels=15):\n",
+    "    levels = numlevels\n",
+    "    extend = \"neither\"\n",
+    "    if diagname == \"ter\":\n",
+    "        levels = numpy.arange(10.,4000.,250.)\n",
+    "    elif diagname == \"avo\":\n",
+    "        levels = numpy.arange(10.,75.,5.)\n",
+    "        extend=\"max\"\n",
+    "    elif diagname == \"eth\":\n",
+    "        levels = numpy.arange(270.,400.,10.)\n",
+    "    elif diagname == \"cape_2d\":\n",
+    "        levels = numpy.arange(200.,4000.,250.)\n",
+    "        extend = \"max\"\n",
+    "    elif diagname == \"cape_3d\":\n",
+    "        levels = numpy.arange(200.,4000.,250.)\n",
+    "        extend = \"max\"\n",
+    "    elif diagname == \"ctt\":\n",
+    "        # Note: The MP scheme doesn't produce cloud water so this\n",
+    "        # is just surface temperature\n",
+    "        levels = numpy.arange(-100, 20., 5.0)\n",
+    "        extend = \"both\"\n",
+    "    elif diagname == \"dbz\":\n",
+    "        #pass\n",
+    "        levels = numpy.arange(15.,75.,5.)\n",
+    "        extend = \"max\"\n",
+    "    elif diagname == \"mdbz\":\n",
+    "        levels = numpy.arange(15.,75.,5.)\n",
+    "        extend = \"max\"\n",
+    "    elif diagname == \"geopt\":\n",
+    "        pass\n",
+    "    elif diagname == \"helicity\":\n",
+    "        levels = numpy.arange(100,500,25)\n",
+    "        extend=\"max\"\n",
+    "    elif diagname == \"lat\":\n",
+    "        pass\n",
+    "    elif diagname == \"lon\":\n",
+    "        pass\n",
+    "    elif diagname == \"omg\":\n",
+    "        pass\n",
+    "    elif diagname == \"p\":\n",
+    "        levels = numpy.arange(95000.,103000.,250.)\n",
+    "        extend = \"min\"\n",
+    "    elif diagname == \"pressure\":\n",
+    "        levels = numpy.arange(950.,1030.,2.5)\n",
+    "        extend = \"min\"\n",
+    "    elif diagname == \"pvo\":\n",
+    "        levels = numpy.arange(.5,5.,.25)\n",
+    "    elif diagname == \"pw\":\n",
+    "        #pass\n",
+    "        levels = numpy.arange(0.1,100,10)\n",
+    "    elif diagname == \"rh2\":\n",
+    "        levels = numpy.arange(50,101,5)\n",
+    "    elif diagname == \"rh\":\n",
+    "        levels = numpy.arange(50.,101.,5)\n",
+    "    elif diagname == \"slp\":\n",
+    "        levels = numpy.arange(950.,1030.,2.5)\n",
+    "        extend = \"min\"\n",
+    "    elif diagname == \"td2\":\n",
+    "        levels = numpy.arange(10,40,5)\n",
+    "        extend = \"max\"\n",
+    "    elif diagname == \"td\":\n",
+    "        levels = numpy.arange(10,40,5)\n",
+    "        extend = \"max\"\n",
+    "    elif diagname == \"tc\":\n",
+    "        levels = numpy.arange(10,40,5)\n",
+    "        extend = \"max\"\n",
+    "    elif diagname == \"theta\":\n",
+    "        levels = numpy.arange(270.,350.,2.5)\n",
+    "    elif diagname == \"tk\":\n",
+    "        levels = numpy.arange(270.,350.,5.)\n",
+    "    elif diagname == \"tv\":\n",
+    "        levels = numpy.arange(270.,350.,5.)\n",
+    "        #levels = numpy.arange(10,40,5)\n",
+    "    elif diagname == \"twb\":\n",
+    "        levels = numpy.arange(270.,350.,5.)    \n",
+    "    elif diagname == \"updraft_helicity\":\n",
+    "        pass\n",
+    "        #levels = numpy.arange(1,100,5)\n",
+    "    elif diagname == \"ua\":\n",
+    "        levels = numpy.arange(-40,40,5)\n",
+    "        extend = \"both\"\n",
+    "    elif diagname == \"va\":\n",
+    "        levels = numpy.arange(-40,40,5)\n",
+    "        extend = \"both\"\n",
+    "    elif diagname == \"wa\":\n",
+    "        pass\n",
+    "        #levels = numpy.arange(0.1,40,5)\n",
+    "    elif diagname == \"uvmet10\":\n",
+    "        levels = numpy.arange(-40,40,5)\n",
+    "        extend = \"both\"\n",
+    "    elif diagname == \"uvmet\":\n",
+    "        levels = numpy.arange(-40,40,5)\n",
+    "        extend = \"both\"\n",
+    "    elif diagname == \"z\":\n",
+    "        levels = numpy.arange(0,200,10)\n",
+    "        extend=\"max\"\n",
+    "    elif diagname == \"cfrac\":\n",
+    "        levels = numpy.arange(0.0,1.1,.2)\n",
+    "    elif diagname == \"wspd_wdir\":\n",
+    "        levels = numpy.arange(-40,40,5)\n",
+    "        extend = \"both\"\n",
+    "    elif diagname == \"wspd_wdir10\":\n",
+    "        levels = numpy.arange(-40,40,5)\n",
+    "        extend = \"both\"\n",
+    "    elif diagname == \"uvmet_wspd_wdir\":\n",
+    "        levels = numpy.arange(-40,40,5)\n",
+    "        extend = \"both\"\n",
+    "    elif diagname == \"uvmet10_wspd_wdir\":\n",
+    "        levels = numpy.arange(-40,40,5)\n",
+    "        extend = \"both\"\n",
+    "    \n",
+    "    return levels, extend"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "import numpy\n",
+    "from matplotlib import pyplot\n",
+    "from matplotlib.cm import get_cmap\n",
+    "from cartopy import crs\n",
+    "from cartopy.feature import NaturalEarthFeature\n",
+    "from netCDF4 import Dataset\n",
+    "from wrf import getvar, to_np, get_cartopy, latlon_coords, cartopy_xlim, cartopy_ylim, Constants\n",
+    "\n",
+    "file_path = single_wrf_file()\n",
+    "wrf_file = Dataset(file_path)\n",
+    "\n",
+    "for diagname in (\"ter\",\"avo\", \"eth\", \"cape_2d\", \"cape_3d\", \"ctt\", \"dbz\", \"mdbz\", \n",
+    "            \"geopt\", \"helicity\", \"lat\", \"lon\", \"omg\", \"p\", \"pressure\", \n",
+    "            \"pvo\", \"pw\", \"rh2\", \"rh\", \"slp\", \"td2\", \"td\", \"tc\",\n",
+    "            \"theta\", \"tk\", \"tv\", \"twb\", \"updraft_helicity\", \"ua\", \"va\", \n",
+    "            \"wa\", \"uvmet10\", \"uvmet\", \"z\", \"wspd_wdir\", \"wspd_wdir10\",\n",
+    "            \"uvmet_wspd_wdir\", \"uvmet10_wspd_wdir\", \"cfrac\"):\n",
+    "    \n",
+    "    # Get the terrain height\n",
+    "    print(diagname)\n",
+    "    kwargs = get_kwargs(diagname)\n",
+    "    diag = getvar(wrf_file, diagname, timeidx=3, **kwargs)\n",
+    "\n",
+    "    if diag.ndim == 3:\n",
+    "        diag = diag[0,:]\n",
+    "    elif diag.ndim == 4:\n",
+    "        diag = diag[0,0,:]\n",
+    "\n",
+    "    # Get the cartopy object and the lat,lon coords\n",
+    "    cart_proj = get_cartopy(diag)\n",
+    "    lats, lons = latlon_coords(diag)\n",
+    "    \n",
+    "\n",
+    "    # Create a figure and get the GetAxes object\n",
+    "    fig = pyplot.figure(figsize=(10, 7.5))\n",
+    "    ax = pyplot.axes(projection=cart_proj)\n",
+    "\n",
+    "    # Download and add the states and coastlines\n",
+    "    # See the cartopy documentation for more on this.\n",
+    "    states = NaturalEarthFeature(category='cultural', \n",
+    "                                 scale='50m', \n",
+    "                                 facecolor='none',\n",
+    "                                 name='admin_1_states_provinces_shp')\n",
+    "\n",
+    "    ax.add_feature(states, linewidth=.5, edgecolor='black', zorder=3)\n",
+    "    ax.coastlines('50m', linewidth=.8, color='black', zorder=4)\n",
+    "\n",
+    "    # Set the contour levels\n",
+    "    levels, extend = contour_levels(diagname, diag, numlevels=15)\n",
+    "\n",
+    "    # Make the contour lines and fill them.\n",
+    "    pyplot.contour(to_np(lons), to_np(lats), \n",
+    "                   to_np(diag), levels=levels, \n",
+    "                   colors=\"black\",\n",
+    "                   transform=crs.PlateCarree())\n",
+    "    pyplot.contourf(to_np(lons), to_np(lats), \n",
+    "                    to_np(diag), levels=levels,\n",
+    "                    transform=crs.PlateCarree(),\n",
+    "                    extend=extend,\n",
+    "                    cmap=get_cmap(\"jet\"))\n",
+    "\n",
+    "    ax.set_xlim(cartopy_xlim(diag))\n",
+    "    ax.set_ylim(cartopy_ylim(diag))\n",
+    "             \n",
+    "    # Add a color bar. The shrink often needs to be set \n",
+    "    # by trial and error.\n",
+    "    cb = pyplot.colorbar(ax=ax, shrink=.99)\n",
+    "    \n",
+    "    pyplot.title(diagname)\n",
+    "    \n",
+    "    # Uncomment this to build the sample images\n",
+    "    #save_fig(diagname)\n",
+    "    \n",
+    "    pyplot.show()\n",
+    "\n",
+    "\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}