Merge pull request #24 from jdavid/master

plotting default map_region is 'default'

examples/small-plot.py

@@ -10,22 +10,17 @@
     http://niflheim.nilu.no/~burkhart/sharing/pflexpart_testdata.tgz
 """
 
-import os
-import sys
 
-import matplotlib as mpl
-import matplotlib.pyplot as plt
-import numpy as np
-from mpl_toolkits.basemap import addcyclic
+from argparse import ArgumentParser
+import os
 
 import reflexible as rf
-import reflexible.mapping as mp
-from reflexible.plotting import plot_sensitivity, plot_totalcolumn, plot_at_level
+from reflexible.plotting import plot_totalcolumn, plot_at_level
 
 
-def plot_backward(SOURCE_FILE, OUTPUT_DIR):
+def plot_backward(SOURCE_FILE, OUTPUT_DIR, region, nested=False, label=''):
     #read the header of a FLEXPART output
-    H = rf.Header(SOURCE_FILE, nested=True)
+    H = rf.Header(SOURCE_FILE, nested=nested)
 
     # # Get the trajectories, these are overlayed when we plot the
     # # sensitivities below.
@@ -37,23 +32,18 @@ def plot_backward(SOURCE_FILE, OUTPUT_DIR):
     FP = None #FP Empty figure object
 
     # iterate over every species and every timestep (k)
-
-    # The region depends on the file, this should be a parameter.
-    map_region = 'north_sea'
-    datainfo_str = 'NORTHSEA'
-
     for s,k in H.C:
         data = H.C[(s,k)]
         # total column
-        TC = plot_totalcolumn(H, data, datainfo_str=datainfo_str, map_region=map_region)
+        TC = plot_totalcolumn(H, data, map_region=region, datainfo_str=label)
         # TC = plot_trajectory(H, T, k, FIGURE=TC)
         filename = '%s_tc_%s.png' % (data.species,
                                      data.timestamp.strftime('%Y%m%dT%H:%M:%S'))
         ofilename = os.path.join(OUTPUT_DIR, filename)
         TC.fig.savefig(ofilename)
 
         # footprint
-        FP = plot_at_level(H, data, datainfo_str=datainfo_str, map_region=map_region)
+        FP = plot_at_level(H, data, map_region=region, datainfo_str=label)
         # FP = plot_trajectory(H, T, k, FIGURE=FP)
         filename = '%s_fp_%s.png' % (data.species,
                                      data.timestamp.strftime('%Y%m%dT%H:%M:%S'))
@@ -65,15 +55,31 @@ def plot_backward(SOURCE_FILE, OUTPUT_DIR):
     """
     Run through plotting routines.
 
-    Expected two params:
+    Expected three params:
     - The netCDF4 file
     - The output dir
+    - The map region
+
+    Optional args:
+    - nested (default False)
+    - label (default empty)
 
     This example is meant to be run with the data in the file
     http://folk.uio.no/johnbur/sharing/stads2_V10.tar
     Use the nested file as, at the time of this writing, the non-nested file is
     broken.
     """
 
-    src_file, out_dir = sys.argv[1:]
-    plot_backward(src_file, out_dir)
+    parser = ArgumentParser(description='Plot the given netCDF4 file')
+    parser.add_argument('src', help='netCDF4 file')
+    parser.add_argument('out', help='output directory')
+    parser.add_argument('region', default='north_sea',
+                        help='region name (e.g. north_sea)')
+    parser.add_argument('--label',
+                        help='label to draw over the plot (e.g. NORTHSEA)')
+    parser.add_argument('--nested', dest='nested', action='store_true',
+                        help='whether the source file is nested or not')
+
+    args = parser.parse_args()
+    plot_backward(args.src, args.out, args.region, nested=args.nested,
+                  label=args.label)

reflexible/__init__.py

@@ -102,13 +102,18 @@ def print_versions():
 
 
 # Import the public functions here
-from .tests.all import test
 from . import conv2netcdf4
-from .scripts.create_ncfile import create_ncfile
-from .data_structures import (Header, Command, Release,
-                              ReleasePoint, Ageclass)
-from .utils import Structure, CacheDict
 from .base_read import read_trajectories, read_partpositions
+from .data_structures import Header, Command, Release, ReleasePoint, Ageclass
 from .flexpart import Flexpart
-from .plotting import (plot_sensitivity, plot_at_level, plot_totalcolumn,
-                       plot_trajectory)
+from .scripts.create_ncfile import create_ncfile
+from .tests.all import test
+from .utils import Structure, CacheDict, curtain_agltoasl, data_range
+
+from .plotting import (
+    plot_sensitivity, plot_at_level, plot_totalcolumn, plot_trajectory,
+    plot_markers, plot_curtain, curtain_for_line)
+
+# XXX Keep this as plot_footprint is referenced in the docs. But maybe remove
+# it in the future, as plot_footprint = plot_at_level
+plot_footprint = plot_at_level

reflexible/conv2netcdf4/flexpart_read.py

@@ -576,56 +576,6 @@ def read_trajectories(H, trajfile='trajectories.txt',
     return RelTraj
 
 
-def groundlevel_for_line(H, X, Y, coords, index=0):
-    """
-    extracts ground level from H.heightnn along a track of lon, lat
-
-    input:  H or H.heightnn (takes the lowest level)
-
-            X, Y ,= H.longitude, H.latitude
-
-            coords = zip(x, y)
-
-    output: groundlevel (a 1-d array with shape (len(flighttrack)
-
-    """
-    try:
-        hgt = H.Heightnn[:, :, 0]
-    except:
-        hgt = H
-
-    # fix for hgt offset
-    hgt = hgt - hgt.min()
-
-    grndlvl = np.zeros(len(coords))
-
-    for i, (x, y) in enumerate(coords):
-
-        I = closest(x, X)
-        J = closest(y, Y)
-        grndlvl[i] = hgt[I, J]
-        grndlvl = np.nan_to_num(grndlvl)
-
-    return grndlvl - grndlvl.min()
-
-
-def curtain_agltoasl(H, curtain_agl, coords, below_gl=0.0):
-    """ converts the agl curtain to asl
-
-        adds .asl_axis attribute to H
-    """
-
-    gl = groundlevel_for_line(H, H.longitude, H.latitude, coords)
-    H.asl_axis = np.linspace(0, H.outheight[-1])
-    xp = H.outheight - H.outheight[0]
-    casl = np.zeros((len(H.asl_axis), len(coords)))
-
-    for i in range(len(coords)):
-        casl[:, i] = np.interp(H.asl_axis, xp + gl[i],
-                               curtain_agl[:, i], left=below_gl)
-    return casl
-
-
 def read_agespectrum(filename, part=False, ndays=20):
     """
     Reads the spectrum.txt files generated from the "make_webpages" scripts.

reflexible/conv2netcdf4/helpers.py

@@ -54,29 +54,6 @@ def get_dir(dir, parent_dir):
     return options_dir, output_dir
 
 
-def data_range(data, min='median'):
-    """
-    return a data range for flexpart data
-
-    optional keyword min = ['median', 'mean', 'min']
-    """
-    dmax = np.nanmax(data)
-    if np.isnan(dmax):
-        dmax = 1e5
-
-    if min == 'mean':
-        dmin = np.mean(data[data.nonzero()])
-    elif min == 'median':
-        dmin = np.median(data[data.nonzero()])
-    else:
-        dmin = np.nanmin(data[data.nonzero()])
-
-    if np.isnan(dmin):
-        dmin = 1e-5
-
-    return [dmin, dmax]
-
-
 def _normdatetime(Y, M, D, h, m, s):
     return datetime.datetime(
         Y, M, D) + datetime.timedelta(hours=h, minutes=m, seconds=s)

reflexible/data_structures.py

@@ -3,6 +3,7 @@
 """
 Definition of the different data structures in reflexible.
 """
+
 import os
 import datetime as dt
 import itertools
@@ -14,6 +15,7 @@
 import xarray as xr
 
 import reflexible
+from reflexible.utils import closest
 
 
 class Header(object):
@@ -346,6 +348,25 @@ def add_trajectory(self):
         """ see :func:`read_trajectories` """
         self.trajectory = reflexible.read_trajectories(self)
 
+    def closest_date(self, dateval, fmt=None):
+        """
+        given a datestring or datetime, tries to find the closest date.
+        if passed a list, assumes it is a list of datetimes
+        """
+
+        if isinstance(dateval, str):
+            if not fmt:
+                if len(dateval) == 8:
+                    fmt = '%Y%m%d'
+                if len(dateval) == 14:
+                    fmt = '%Y%m%d%H%M%S'
+                else:
+                    raise IOError("no format provided for datestring")
+            print("Assuming date format: {0}".format(fmt))
+            dateval = dt.datetime.strptime(dateval, fmt)
+
+        return closest(dateval, self['available_dates_dt'])
+
     @property
     def options(self):
         return {'readp': self.readp,
@@ -509,6 +530,7 @@ def __getitem__(self, item):
 
         return c
 
+
 # TODO: Following John, the get_slabs function should be deprecated
 def get_slabs(Heightnn, data_cube):
     """Preps data_cube for plotting.
@@ -1049,7 +1071,7 @@ def to_file(self, rfile):
             outf.close()
 
 
-class Trajectory(dict, object):
+class Trajectory(dict):
     """ a dictionary-like container for the trajectories.
     This should become a pandas dataframe in the future. """