Merge 97d5b61 into db8faa7

doc/source/getting_started.rst

@@ -83,7 +83,7 @@ possible way (i.e. passing to it the unique required argument)::
 
 note that the directory containing the FLEXPART output is the
 `create_ncfile` input. The output of this command is rather verbose,
-but if everything goes fine, we will see something like this at the
+but if everything goes fine, you will see something like this at the
 end::
 
   getting grid for:  ['20100531210000']

reflexible/__init__.py

@@ -21,6 +21,7 @@
 ======
 
     JFB: John F Burkhart <jfburkhart@gmail.com>
+    FA: Francesc Alted <faltet@gmail.com>
 
 CONTRIBUTORS
 ============
@@ -89,15 +90,14 @@ def print_versions():
 
 # Some data sources (for testing purposes mainly)
 datasets = {
-    'Fwd1_V10.0': os.path.join(this_dir, "uio_examples/Fwd1_V10.0"),
-    'Fwd1_V10.1': os.path.join(this_dir, "uio_examples/Fwd1_V10.1/output.mpi"),
-    'Bwd1_V9.02': os.path.join(this_dir, "uio_examples/Bwd1_V9.02/outputs"),
-    'Bwd2_V9.2beta': os.path.join(this_dir,
-                                  "uio_examples/Bwd2_V9.2beta/outputs"),
-    'Fwd1_V9.02': os.path.join(this_dir, "uio_examples/Fwd1_V9.02/outputs"),
-    'Fwd2_V9.02': os.path.join(this_dir, "uio_examples/Fwd2_V9.02/outputs"),
-    'HelloWorld_V9.02': os.path.join(this_dir,
-                                     "uio_examples/HelloWorld_V9.02/outputs"),
+    'Fwd1_V10.1': os.path.join(this_dir, "uio_examples/Fwd1_V10.1"),
+    'Bwd1_V9.02': os.path.join(this_dir, "uio_examples/Bwd1_V9.02"),
+    'Bwd2_V9.2beta': os.path.join(this_dir, "uio_examples/Bwd2_V9.2beta"),
+    'Fwd1_V9.02': os.path.join(this_dir, "uio_examples/Fwd1_V9.02"),
+    'Fwd2_V9.02': os.path.join(this_dir, "uio_examples/Fwd2_V9.02"),
+    'HelloWorld_V9.02': os.path.join(this_dir, "uio_examples/HelloWorld_V9.02"),
+    'Only_Outputs_V9.02': os.path.join(
+        this_dir, "uio_examples/HelloWorld_V9.02/outputs"),
 }
 
 
@@ -108,5 +108,7 @@ def print_versions():
 from .data_structures import (Header, Command, Release,
                               ReleasePoint, Ageclass)
 from .utils import Structure, CacheDict
-
 from .base_read import read_trajectories, read_partpositions
+from .flexpart import Flexpart
+from .plotting import (plot_sensitivity, plot_at_level, plot_totalcolumn,
+                       plot_trajectory)

reflexible/base_read.py

@@ -320,3 +320,4 @@ def read_shortpositions(filename, nspec, xlon0, ylat0, dx, dy, dataframe=False,
     """
     #TODO: read_shortposit see partoutput_short.f90
     pass
+

reflexible/conv2netcdf4/__init__.py

@@ -1,4 +1,6 @@
 # Import the public functions here
 from .legacy_structures import Header, Structure, BinaryFile, FDC
-from .flexpart_read import read_header, read_trajectories, read_command
+from .flexpart_read import (
+    read_header, read_trajectories, read_command, read_releases)
 from .grid_read import read_grid, get_slabs, fill_grids
+from .helpers import get_fpdirs

reflexible/conv2netcdf4/flexpart_read.py

@@ -3,29 +3,227 @@
 from __future__ import print_function
 
 import os
-import re
 import datetime
-from math import pi, sqrt, cos
-
+import re
 import numpy as np
+import bcolz
+from math import pi, sqrt, cos
 
 import reflexible.conv2netcdf4
 from .helpers import closest
 
 
+def read_releases(pathname):
+    """
+    Parses release file in `pathname` and return its contents.
+
+    """
+    signature = open(pathname).read(1)
+    if signature == "&":
+        return read_releases_v10(pathname)
+    else:
+        # This is mainly a placeholder, but we should worry only about v10
+        return read_releases_v9(pathname)
+
+
+def read_releases_v10(pathname):
+    """
+    Parses release file in `pathname` and return a ctable with its contents.
+
+    This is only suited for files in Fortran90 namelist format (FP v10).
+
+    Parameters
+    ----------
+    pathname : pathname
+      Release file name (in Fortran90 namelist format).
+
+    Returns
+    -------
+    A ctable object from bcolz package.
+    """
+    # Setup the container for the data
+    dtype = [('IDATE1', np.int32), ('ITIME1', np.int32),
+             ('IDATE2', np.int32), ('ITIME2', np.int32),
+             ('LON1', np.float32), ('LON2', np.float32),
+             ('LAT1', np.float32), ('LAT2', np.float32),
+             ('Z1', np.float32), ('Z2', np.float32),
+             ('ZKIND', np.int8), ('MASS', np.float32),
+             ('PARTS', np.int32), ('COMMENT', 'S32')]
+    cparams = bcolz.cparams(cname="lz4", clevel=6, shuffle=1)
+    ctable = bcolz.zeros(0, dtype=dtype, cparams=cparams)
+    nrecords = ctable['IDATE1'].chunklen
+    releases = np.zeros(nrecords, dtype=dtype)
+
+    # Prepare for reading the input
+    input_str = open(pathname, 'r').read()
+    marker = "&RELEASE\n"
+    len_marker = len(marker)
+    release_re = r'\S+=\s+[\"|\s](\S+)[,|\"|\w]'
+
+    # Loop over all the marker groups
+    i, n = 0, 0
+    while True:
+        i = input_str.find(marker, i)
+        j = input_str.find(marker, i + 1)
+        n += 1
+        group_block = input_str[i + len_marker: j]
+        i = j
+        values = tuple(re.findall(release_re, group_block))
+        try:
+            releases[(n - 1) % nrecords] = values
+        except ValueError:
+            print("Problem at: group: %d, %s" % (n, group_block))
+            print("values:", values)
+            raise
+        if (n % nrecords) == 0:
+            ctable.append(releases)
+        if (i == -1) or (j == -1):
+            break  # marker is not found anymore
+    # Remainder
+    ctable.append(releases[:n % nrecords])
+    ctable.flush()
+
+    return ctable
+
+
+def read_releases_v9(path):
+    """Read metadata from a RELEASES path and return it as a dict.
+
+    Only 'release_point_names' entry returned.
+    """
+    rpnames = []
+    with open(path) as f:
+        prev_line = None
+        for line in f:
+            if prev_line is not None and "comment" in line:
+                rpnames.append(prev_line.strip())
+            prev_line = line
+    # Return just the release point names for now
+    return {"release_point_names": np.array(rpnames, dtype="S45")}
+
+
+def read_releases_old(path, headerrows=11):
+    """
+    Reads a FLEXPART releases file.
+
+    .. note::
+        Assumes releases file has a header of 11 lines. Use
+        option "headerrows" to override this.
+
+    USAGE::
+
+        > A = read_releases_v9(path)
+
+    where filepath is either a file object or a path.
+
+
+    Returns
+      a record array with fields:
+
+      ============      ==========================
+      fields            description
+      ============      ==========================
+      start_time        datetime start
+      end_time          datetime end
+      lllon             lower left longitude
+      llat              lower left latitude
+      urlon             upper right longitude
+      urlat             upper right latitude
+      altunit           1=magl, 2=masl, 3=hPa
+      elv1              lower z level
+      elv2              upper z level
+      numpart           numparticles
+      mass              mass for each spec, so the
+                        array will actually have
+                        fields: mass0, mass1,..
+      id                ID for each release
+      ============      ==========================
+
+    """
+
+    def getfile_lines(infile):
+        """ returns all lines from a file or file string
+        reverts to beginning of file."""
+
+        if isinstance(infile, str):
+            return open(infile, 'r').readlines()
+        else:
+            infile.seek(0)
+            return infile.readlines()
+
+    lines = getfile_lines(path)
+    lines = [i.strip() for i in lines]  # clean line ends
+
+    # we know nspec is at line 11
+    nspec = int(lines[headerrows])
+    blocklength = headerrows + nspec
+    spec = []
+    for i in range(nspec):
+        spec.append(int(lines[headerrows + 1 + i]))
+    indx = headerrows + 1 + (i + 2)
+    blocks = []
+    for i in range(indx, len(lines), blocklength + 1):
+        blocks.append(lines[i:i + blocklength])
+    for b in blocks:
+        b[0] = datetime.datetime.strptime(b[0], '%Y%m%d %H%M%S')
+        b[1] = datetime.datetime.strptime(b[1], '%Y%m%d %H%M%S')
+        b[2:6] = [np.float(i) for i in b[2:6]]
+        b[6] = int(b[6])
+        b[7] = float(b[7])
+        b[8] = float(b[8])
+        b[9] = int(b[9])
+        for i in range(nspec):
+            b[10 + i] = float(b[10 + i])
+        b = tuple(b)
+
+    names = ['start_time', 'end_time', 'lllon', 'lllat', 'urlon', 'urlat',
+             'altunit', 'elv1', 'elv2', 'numpart']
+    # formats = [object, object, np.float, np.float, np.float, np.float,\
+    #                      int, np.float, np.float, int]
+    for i in range(nspec):
+        names.append('mass%s' % i)
+        # formats.append(np.float)
+    names.append('id')
+    # formats.append('S30')
+
+    # dtype = {'names':names, 'formats':formats}
+    # RELEASES = np.rec.array(blocks,dtype=dtype)
+    return np.rec.fromrecords(blocks, names=names)
+
+
 def read_command(path, headerrows=7):
     with open(path, 'r') as comfile:
-        for i in range(headerrows):
+        if "&COMMAND" in comfile.readline():
+            return read_command_v10(path, 1)
+        for i in range(headerrows - 1):
             comfile.readline()
         signature = comfile.readline().strip()
         if signature[:5] == "1. __":
-            command = read_command_new(path, headerrows)
+            return read_command_v9(path, headerrows)
         else:
-            command = read_command_old(path, headerrows)
-    return command
+            return read_command_old(path, headerrows)
 
 
-def read_command_new(path, headerrows):
+def read_command_v10(path, headerrows):
+    lines = open(path, 'r').readlines()
+    command_vals = [i.strip() for i in lines[headerrows:]]  # clean line ends
+    commands = {}
+    for command in command_vals:
+        if '=' not in command:
+            break  # this is the end of commands
+        key, val = command.split("=")
+        val = val[:-1]   # get rid of the trailing ,
+        try:
+            commands[key] = int(val)
+        except ValueError:
+            try:
+                commands[key] = float(val)
+            except ValueError:
+                commands[key] = val
+    return commands
+
+
+def read_command_v9(path, headerrows):
     """Quick and dirty approach for reading COMMAND file for FP V9"""
     COMMAND_KEYS = (
         'SIM_DIR',
@@ -209,95 +407,6 @@ def read_command_old(path, headerrows):
     return COMMAND
 
 
-def read_releases(path, headerrows=11):
-    """
-    Reads a FLEXPART releases file.
-
-    .. note::
-        Assumes releases file has a header of 11 lines. Use
-        option "headerrows" to override this.
-
-    USAGE::
-
-        > A = read_releases(filepath)
-
-    where filepath is either a file object or a path.
-
-
-    Returns
-      a record array with fields:
-
-      ============      ==========================
-      fields            description
-      ============      ==========================
-      start_time        datetime start
-      end_time          datetime end
-      lllon             lower left longitude
-      llat              lower left latitude
-      urlon             upper right longitude
-      urlat             upper right latitude
-      altunit           1=magl, 2=masl, 3=hPa
-      elv1              lower z level
-      elv2              upper z level
-      numpart           numparticles
-      mass              mass for each spec, so the
-                        array will actually have
-                        fields: mass0, mass1,..
-      id                ID for each release
-      ============      ==========================
-
-    """
-
-    def getfile_lines(infile):
-        """ returns all lines from a file or file string
-        reverts to beginning of file."""
-
-        if isinstance(infile, str):
-            return open(infile, 'r').readlines()
-        else:
-            infile.seek(0)
-            return infile.readlines()
-
-    lines = getfile_lines(path)
-    lines = [i.strip() for i in lines]  # clean line ends
-
-    # we know nspec is at line 11
-    nspec = int(lines[headerrows])
-    blocklength = headerrows + nspec
-    spec = []
-    for i in range(nspec):
-        spec.append(int(lines[headerrows + 1 + i]))
-    indx = headerrows + 1 + (i + 2)
-    blocks = []
-    for i in range(indx, len(lines), blocklength + 1):
-        blocks.append(lines[i:i + blocklength])
-    for b in blocks:
-        b[0] = datetime.datetime.strptime(b[0], '%Y%m%d %H%M%S')
-        b[1] = datetime.datetime.strptime(b[1], '%Y%m%d %H%M%S')
-        b[2:6] = [np.float(i) for i in b[2:6]]
-        b[6] = int(b[6])
-        b[7] = float(b[7])
-        b[8] = float(b[8])
-        b[9] = int(b[9])
-        for i in range(nspec):
-            b[10 + i] = float(b[10 + i])
-        b = tuple(b)
-
-    names = ['start_time', 'end_time', 'lllon', 'lllat', 'urlon', 'urlat',
-             'altunit', 'elv1', 'elv2', 'numpart']
-    # formats = [object, object, np.float, np.float, np.float, np.float,\
-    #                      int, np.float, np.float, int]
-    for i in range(nspec):
-        names.append('mass%s' % i)
-        # formats.append(np.float)
-    names.append('id')
-    # formats.append('S30')
-
-    # dtype = {'names':names, 'formats':formats}
-    # RELEASES = np.rec.array(blocks,dtype=dtype)
-    return np.rec.fromrecords(blocks, names=names)
-
-
 def read_trajectories(H, trajfile='trajectories.txt',
                       ncluster=5,
                       ageclasses=20):