source: issm/oecreview/Archive/24684-25833/ISSM-25064-25065.diff

Index: ../trunk-jpl/src/m/shp/shpwrite.py
===================================================================
--- ../trunk-jpl/src/m/shp/shpwrite.py  (nonexistent)
+++ ../trunk-jpl/src/m/shp/shpwrite.py  (revision 25065)
@@ -0,0 +1,60 @@
+try:
+    import shapefile
+except ImportError:
+    print("could not import shapefile, PyShp has not been installed, no shapefile reading capabilities enabled")
+
+def shpwrite(shp, filename): #{{{
+    '''
+    SHPREAD - write a shape file from a contour structure
+
+    The current implementation of shpwrite depends on PyShp.
+
+    Usage:
+        shpwrite(shp, filename)
+
+    Example:
+        shpwrite(shp, 'domainoutline.shp')
+
+    See also SHPREAD
+
+    Sources:
+    - https://github.com/GeospatialPython/pyshp
+
+    TODO:
+    - Should we check if there is only one element (see how MATLAB's shaperead 
+    and shapewrite handle single shapes versus multiple shapes)?
+    '''
+
+    # Sanity checks
+    for shape in shp:
+        print(shape)
+
+    if shp[0].Geometry == 'Point':
+        shapeType = 1
+    elif shp[0].Geometry == 'Line':
+        shapeType = 3
+    elif shp[0].Geometry == 'Polygon':
+        shapeType = 5
+    else:
+        raise Exception('shpwrite error: geometry \'{}\' is not currently supported'.format(shp[0].Geometry))
+
+    sf = shapefile.Writer(filename, shapeType=shapeType)
+
+    for i in range(len(shp)):
+        sf.field('name', 'C') # TODO: Allow shape ids/names to be passed through
+        if shapeType == 1: # POINT
+            sf.point(shp[i].x, shp[i].y)
+        elif shapeType == 3: # POLYLINE
+            points = []
+            for j in range(len(shp[i].x)):
+                points.append([shp[i].x[j], shp[i].y[j]])
+            sf.line(line)
+        elif shapeType == 5: # POLYGON
+            points = []
+            for j in range(len(shp[i].x)):
+                points.append([shp[i].x[j], shp[i].y[j]])
+            sf.poly(points)
+        sf.null()
+        sf.record(str(i))
+    sf.close()
+#}}}
Index: ../trunk-jpl/src/m/shp/shpwrite.m
===================================================================
--- ../trunk-jpl/src/m/shp/shpwrite.m   (revision 25064)
+++ ../trunk-jpl/src/m/shp/shpwrite.m   (revision 25065)
@@ -9,10 +9,6 @@
 %
 %   See also SHPREAD
 
-
-%initialize number of profile
-count=0;
-
 contours=struct([]);
 for i=1:length(shp),
        if strcmpi(shp(i).Geometry,'Point'),
@@ -21,11 +17,13 @@
                contours(i).Geometry='Polygon';
        elseif strcmpi(shp(i).Geometry,'Line'),
                contours(i).Geometry='Line';
+       else,
+               error(['shpwrite error: geometry ' shp(i).Geometry ' not currently supported']);
        end
        contours(i).id=i;
        contours(i).X=shp(i).x;
        contours(i).Y=shp(i).y;
 end
-       
+
 shapewrite(contours,filename);
 end
Index: ../trunk-jpl/src/m/mesh/meshintersect3d.py
===================================================================
--- ../trunk-jpl/src/m/mesh/meshintersect3d.py  (revision 25064)
+++ ../trunk-jpl/src/m/mesh/meshintersect3d.py  (revision 25065)
@@ -6,9 +6,11 @@
 
 from pairoptions import pairoptions
 
+
 def meshintersect3d(x, y, z, xs, ys, zs, *args): #{{{
     '''
-    MESHINTERSECT - returns indices (into x, y, and z) of common values between (x, y, z) and (xs, ys, zs) (i.e. x(index) = xs; y(index) = ys).
+    MESHINTERSECT - returns indices (into x, y, and z) of common values between 
+    (x, y, z) and (xs, ys, zs) (i.e. x(index) = xs; y(index) = ys).
     '''
 
     # Process options
Index: ../trunk-jpl/src/m/coordsystems/laea.py
===================================================================
--- ../trunk-jpl/src/m/coordsystems/laea.py     (nonexistent)
+++ ../trunk-jpl/src/m/coordsystems/laea.py     (revision 25065)
@@ -0,0 +1,15 @@
+def laea(lat, long): #{{{
+    '''
+    LAEA - Lambert Azimuthal Equal Area projection at lat, long projection 
+    center.
+
+        Usage:
+            string = laea(45, -90)
+
+        Example:
+            string = laea(45, -90)
+            return string = '+proj=laea +lat_0=45 +lon_0=-90 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs'
+    '''
+
+    return '+proj=laea +lat_0={} +lon_0={} +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs'.format(lat, long)
+#}}}
Index: ../trunk-jpl/src/m/coordsystems/gdaltransform.m
===================================================================
--- ../trunk-jpl/src/m/coordsystems/gdaltransform.m     (revision 25064)
+++ ../trunk-jpl/src/m/coordsystems/gdaltransform.m     (revision 25065)
@@ -19,9 +19,9 @@
 %              Bamber's Greeland projection
 %                      +proj=stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
 %
-%              To get proj.4 string from EPSG, use gdalsrsinfo. Example:
+%      To get proj.4 string from EPSG, use gdalsrsinfo. Example:
 %
-%                      gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //"
+%              gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //"
 
        %give ourselves unique file names
        filename_in  = tempname();
Index: ../trunk-jpl/src/m/geometry/AboveGround.py
===================================================================
--- ../trunk-jpl/src/m/geometry/AboveGround.py  (nonexistent)
+++ ../trunk-jpl/src/m/geometry/AboveGround.py  (revision 25065)
@@ -0,0 +1,8 @@
+import numpy as np
+
+def function AboveGround(lat, long, r, height): #{{{
+    r = r + height
+    x = r * np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(long))
+    y = r * np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(long))
+    z = r * np.sin(np.deg2rad(lat))
+#}}}
Index: ../trunk-jpl/src/m/classes/plotoptions.py
===================================================================
--- ../trunk-jpl/src/m/classes/plotoptions.py   (revision 25064)
+++ ../trunk-jpl/src/m/classes/plotoptions.py   (revision 25065)
@@ -1,4 +1,5 @@
 from collections import OrderedDict
+
 import pairoptions
 
 
@@ -6,8 +7,8 @@
     '''
     PLOTOPTIONS class definition
 
-        Usage:
-            plotoptions = plotoptions(*arg)
+    Usage:
+        plotoptions = plotoptions(*arg)
     '''
 
     def __init__(self, *arg):  # {{{
@@ -35,7 +36,7 @@
     def buildlist(self, *arg):  #{{{
         #check length of input
         if len(arg) % 2:
-            raise TypeError('Invalid parameter / value pair arguments')
+            raise TypeError('Invalid parameter/value pair arguments')
 
         #go through args and build list (like pairoptions)
         rawoptions = pairoptions.pairoptions(*arg)
@@ -102,7 +103,7 @@
                         if len(nums) != 2:
                             continue
                         if False in [x.isdigit() for x in nums]:
-                            raise ValueError('error: in option i - j both i and j must be integers')
+                            raise ValueError('error: in option i-j both i and j must be integers')
                         for j in range(int(nums[0]) - 1, int(nums[1])):
                             self.list[j].addfield(field, rawlist[i][1])
 
Index: ../trunk-jpl/src/m/classes/basin.m
===================================================================
--- ../trunk-jpl/src/m/classes/basin.m  (revision 25064)
+++ ../trunk-jpl/src/m/classes/basin.m  (revision 25065)
@@ -18,42 +18,39 @@
        end
        methods
                function self = basin(varargin) % {{{
-                       switch nargin
-                               case 0
-                                       self=setdefaultparameters(self);
-                               otherwise
+                       self=setdefaultparameters(self);
 
-                                       self=setdefaultparameters(self);
-                                       options=pairoptions(varargin{:}); 
-                       
-                                       %recover field values: 
-                                       self.boundaries=getfieldvalue(options,'boundaries',{});
-                                       self.name=getfieldvalue(options,'name','');
-                                       self.continent=getfieldvalue(options,'continent','');
+                       if nargin>0,
+                               options=pairoptions(varargin{:}); 
+               
+                               %recover field values: 
+                               self.boundaries=getfieldvalue(options,'boundaries',{});
+                               self.name=getfieldvalue(options,'name','');
+                               self.continent=getfieldvalue(options,'continent','');
 
-                                       %figure out projection string: 
+                               %figure out projection string: 
+                               if exist(options,'epsg'),
+                                       if exist(options,'proj'),
+                                               error(['Error in constructor for basin ' self.name '. Cannot supply epsg and proj at the same time!']);
+                                       end
+                                       epsg=getfieldvalue(options,'epsg');
+                                       proj=epsg2proj(epsg); %convert to PROJ.4 format
+                               elseif exist(options,'proj'),
                                        if exist(options,'epsg'),
-                                               if exist(options,'proj'),
-                                                       error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
-                                               end
-                                               epsg=getfieldvalue(options,'epsg');
-                                               %convert into proj4: 
-                                               proj=epsg2proj(epsg);
-                                       elseif exist(options,'proj'),
-                                               if exist(options,'epsg'),
-                                                       error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
-                                               end
-                                               proj=getfieldvalue(options,'proj');
-                                       else
-                                               proj= '+proj=longlat +datum=WGS84 +no_defs';
+                                               error(['Error in constructor for basin ' self.name '. Cannot supply proj and epsg at the same time!']);
                                        end
-                                       self.proj=proj;
+                                       proj=getfieldvalue(options,'proj');
+                               else
+                                       proj='+proj=longlat +datum=WGS84 +no_defs';
+                               end
+
+                               self.proj=proj;
                        end
                end % }}}
                function self = setdefaultparameters(self) % {{{
                        self.name='';
                        self.continent='';
-                       self.proj='+proj=longlat +datum=WGS84 +no_defs '; %EPSG 4326
+                       self.proj='+proj=longlat +datum=WGS84 +no_defs'; %EPSG 4326
                        self.boundaries={};
 
                end % }}}
@@ -92,11 +89,11 @@
                        options=pairoptions(varargin{:});
                        extension=getfieldvalue(options,'extension',1);
 
-                       [path,nme,ext]=fileparts(self.name);
+                       [path,name,ext]=fileparts(self.name);
                        if extension,
-                               output=[nme ext];
+                               output=[name ext];
                        else
-                               output=nme;
+                               output=name;
                        end
                end % }}}
                function plot(self,varargin) % {{{
@@ -173,7 +170,7 @@
                        inshapefile=getfieldvalue(options,'shapefile');
                        outputshapefile=getfieldvalue(options,'outputshapefile','');
 
-                       if (exist(options,'epsgshapefile') &  exist(options,'projshapefile')), 
+                       if (exist(options,'epsgshapefile') & exist(options,'projshapefile')), 
                                error('Basin shapefilecrop error message: cannot specify epsgshapefile and projshapefile at the same time');
                        end
                        if exist(options,'epsgshapefile'),
@@ -185,7 +182,7 @@
                        end
 
 
-                       %create list of contours that have critical length > threshold:  (in lat,long)
+                       %create list of contours that have critical length > threshold (in lat,long):
                        contours=shpread(inshapefile);
                        llist=[];
                        for i=1:length(contours),
@@ -211,8 +208,8 @@
                        flags=zeros(length(contours),1);
                        for i=1:length(contours),
                                h=contours(i); 
-                               in=inpolygon(h.x,h.y,ba.x,ba.y); 
-                               if ~isempty(find(in==0)),
+                               isin=inpolygon(h.x,h.y,ba.x,ba.y); 
+                               if ~isempty(find(isin==0)),
                                        flags(i)=1;
                                end
                        end
Index: ../trunk-jpl/src/m/qmu/helpers.py
===================================================================
--- ../trunk-jpl/src/m/qmu/helpers.py   (revision 25064)
+++ ../trunk-jpl/src/m/qmu/helpers.py   (revision 25065)
@@ -1,8 +1,9 @@
-import numpy as np
 from collections import OrderedDict
 from copy import deepcopy
 
+import numpy as np
 
+
 class struct(object):
     '''An empty struct that can be assigned arbitrary attributes'''
     pass
@@ -10,40 +11,41 @@
 
 class Lstruct(list):
     '''
-An empty struct that can be assigned arbitrary attributes
-    but can also be accesed as a list. Eg. x.y = 'hello', x[:] = ['w', 'o', 'r', 'l', 'd']
+    An empty struct that can be assigned arbitrary attributes but can also be 
+    accesed as a list. Eg. x.y = 'hello', x[:] = ['w', 'o', 'r', 'l', 'd']
 
-Note that 'x' returns the array and x.__dict__ will only return
-    attributes other than the array
+    Note that 'x' returns the array and x.__dict__ will only return attributes 
+    other than the array
 
-List-based and struct-based behaviors work normally, however they are referenced
-    as if the other does not exist; len(x) corresponds only to
-    the list component of x, len(x.a) corresponds to x.a, x.__dict__
-    corresponds only to the non-x-list attributes
+    List-based and struct-based behaviors work normally, however they are 
+    referenced as if the other does not exist; len(x) corresponds only to the 
+    list component of x, len(x.a) corresponds to x.a, x.__dict__ corresponds 
+    only to the non-x-list attributes
 
-Example uses:
+    Examples:
 
-x = Lstruct(1, 2, 3, 4) -> [1, 2, 3, 4]
-    x.a = 'hello'
-    len(x) -> 4
-    x.append(5)
-    len(x) -> 5
-    x[2] -> 3
-    x.a -> 'hello'
-    print x -> [1, 2, 3, 4, 5]
-    x.__dict__ -> {'a':'hello'}
-    x.b = [6, 7, 8, 9]
-    x.b[-1] -> 9
-    len(x.b) -> 4
+        x = Lstruct(1, 2, 3, 4) -> [1, 2, 3, 4]
+        x.a = 'hello'
+        len(x) -> 4
+        x.append(5)
+        len(x) -> 5
+        x[2] -> 3
+        x.a -> 'hello'
+        print x -> [1, 2, 3, 4, 5]
+        x.__dict__ -> {'a':'hello'}
+        x.b = [6, 7, 8, 9]
+        x.b[-1] -> 9
+        len(x.b) -> 4
 
-Other valid constructors:
-          x = Lstruct(1, 2, 3, a = 'hello') -> x.a -> 'hello', x -> [1, 2, 3]
-          x = Lstruct(1, 2, 3)(a = 'hello')
-          x = Lstruct([1, 2, 3], x = 'hello')
-          x = Lstruct((1, 2, 3), a = 'hello')
+    Other valid constructors:
+        x = Lstruct(1, 2, 3, a = 'hello') -> x.a -> 'hello', x -> [1, 2, 3]
+        x = Lstruct(1, 2, 3)(a = 'hello')
+        x = Lstruct([1, 2, 3], x = 'hello')
+        x = Lstruct((1, 2, 3), a = 'hello')
 
-Credit: https://github.com/Vectorized/Python-Attribute-List
-'''
+    Sources:
+    -https://github.com/Vectorized/Python-Attribute-List
+    '''
 
     def __new__(self, *args, **kwargs):
         return super(Lstruct, self).__new__(self, args, kwargs)
@@ -62,36 +64,34 @@
 
 class OrderedStruct(object):
     '''
-A form of dictionary-like structure that maintains the
-    ordering in which its fields/attributes and their
-    corresponding values were added.
+    A form of dictionary-like structure that maintains the ordering in which 
+    its fields/attributes and their corresponding values were added.
 
-OrderedDict is a similar device, however this class
-    can be used as an "ordered struct/class" giving
-    it much more flexibility in practice. It is
+    OrderedDict is a similar device, however this class can be used as an 
+    "ordered struct/class" giving it much more flexibility in practice. It is
     also easier to work with fixed valued keys in-code.
 
-Eg:
-OrderedDict:  # a bit clumsy to use and look at
-    x['y'] = 5
+    Example:
+        OrderedDict:  # a bit clumsy to use and look at
+            x['y'] = 5
 
-OrderedStruct:  # nicer to look at, and works the same way
-    x.y = 5
-    OR
-    x['y'] = 5  # supports OrderedDict-style usage
+        OrderedStruct:  # nicer to look at, and works the same way
+            x.y = 5
+            OR
+            x['y'] = 5  # supports OrderedDict-style usage
 
-Supports: len(x), str(x), for-loop iteration.
-Has methods: x.keys(), x.values(), x.items(), x.iterkeys()
+    Supports: len(x), str(x), for-loop iteration.
+    Has methods: x.keys(), x.values(), x.items(), x.iterkeys()
 
-Usage:
-    x = OrderedStruct()
-    x.y = 5
-    x.z = 6
-    OR
-    x = OrderedStruct('y', 5, 'z', 6)
+    Usage:
+        x = OrderedStruct()
+        x.y = 5
+        x.z = 6
+        OR
+        x = OrderedStruct('y', 5, 'z', 6)
 
-    # note below that the output fields as iterables are always
-    #    in the same order as the inputs
+    note below that the output fields as iterables are always in the same 
+    order as the inputs
 
     x.keys() -> ['y', 'z']
     x.values() -> [5, 6]
@@ -110,16 +110,18 @@
     ('x', 5)
     ('y', 6)
 
-    Note: to access internal fields use dir(x)
-        (input fields will be included, but
-        are not technically internals)
+    Note: to access internal fields use dir(x) (input fields will be included, 
+    but are not technically internals)
     '''
 
     def __init__(self, *args):
-        '''Provided either nothing or a series of strings, construct a structure that will,
-when accessed as a list, return its fields in the same order in which they were provided'''
+        '''
+        Provided either nothing or a series of strings, construct a structure 
+        that will, when accessed as a list, return its fields in the same order 
+        in which they were provided
+        '''
 
-    # keys and values
+        # keys and values
         self._k = []
         self._v = []
 
@@ -184,9 +186,11 @@
             yield(a, b)
 
     def __copy__(self):
-        # shallow copy, hard copies of trivial attributes,
-        # references to structures like lists/OrderedDicts
-        # unless redefined as an entirely different structure
+        '''
+        shallow copy, hard copies of trivial attributes,
+        references to structures like lists/OrderedDicts
+        unless redefined as an entirely different structure
+        '''
         newInstance = type(self)()
         for k, v in list(self.items()):
             exec(('newInstance.%s = v') % (k))
@@ -193,11 +197,13 @@
         return newInstance
 
     def __deepcopy__(self, memo=None):
-        # hard copy of all attributes
-        # same thing but call deepcopy recursively
-        # technically not how it should be done,
-        # (see https://docs.python.org/2/library/copy.html  #copy.deepcopy )
-        # but will generally work in this case
+        '''
+        hard copy of all attributes
+        same thing but call deepcopy recursively
+        technically not how it should be done,
+        (see https://docs.python.org/2/library/copy.html  #copy.deepcopy )
+        but will generally work in this case
+        '''
         newInstance = type(self)()
         for k, v in list(self.items()):
             exec(('newInstance.%s = deepcopy(v)') % (k))
@@ -211,7 +217,7 @@
         i = self._k.index(key)
         k = self._k.pop(i)
         v = self._v.pop(i)
-    #exec('del self.%s')%(key)
+        #exec('del self.%s')%(key)
         return (k, v)
 
     def keys(self):
@@ -226,8 +232,10 @@
 
 def isempty(x):
     '''
-    returns true if object is +  -infinity, NaN, None, '', has length 0, or is an
-    array/matrix composed only of such components (includes mixtures of "empty" types)'''
+    returns true if object is +/-infinity, NaN, None, '', has length 0, or is 
+    an array/matrix composed only of such components (includes mixtures of 
+    "empty" types)
+    '''
 
     if type(x) in [list, np.ndarray, tuple]:
         if np.size(x) == 0:
@@ -261,8 +269,11 @@
 
 
 def fieldnames(x, ignore_internals=True):
-    '''returns a list of fields of x
-    ignore_internals ignores all fieldnames starting with '_' and is True by default'''
+    '''
+    returns a list of fields of x
+    ignore_internals ignores all fieldnames starting with '_' and is True by 
+    default
+    '''
     result = list(vars(x).keys())
 
     if ignore_internals:
@@ -272,8 +283,11 @@
 
 
 def isfield(x, y, ignore_internals=True):
-    '''is y is a field of x?
-    ignore_internals ignores all fieldnames starting with '_' and is True by default'''
+    '''
+    is y is a field of x?
+    ignore_internals ignores all fieldnames starting with '_' and is True by 
+    default
+    '''
     return str(y) in fieldnames(x, ignore_internals)
 
 
@@ -280,7 +294,8 @@
 def fileparts(x):
     '''
     given:   "path/path/.../file_name.ext"
-    returns: [path, file_name, ext] (list of strings)'''
+    returns: [path, file_name, ext] (list of strings)
+    '''
     try:
         a = x[:x.rindex('/')]  #path
         b = x[x.rindex('/') + 1:]  #full filename
@@ -296,17 +311,17 @@
 
 def fullfile(*args):
     '''
-    use:
+    usage:
+        fullfile(path, path, ... , file_name + ext)
 
-    fullfile(path, path, ... , file_name + ext)
     returns: "path/path/.../file_name.ext"
 
     with all arguments as strings with no "/"s
 
     regarding extensions and the '.':
-    as final arguments ('file.doc') or ('file' + '.doc') will work
-    ('final', '.doc'), and the like, will not (you'd get 'final/.doc')
-'''
+        as final arguments ('file.doc') or ('file' + '.doc') will work
+        ('final', '.doc'), and the like, will not (you'd get 'final/.doc')
+    '''
     result = str(args[0])
     for i in range(len(args[1:])):
         # if last argument wasn't empty, add a '/' between it and the next argument
@@ -320,8 +335,11 @@
 def findline(fidi, s):
     '''
     returns full first line containing s (as a string), or None
+
     Note: will include any newlines or tabs that occur in that line,
-    use str(findline(f, s)).strip() to remove these, str() in case result is None'''
+    use str(findline(f, s)).strip() to remove these, str() in case result is 
+    None
+    '''
     for line in fidi:
         if s in line:
             return line
@@ -330,20 +348,21 @@
 
 def empty_nd_list(shape, filler=0., as_numpy_ndarray=False):
     '''
-returns a python list of the size/shape given (shape must be int or tuple)
+    returns a python list of the size/shape given (shape must be int or tuple)
     the list will be filled with the optional second argument
 
     filler is 0.0 by default
 
-    as_numpy_ndarray will return the result as a numpy.ndarray and is False by default
+    as_numpy_ndarray will return the result as a numpy.ndarray and is False by 
+    default
 
-    Note: the filler must be either None/np.nan/float('NaN'), float/double, or int
-        other numpy and float values such as +/- np.inf will also work
+    Note: the filler must be either None/np.nan/float('NaN'), float/double, or 
+    int. other numpy and float values such as +/- np.inf will also work
 
-use:
-    empty_nd_list((5, 5), 0.0)  # returns a 5x5 matrix of 0.0's
-    empty_nd_list(5, None)  # returns a 5 long array of NaN
-'''
+    Usage:
+        empty_nd_list((5, 5), 0.0)  # returns a 5x5 matrix of 0.0's
+        empty_nd_list(5, None)  # returns a 5 long array of NaN
+    '''
     result = np.empty(shape)
     result.fill(filler)
     if not as_numpy_ndarray:
Index: ../trunk-jpl/src/m/mech/analyticaldamage.py
===================================================================
--- ../trunk-jpl/src/m/mech/analyticaldamage.py (revision 25064)
+++ ../trunk-jpl/src/m/mech/analyticaldamage.py (revision 25065)
@@ -1,6 +1,6 @@
 import numpy as np
+
 from averaging import averaging
-#from plotmodel import plotmodel
 from thomasparams import thomasparams
 
 
Index: ../trunk-jpl/src/m/plot/applyoptions.m
===================================================================
--- ../trunk-jpl/src/m/plot/applyoptions.m      (revision 25064)
+++ ../trunk-jpl/src/m/plot/applyoptions.m      (revision 25065)
@@ -307,7 +307,6 @@
        end
 end
 
-
 %text (default value is empty, not NaN...)
 if exist(options,'text');
        textstring=getfieldvalue(options,'text');
@@ -351,7 +350,7 @@
        scaleruler(options);
 end
 
-%streamliness
+%streamlines
 if exist(options,'streamlines'),
        plot_streamlines(md,options);
 end
Index: ../trunk-jpl/src/m/plot/plot_manager.py
===================================================================
--- ../trunk-jpl/src/m/plot/plot_manager.py     (revision 25064)
+++ ../trunk-jpl/src/m/plot/plot_manager.py     (revision 25065)
@@ -1,14 +1,3 @@
-
-from checkplotoptions import checkplotoptions
-from plot_mesh import plot_mesh
-from plot_BC import plot_BC
-from plot_elementnumbering import plot_elementnumbering
-from plot_vertexnumbering import plot_vertexnumbering
-from processmesh import processmesh
-from processdata import processdata
-from plot_unit import plot_unit
-from applyoptions import applyoptions
-
 try:
     from osgeo import gdal
     overlaysupport = True
@@ -16,8 +5,17 @@
     print('osgeo/gdal for python not installed, overlay plots are not enabled')
     overlaysupport = False
 
+from applyoptions import applyoptions
+from checkplotoptions import checkplotoptions
+from plot_BC import plot_BC
+from plot_mesh import plot_mesh
+from plot_elementnumbering import plot_elementnumbering
 if overlaysupport:
     from plot_overlay import plot_overlay
+from plot_unit import plot_unit
+from plot_vertexnumbering import plot_vertexnumbering
+from processdata import processdata
+from processmesh import processmesh
 
 
 def plot_manager(md, options, fig, axgrid, gridindex):
@@ -26,9 +24,12 @@
 
     'fig' is a handle to the figure instance created by plotmodel.
 
-    'ax' is a handle to the axes instance created by plotmodel.  This is
+    'axgrid' is a handle to the axes instance created by plotmodel. This is
     currently generated using matplotlib's AxesGrid toolkit.
 
+    'gridindex' is passed through to specialized plot* functions and 
+    applyoptions.
+
     Usage:
         plot_manager(md, options, fig, ax)
 
@@ -36,9 +37,14 @@
     '''
     #parse options and get a structure of options
     options = checkplotoptions(md, options)
+
     #get data to be plotted
     data = options.getfieldvalue('data')
+
     #add ticklabel has a default option
+    #
+    # TODO: Check why we are doing this and if it is absolutely necessary (see
+    #       also src/m/plot/plot_mesh.py, src/m/plot/applyoptions.py)
     options.addfielddefault('ticklabels', 'on')
 
     ax = axgrid[gridindex]
@@ -57,8 +63,7 @@
     if isinstance(data, str):
         if data == 'mesh':
             plot_mesh(md, options, fig, axgrid, gridindex)
-
-    #fig.delaxes(fig.axes[1])  # hack to remove colorbar after the fact
+            #fig.delaxes(fig.axes[1])  # hack to remove colorbar after the fact
             return
         elif data == 'BC':
             plot_BC(md, options, fig, axgrid, gridindex)
Index: ../trunk-jpl/src/m/plot/plot_overlay.py
===================================================================
--- ../trunk-jpl/src/m/plot/plot_overlay.py     (revision 25064)
+++ ../trunk-jpl/src/m/plot/plot_overlay.py     (revision 25065)
@@ -1,24 +1,31 @@
-import numpy as np
-from processmesh import processmesh
-from processdata import processdata
-from xy2ll import xy2ll
+import os
+
 import matplotlib.pyplot as plt
 import matplotlib as mpl
-import os
 try:
     from mpl_toolkits.basemap import Basemap
 except ImportError:
     print('Basemap toolkit not installed')
+import numpy as np
 try:
     from osgeo import gdal
 except ImportError:
     print('osgeo/gdal for python not installed, plot_overlay is disabled')
 
+from processdata import processdata
+from processmesh import processmesh
+from xy2ll import xy2ll
 
+
 def plot_overlay(md, data, options, ax):
     '''
-    Function for plotting a georeferenced image.  This function is called
-    from within the plotmodel code.
+    Function for plotting a georeferenced image. Called from plot_manager by 
+    call to plotmodel.
+
+    Usage:
+        plot_overlay(md, data, options, ax)
+
+    See also: PLOTMODEL
     '''
 
     x, y, z, elements, is2d, isplanet = processmesh(md, [], options)
@@ -80,7 +87,7 @@
         plt.hist(arr.flatten(), bins=256, range=(0., 1.))
         plt.title('histogram of overlay image, use for setting overlaylims')
         plt.show()
-        plt.sca(ax)  # return to original axes / figure
+        plt.sca(ax) # return to original axes/figure
 
     # get parameters from cropped geotiff
     trans = gtif.GetGeoTransform()
Index: ../trunk-jpl/src/m/plot/plotmodel.m
===================================================================
--- ../trunk-jpl/src/m/plot/plotmodel.m (revision 25064)
+++ ../trunk-jpl/src/m/plot/plotmodel.m (revision 25065)
@@ -9,7 +9,7 @@
 numberofplots=options.numberofplots;
 
 %get number of subplots
-subplotwidth=ceil(sqrt(options.numberofplots));
+subplotwidth=ceil(sqrt(numberofplots));
 
 %if nlines and ncols specified, then bypass.
 if ~exist(options.list{1},'nlines') & ~exist(options.list{1},'ncols')
@@ -103,7 +103,7 @@
                end
        end % }}}
 
-       %Use zbuffer renderer (snoother colors) and white background
+       %Use zbuffer renderer (smoother colors) and white background
        set(f,'Renderer','zbuffer','color',getfieldvalue(options.list{1},'figurebackgroundcolor','w'));
 
        %Go through all data plottable and close window if an error occurs
@@ -124,5 +124,5 @@
                rethrow(me);
        end
 else
-       error('plotmodel error message: no output data found. ');
+       error('plotmodel error message: no output data found.');
 end
Index: ../trunk-jpl/src/m/plot/plot_manager.m
===================================================================
--- ../trunk-jpl/src/m/plot/plot_manager.m      (revision 25064)
+++ ../trunk-jpl/src/m/plot/plot_manager.m      (revision 25065)
@@ -14,9 +14,7 @@
 
 %figure out if this is a special plot
 if ischar(data),
-
        switch data,
-
                case 'boundaries',
                        plot_boundaries(md,options,subplotwidth,i);
                        return;
@@ -32,23 +30,18 @@
                case 'highlightelements',
                        plot_highlightelements(md,options,subplotwidth,i);
                        return;
-
                case 'qmumean',
                        plot_qmumean(md,options,nlines,ncols,i);
                        return;
-
                case 'qmustddev',
                        plot_qmustddev(md,options,nlines,ncols,i);
                        return;
-
                case 'qmuhistnorm',
                        plot_qmuhistnorm(md,options,nlines,ncols,i);
                        return;
-
                case 'qmu_mass_flux_segments',
                        plot_qmu_mass_flux_segments(md,options,nlines,ncols,i);
                        return;
-
                case 'part_hist',
                        plot_parthist(md,options,nlines,ncols,i);
                        return;
@@ -120,10 +113,8 @@
                case 'segments'
                        plot_segments(md,options,subplotwidth,i,data)
                        return
-
                case 'quiver'
                        data=[md.initialization.vx md.initialization.vy]; %Go ahead and try plot_unit
-
                case {'strainrate_tensor','strainrate','strainrate_principal','strainrate_principalaxis1','strainrate_principalaxis2','strainrate_principalaxis3',...
                                'stress_tensor','stress','stress_principal','stress_principalaxis1','stress_principalaxis2','stress_principalaxis3',...
                                'deviatoricstress_tensor','deviatoricstress','deviatoricstress_principal','deviatoricstress_principalaxis1','deviatoricstress_principalaxis2','deviatoricstress_principalaxis3'},
@@ -137,13 +128,10 @@
                        return;
                case 'transient_results',
                        plot_transient_results(md,options,subplotwidth,i);
-
                case 'transient_field',
                        plot_transient_field(md,options,subplotwidth,i);
                        return;
-
        otherwise,
-
                if ismember(data,properties('model')),
                        data=eval(['md.' data ';']);
                else
@@ -158,19 +146,19 @@
        return;
 end
 
-%Figure out if this is a semi-transparent plot.
+%Figure out if this is a Google Maps plot.
 if exist(options,'googlemaps'),
        plot_googlemaps(md,data,options,nlines,ncols,i);
        return;
 end
 
-%Figure out if this is a semi-transparent plot.
+%Figure out if this is a landsat plot.
 if getfieldvalue(options,'landsat',0),
        plot_landsat(md,data,options,nlines,ncols,i);
        return;
 end
 
-%Figure out if this is a semi-transparent plot.
+%Figure out if this is a gridded plot.
 if exist(options,'gridded'),
        plot_gridded(md,data,options,nlines,ncols,i);
        return;
Index: ../trunk-jpl/src/m/shp/shpread.py
===================================================================
--- ../trunk-jpl/src/m/shp/shpread.py   (nonexistent)
+++ ../trunk-jpl/src/m/shp/shpread.py   (revision 25065)
@@ -0,0 +1,136 @@
+from os import path
+try:
+    import shapefile
+except ImportError:
+    print("could not import shapefile, PyShp has not been installed, no shapefile reading capabilities enabled")
+
+from pairoptions import pairoptions
+
+
+def shpread(filename, *args): #{{{
+    '''
+    SHPREAD - read a shape file and build a dict
+
+    This routine reads a shape file .shp and builds a dict array containing the 
+    fields x and y corresponding to the coordinates, one for the filename 
+    of the shp file, for the density, for the nodes, and a field closed to 
+    indicate if the domain is closed. If this initial shapefile is point only, 
+    the fields closed and points are ommited.
+    The first argument is the .shp file to be read and the second one
+    (optional) indicates if the last point shall be read (1 to read it, 0 not 
+    to).
+
+    The current implementation of shpread depends on PyShp.
+
+    Usage:
+        dict = shpread(filename)
+
+    Example:
+        From underling PyShp implementation, "The shapefile format is actually 
+        a collection of three files. You specify the base filename of the 
+        shapefile or the complete filename of any of the shapefile component 
+        files.""
+
+        dict = shpread('domainoutline.shp')
+        OR
+        dict = shpread('domainoutline.dbf')
+        OR
+        dict = shpread('domainoutline')
+
+        "OR any of the other 5+ formats which are potentially part of a 
+        shapefile. The library does not care about file extensions". We do, 
+        however, check that a file with the base filename or base filename with 
+        .shp extension exists.
+
+    See also EXPDOC, EXPWRITEASVERTICES
+
+    Sources:
+    - https://github.com/GeospatialPython/pyshp
+    '''
+
+    #recover options
+    options = pairoptions(*args)
+
+    #some checks
+    if not (path.isfile(filename) or path.isfile(filename + '.shp')):
+        raise RuntimeError('shpread error message: file {} or {}.shp not found!'.format(filename, filename))
+
+    #read shapefile
+    sf = shapefile.Reader(filename)
+
+    Dicts = []
+
+    #retrieve ID (if it exists)
+    fields = sf.fields
+    for i in range(1, len(fields)): # skip over first field, which is "DeletionFlag"
+        if fields[i][0] == 'id':
+            name = str(sf.record(i - 1)[0]) # record index is offset by one, again, because of "DeletionFlag"
+            break
+
+    shapes = sf.shapes()
+    for i in range(len(shapes)):
+        Dict = {}
+        shape = shapes[i]
+        if shape.shapeType == shapefile.POINT:
+            Dict['x'] = shape.points[0][0]
+            Dict['y'] = shape.points[0][1]
+            Dict['density'] = 1
+            Dict['Geometry'] = 'Point'
+        elif shape.shapeType == shapefile.POLYLINE:
+            num_points = len(shape.points)
+            x = []
+            y = []
+            for j in range(num_points):
+                point = shape.points[j]
+                x.append(point[0])
+                y.append(point[1])
+            Dict['x'] = x
+            Dict['y'] = y
+            Dict['nods'] = num_points
+            Dict['density'] = 1
+            Dict['closed'] = 1
+            Dict['BoundingBox'] = shape.bbox
+            Dict['Geometry'] = 'Line'
+        elif shape.shapeType == shapefile.POLYGON:
+            num_points = len(shape.points)
+            x = []
+            y = []
+            for j in range(num_points):
+                point = shape.points[j]
+                x.append(point[0])
+                y.append(point[1])
+            Dict['x'] = x
+            Dict['y'] = y
+            Dict['nods'] = num_points
+            Dict['density'] = 1
+            Dict['closed'] = 1
+            Dict['BoundingBox'] = shape.bbox
+            Dict['Geometry'] =  'Polygon'
+        else:
+            # NOTE: We could do this once before looping over shapes as all 
+            #       shapes in the file must be of the same type, but we would 
+            #       need to have a second check anyway in order to know how to 
+            #       parse the points. So, let's just assume the file is not 
+            #       malformed.
+            #
+            raise Exception('shpread error: geometry {} is not currently supported'.format(shape.shapeTypeName))
+        name = ''
+        fields = sf.fields
+        for j in range(1, len(fields)): # skip over first field, which is "DeletionFlag"
+            fieldname = fields[j][0]
+            # 'id' field gets special treatment
+            if fieldname == 'id':
+                name = str(sf.record(j - 1)[0]) # record index is offset by one, again, because of "DeletionFlag"
+            else:
+                Dict[fieldname] = sf.record(j - 1)[0]
+        Dict['name'] = name
+        Dicts.append(Dict)
+
+    invert = options.getfieldvalue('invert', 0)
+    if invert:
+        for i in range(len(Dicts)):
+            Dicts[i].x = np.flipud(Dicts[i].x)
+            Dicts[i].y = np.flipud(Dicts[i].y)
+
+    return Dicts
+#}}}
Index: ../trunk-jpl/src/m/shp/shpread.m
===================================================================
--- ../trunk-jpl/src/m/shp/shpread.m    (revision 25064)
+++ ../trunk-jpl/src/m/shp/shpread.m    (revision 25065)
@@ -1,22 +1,22 @@
 function Struct=shpread(filename,varargin)
-%SHPREAD - read a shape file and build a Structure
+%SHPREAD - read a shape file and build a struct
 %
-%   This routine reads a shape file .shp and builds a Structure containing the 
-%   fields x and y corresponding to the coordinates, one for the filename of
-%   the shp file, for the density, for the nodes, and a field closed to 
-%   indicate if the domain is closed. 
-%   If this initial shapefile is point only, the fields closed and
-%   points are ommited
-%   The first argument is the .shp file to be read and the second one (optional) 
-%   indicates if the last point shall be read (1 to read it, 0 not to).
+%      This routine reads a shape file .shp and builds a structure array 
+%      containing the fields x and y corresponding to the coordinates, one for the 
+%      filename of the shp file, for the density, for the nodes, and a field 
+%      closed to indicate if the domain is closed. If this initial shapefile is 
+%      point only, the fields closed and points are omitted.
+%      The first argument is the .shp file to be read and the second one 
+%      (optional) indicates if the last point shall be read (1 to read it, 0 not 
+%      to).
 %
-%   Usage:
-%      Struct=shpread(filename)
+%      Usage:
+%              Struct=shpread(filename)
 %
-%   Example:
-%      Struct=shpread('domainoutline.shp')
+%      Example:
+%              Struct=shpread('domainoutline.shp')
 %
-%   See also EXPDOC, EXPWRITEASVERTICES
+%      See also EXPDOC, EXPWRITEASVERTICES
 
 %recover options
 options=pairoptions(varargin{:});
@@ -26,9 +26,6 @@
        error(['shpread error message: file ' filename ' not found!']);
 end
 
-%initialize number of profile
-count=0;
-
 %read shapefile
 shp=shaperead(filename);
 
@@ -35,7 +32,7 @@
 Struct=struct([]);
 fields=fieldnames(shp);
 for i=1:length(shp),
-       if strcmpi(shp(i).Geometry,'Polygon'),
+       if strcmpi(shp(i).Geometry,'Point'),
                x=shp(i).X'; y=shp(i).Y';
                ids=find(isnan(x));
                x(ids)=[]; y(ids)=[];
@@ -42,9 +39,7 @@
 
                Struct(end+1).x=x;
                Struct(end).y=y;
-               Struct(end).nods=length(x);
                Struct(end).density=1;
-               Struct(end).closed=1;
                if isfield(shp,'id'),
                        Struct(end).name=num2str(shp(i).id);
                else
@@ -56,9 +51,7 @@
                                Struct(end).(field)=shp(i).(field);
                        end
                end
-       end
-       
-       if strcmpi(shp(i).Geometry,'Line'),
+       elseif strcmpi(shp(i).Geometry,'Line'),
                x=shp(i).X'; y=shp(i).Y';
                ids=find(isnan(x));
                x(ids)=[]; y(ids)=[];
@@ -79,10 +72,7 @@
                                Struct(end).(field)=shp(i).(field);
                        end
                end
-       end
-
-
-       if strcmpi(shp(i).Geometry,'Point'),
+       elseif strcmpi(shp(i).Geometry,'Polygon'),
                x=shp(i).X'; y=shp(i).Y';
                ids=find(isnan(x));
                x(ids)=[]; y(ids)=[];
@@ -89,7 +79,9 @@
 
                Struct(end+1).x=x;
                Struct(end).y=y;
+               Struct(end).nods=length(x);
                Struct(end).density=1;
+               Struct(end).closed=1;
                if isfield(shp,'id'),
                        Struct(end).name=num2str(shp(i).id);
                else
Index: ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py
===================================================================
--- ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py   (revision 25064)
+++ ../trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py   (revision 25065)
@@ -1,4 +1,3 @@
-import os
 import subprocess
 
 import numpy as np
@@ -9,7 +8,7 @@
 from pairoptions import *
 
 
-def gmshplanet(*varargin):
+def gmshplanet(*args):
     '''
     GMSHPLANET - mesh generation for a sphere. Very specific code for gmsh from $ISSM_DIR/src/demos/simple_geo/sphere.geo
 
@@ -33,7 +32,7 @@
         raise RuntimeError('gmshplanet: Gmsh major version %s not supported!' % gmshmajorversion)
 
     #process options
-    options = pairoptions(*varargin)
+    options = pairoptions(*args)
     #options = deleteduplicates(options, 1)
 
     #recover parameters:
Index: ../trunk-jpl/src/m/coordsystems/gdaltransform.py
===================================================================
--- ../trunk-jpl/src/m/coordsystems/gdaltransform.py    (revision 25064)
+++ ../trunk-jpl/src/m/coordsystems/gdaltransform.py    (revision 25065)
@@ -5,31 +5,35 @@
 
 from loadvars import *
 
+
 def gdaltransform(x, y, proj_in, proj_out): #{{{
     '''
     GDALTRANSFORM - switch from one projection system to another
 
-        Usage:
-            [x, y] = gdaltransform(x1, y1, epsg_in, epsg_out)
+    Usage:
+        [x, y] = gdaltransform(x1, y1, epsg_in, epsg_out)
 
-        Example:
-            [x, y] = gdaltranform(md.mesh.long, md.mesh.lat, 'EPSG:4326', 'EPSG:3031')
+    Example:
+        [x, y] = gdaltranform(md.mesh.long, md.mesh.lat, 'EPSG:4326', 'EPSG:3031')
 
-        For reference:
-            EPSG: 4326 (lat, long)
-            EPSG: 3341 (Greenland,  UPS 45W, 70N)
-            EPSG: 3031 (Antarctica, UPS 0E,  71S)
+    For reference:
+        EPSG: 4326 (lat, long)
+        EPSG: 3341 (Greenland,  UPS 45W, 70N)
+        EPSG: 3031 (Antarctica, UPS 0E,  71S)
 
-        ll2xy default projection Antarctica:
-            +proj = stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
-        ll2xy default projection Greenland:
-            +proj = stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
-        Bamber's Greenland projection
-            +proj = stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
+    ll2xy default projection Antarctica:
+        +proj = stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
+    ll2xy default projection Greenland:
+        +proj = stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
+    Bamber's Greenland projection
+        +proj = stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs
 
-        To get proj.4 string form EPSG, use gdalsrsinfo. Example:
+    To get proj.4 string form EPSG, use gdalsrsinfo. Example:
 
-            gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //"
+        gdalsrsinfo epsg:4326 | grep "PROJ.4" | sed "s/PROJ.4 : //"
+
+    TODO:
+    - Remove shlex and follow subprocess pattern implemented in src/m/coordsystems/epsg2proj.py
     '''
 
     # Give ourselves unique file names
@@ -41,7 +45,7 @@
     file_in.write(b'%8g %8g\n' % (x.flatten(1) y.flatten(1)))
     file_in.close()
 
-    args = shlex.split('gdaltransform -s_srs %s -t_srs %s < %s > %s' % (proj_in, proj_out, filename_in, filename_out))
+    args = shlex.split('gdaltransform -s_srs {} -t_srs {} < {} > {}'.format(proj_in, proj_out, filename_in, filename_out))
     subprocess.check_call(args) # Will raise CalledProcessError if return code is not 0
 
     A = loadvars(filename_out)
@@ -52,4 +56,6 @@
 
     os.remove(filename_in)
     os.remove(filename_out)
+
+    return [xout, yout]
 #}}}
Index: ../trunk-jpl/src/m/geometry/polyarea.py
===================================================================
--- ../trunk-jpl/src/m/geometry/polyarea.py     (nonexistent)
+++ ../trunk-jpl/src/m/geometry/polyarea.py     (revision 25065)
@@ -0,0 +1,28 @@
+import math
+
+import numpy as np
+
+
+def polyarea: #{{{
+    '''
+    POLYAREA - returns the area of the 2-D polygon defined by the vertices in 
+    lists x and y
+    
+    Partial implementation of MATLAB's polyarea function. If x and y are 
+    lists of the same length, then polyarea returns the scalar area of the 
+    polygon defined by x and y.
+
+    Usage:
+        a = polyarea(x, y)
+
+    Sources:
+    - https://www.mathworks.com/help/matlab/ref/polyarea.html
+    - https://stackoverflow.com/questions/24467972/calculate-area-of-polygon-given-x-y-coordinates
+    - https://en.wikipedia.org/wiki/Shoelace_formula
+
+    TODO:
+    - Test that output falls within some tolerance of MATLAB's polyarea 
+    function.
+    '''
+    return 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1)))
+#}}}
Index: ../trunk-jpl/src/m/classes/boundary.py
===================================================================
--- ../trunk-jpl/src/m/classes/boundary.py      (nonexistent)
+++ ../trunk-jpl/src/m/classes/boundary.py      (revision 25065)
@@ -0,0 +1,176 @@
+import math
+
+import numpy as np
+
+from epsg2proj import epsg2proj
+from pairoptions import pairoptions
+from shpread import shpread
+
+
+class boundary(object): #{{{
+    '''
+    BOUNDARY class definition
+
+    Usage:
+        boundary = boundary()
+    '''
+
+    def __init__(self, *args): #{{{
+        self.boundaries = []
+        self.name       = ''
+        self.continent  = ''
+        self.proj       = epsg2proj(4326)
+
+        self.setdefaultparameters()
+
+        if len(args):
+            options = pairoptions(*args)
+
+            #recover field values
+            self.shppath = options.getfieldvalue('shppath', '')
+            self.shpfilename = options.getfieldvalue('shpfilename', '')
+            self.orientation = options.getfieldvalue('orientation', 'normal')
+
+            #figure out projection string:
+            if options.exist('epsg'):
+                if options.exist('proj'):
+                    raise RuntimeError('Error in constructor for boundary {}. Cannot supply epsg and proj at the same time!'.format(self.shppath))
+                epsg = options.getfieldvalue('epsg')
+                proj = epsg2proj(epsg) # convert to PROJ.4 format
+            elif options.exist('proj'):
+                if options.exist('epsg'):
+                    raise RuntimeError('Error in constructor for boundary {}. Cannot supply proj and epsg at the same time!'.format(self.shppath))
+                proj = options.getfieldvalue('proj')
+            else:
+                proj = '+proj=longlat +datum=WGS84 +no_defs'
+
+            self.proj = proj
+    #}}}
+
+    def __repr__(self): #{{{
+        s = '   boundary parameters:\n'
+        s += '{}\n'.format(fielddisplay(self, 'shppath', 'path to filename for this boundary'))
+        s += '{}\n'.format(fielddisplay(self, 'shpfilename', 'shape file name'))
+        s += '{}\n'.format(fielddisplay(self, 'orientation', 'orientation (default is \'normal\', can be \'reverse\')'))
+        s += '{}\n'.format(fielddisplay(self, 'proj', 'shape file projection string (follows PROJ.4 standard)'))
+
+        return s
+    #}}}
+
+    def setdefaultparameters(self): #{{{
+        self.shppath = ''
+        self.shpfilename = ''
+        self.orientation = 'normal'
+        self.proj = '+proj=longlat +datum=WGS84 +no_defs' #EPSG 4326
+
+        return self
+    #}}}
+
+    def name(self): #{{{
+        output = self.shpfilename
+
+        return output
+    #}}}
+
+    def edges(self): #{{{
+        #read domain
+        path, name, ext = fileparts(shp.shpfilename)
+        if ext != 'shp':
+            ext = 'shp'
+        output = shpread('{}/{}.{}'.format(self.shppath, name, ext))
+
+        #do we reverse?
+        if self.orientation == 'reverse':
+            output.x = np.flipud(output.x)
+            output.y = np.flipud(output.y)
+    #}}}
+
+    def plot(self, *args): #{{{
+        #recover options
+        options = pairoptions(*args)
+
+        #parse input
+        # TODO: Sort out what options we should set (see 
+        # src/m/classes/boundary.m)
+        
+        if options.exist('epsg'):
+            if options.exist('proj'):
+                raise RuntimeError('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection')
+            proj = epsg2proj(options.getfieldvalue('epsg'))
+        elif options.exist('proj'):
+            proj = options.getfieldvalue('proj')
+        else:
+            proj = epsg2proj(4326)
+
+        #read domain
+        path, name, ext = fileparts(shp.shpfilename)
+        if ext != 'shp':
+            ext = 'shp'
+        domain = shpread('{}/{}.{}'.format(self.shppath, name, ext))
+
+        #convert boundary to another reference frame #{{{
+        for i in range(len(domain)):
+            try:
+                x, y = gdaltransform(domain[i].x, domain[i].y, self.proj, proj)
+            except error as e:
+                print(e)
+                print(self)
+
+        # TODO: Figure out how to recover figure here: do we pass 'fig' and 
+        # 'ax' in args?
+        #for i in range(len(domain)):
+            # x = domain[i].x * unitmultiplier
+            # y = domain[i].y * unitmultiplier
+            # if len(x) == 1:
+        #}}}
+
+        #TODO: Finish translating from MATLAB after test2004.py runs without plot
+    #}}}
+
+    def checkconsistency(self, *args): #{{{
+        #recover options
+        options = pairoptions(*args)
+        tolerance = getfieldvalue(options, 'tolerance', 1e-5)
+
+        #recover edges
+        edges = self.edges()
+
+        if edges.Geometry == 'Point':
+            return
+        else:
+            #check that we don't have two identical vertices
+            x = edges.x
+            y = edges.y
+            distance = math.sqrt((x[:-2] - x[1:-1]) ** 2 + (y[:-2] - y[1:-1]) ** 2)
+            dmax = distance.max()
+            isidentical = np.where(np.asarray(distance) < dmax * tolerance)
+            for elem in isidentical: # distance should be a 1D array, but return from np.where is a tuple of arrays
+                if len(elem) != 0:
+                    raise Exception('boundary {} has two vertices extermely close to one another'.format(shpfilename))
+
+    def plot3d(self, *args): #{{{
+        #recover options
+        options = pairoptions(*args)
+
+        #parse input
+        # TODO: Sort out what options we should set (see 
+        # src/m/classes/boundary.m)
+        
+        if options.exist('epsg'):
+            if options.exist('proj'):
+                raise RuntimeError('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection')
+            proj = epsg2proj(options.getfieldvalue('epsg'))
+        elif options.exist('proj'):
+            proj = options.getfieldvalue('proj')
+        else:
+            proj = epsg2proj(4326)
+
+        #read domain
+        path, name, ext = fileparts(shp.shpfilename)
+        if ext != 'shp':
+            ext = 'shp'
+        domain = shpread('{}/{}.{}'.format(self.shppath, name, ext))
+
+        #TODO: Finish translating from MATLAB after test2004.py runs without plot
+    #}}}
+#}}}
Index: ../trunk-jpl/src/m/classes/sealevelmodel.m
===================================================================
--- ../trunk-jpl/src/m/classes/sealevelmodel.m  (revision 25064)
+++ ../trunk-jpl/src/m/classes/sealevelmodel.m  (revision 25065)
@@ -162,10 +162,10 @@
                        list=unique(list);
                end % }}}
                function addbasin(self,bas) % {{{
-               if ~strcmpi(class(bas),'basin')
-                       error('addbasin method only takes a ''basin'' class object as input');
-               end;
-               self.basins{end+1}=bas;
+                       if ~strcmpi(class(bas),'basin')
+                               error('addbasin method only takes a ''basin'' class object as input');
+                       end;
+                       self.basins{end+1}=bas;
                end % }}}
                function intersections(self,varargin) % {{{
 
Index: ../trunk-jpl/src/m/classes/sealevelmodel.py
===================================================================
--- ../trunk-jpl/src/m/classes/sealevelmodel.py (revision 25064)
+++ ../trunk-jpl/src/m/classes/sealevelmodel.py (revision 25065)
@@ -10,21 +10,22 @@
 from private import private
 from TwoDToThreeD import TwoDToThreeD
 
+
 class sealevelmodel(object):
     '''
     SEALEVELMODEL class definition
 
-        Usage:
-            slm = sealevelmodel(*args)
+    Usage:
+        slm = sealevelmodel(*args)
 
-            where args is a variable list of options
+        where args is a variable list of options
 
-        Example:
-            slm = sealevelmodel(
-                'icecap', md_greenland,
-                'icecap', md_antarctica,
-                'earth', md_earth
-            )
+    Example:
+        slm = sealevelmodel(
+            'icecap', md_greenland,
+            'icecap', md_antarctica,
+            'earth', md_earth
+        )
     '''
 
     def __init__(self, *args): #{{{
@@ -58,11 +59,11 @@
     #}}}
 
     def __repr__(self): # {{{
-        s = '%s\n' % fielddisplay(self, 'icecaps', 'ice caps')
-        s += '%s\n' % fielddisplay(self, 'earth', 'earth')
-        s += '%s\n' % fielddisplay(self, 'settings', 'settings properties')
-        s += '%s\n' % fielddisplay(self, 'cluster', 'cluster parameters (number of cpus...')
-        s += '%s\n' % fielddisplay(self, 'miscellaneous', 'miscellaneous fields')
+        s = '{}\n'.format(fielddisplay(self, 'icecaps', 'ice caps'))
+        s += '{}\n'.format(fielddisplay(self, 'earth', 'earth'))
+        s += '{}\n'.format(fielddisplay(self, 'settings', 'settings properties'))
+        s += '{}\n'.format(fielddisplay(self, 'cluster', 'cluster parameters (number of cpus...'))
+        s += '{}\n'.format(fielddisplay(self, 'miscellaneous', 'miscellaneous fields'))
     #}}}
 
     def setdefaultparameters(self): # {{{
@@ -84,7 +85,7 @@
         # Is the coupler turned on?
         for i in range(len(slm.icecaps)):
             if slm.icecaps[i].transient.iscoupler == 0:
-                warnings.warn('sealevelmodel checkconsistency error: icecap model %s should have the transient coupler option turned on!' % slm.icecaps[i].miscellaneous.name)
+                warnings.warn('sealevelmodel checkconsistency error: icecap model {} should have the transient coupler option turned on!'.format(slm.icecaps[i].miscellaneous.name))
 
         if slm.earth.transient.iscoupler == 0:
             warnings.warn('sealevelmodel checkconsistency error: earth model should have the transient coupler option turned on!')
@@ -92,18 +93,18 @@
         # Check that the transition vectors have the right size
         for i in range(len(slm.icecaps)):
             if slm.icecaps[i].mesh.numberofvertices != len(slm.earth.slr.transitions[i]):
-                raise RuntimeError('sealevelmodel checkconsistency issue with size of transition vector for ice cap: %i name: %s' % (i, slm.icecaps[i].miscellaneous.name))
+                raise RuntimeError('sealevelmodel checkconsistency issue with size of transition vector for ice cap: {} name: {}'.format(i, slm.icecaps[i].miscellaneous.name))
 
         # Check that run frequency is the same everywhere
         for i in range(len(slm.icecaps)):
             if slm.icecaps[i].slr.geodetic_run_frequency != slm.earth.geodetic_run_frequency:
-                raise RuntimeError('sealevelmodel checkconsistency error: icecap model %s should have the same run frequency as earth!' % slm.icecaps[i].miscellaneous.name)
+                raise RuntimeError('sealevelmodel checkconsistency error: icecap model {} should have the same run frequency as earth!'.format(slm.icecaps[i].miscellaneous.name))
 
         # Make sure steric_rate is the same everywhere
         for i in range(len(slm.icecaps)):
             md = slm.icecaps[i]
             if np.nonzero(md.slr.steric_rate - slm.earth.slr.steric_rate[slm.earth.slr.transitions[i]]) != []:
-                raise RuntimeError('steric rate on ice cap %s is not the same as for the earth' % md.miscellaneous.name)
+                raise RuntimeError('steric rate on ice cap {} is not the same as for the earth'.format(md.miscellaneous.name))
     #}}}
 
     def mergeresults(self): # {{{
@@ -137,7 +138,7 @@
     
     def listcaps(self): # {{{
         for i in range(len(self.icecaps)):
-            print('%i: %s' % (i, self.icecaps[i].miscellaneous.name))
+            print('{}: {}'.format(i, self.icecaps[i].miscellaneous.name))
     #}}}
 
     def continents(self): # {{{
@@ -184,7 +185,7 @@
             yei = mdi.mesh.y[mdi.mesh.elements] * onesmatrix / 3
             zei = mdi.mesh.z[mdi.mesh.elements] * onesmatrix / 3
 
-            print('Computing vertex intersections for basin %s' % self.basins[i].name)
+            print('Computing vertex intersections for basin {}'.format(self.basins[i].name))
 
             self.transitions.append(meshintersect3d(self.earth.mesh.x, self.earth.mesh.y, self.earth.mesh.z, mdi.mesh.x, mdi.mesh.y, mdi.mesh.z, 'force', force))
             self.eltransitions.append(meshintersect3d(xe, ye, ze, xei, yei, zei, 'force', force))
Index: ../trunk-jpl/src/m/classes/basin.py
===================================================================
--- ../trunk-jpl/src/m/classes/basin.py (nonexistent)
+++ ../trunk-jpl/src/m/classes/basin.py (revision 25065)
@@ -0,0 +1,214 @@
+import math
+
+from boundary import boundary
+from epsg2proj import epsg2proj
+from helpers import fileparts
+from laea import laea
+from pairoptions import pairoptions
+from polyarea import polyarea
+from shpread import shpread
+
+
+class basin(object): #{{{
+    '''
+    BASIN class definition
+
+    Usage:
+        basin = basin()
+    '''
+
+    def __init__(self, *args): #{{{
+        self.boundaries = []
+        self.name       = ''
+        self.continent  = ''
+        self.proj       = epsg2proj(4326)
+
+        self.setdefaultparameters()
+
+        if len(args):
+            options = pairoptions(*args)
+
+            #recover field values
+            self.boundaries = options.getfieldvalue('boundaries', [])
+            self.name = options.getfieldvalue('name', '')
+            self.continent = options.getfieldvalue('continent', '')
+
+            # figure out projection string
+            if options.exist('epsg'):
+                if options.exist('proj'):
+                    raise RuntimeError('Error in constructor for basin %s. Cannot supply epsg and proj at the same time!'.format(self.name))
+                epsg = options.getfieldvalue('epsg', '')
+                proj = epsg2proj(epsg)#convert to PROJ.4 format
+            elif options.exist('proj'):
+                if options.exist('epsg'):
+                    raise RuntimeError('Error in constructor for basin {}. Cannot supply proj and epsg at the same time!'.format(self.name))
+                proj = options.getfieldvalue('proj', '')
+            else:
+                proj = '+proj=longlat +datum=WGS84 +no_defs'
+
+            self.proj = proj
+    #}}}
+
+    def __repr__(self): # {{{
+        s = '   basin parameters:\n'
+        s += '{}\n'.format(fielddisplay(self, 'continent', 'continent name'))
+        s += '{}\n'.format(fielddisplay(self, 'name', 'basin name'))
+        s += '{}\n'.format(fielddisplay(self, 'proj', 'basin projection string (follows PROJ.4 standard'))
+        s += '{}\n'.format(fielddisplay(self, 'boundaries','list of boundary objects'))
+        for i in range(len(self.boundaries)):
+            s += '             boundary #{}: {}\n'.format((i + 1), self.boundaries[i])
+
+        return s
+    #}}}
+
+    def setdefaultparameters(self): # {{{
+        self.name       = ''
+        self.continent  = ''
+        self.proj       = '+proj=longlat +datum=WGS84 +no_defs' #EPSG 4326
+        self.boundaries = []
+
+        return self
+    #}}}
+
+    def isnameany(self, *args): #{{{
+        for arg in args:
+            if arg == self.name:
+                return 1
+        return 0
+    #}}}
+
+    def iscontinentany(self, *args): #{{{
+        for arg in args:
+            if arg == self.continent:
+                return 1
+        return 0
+    #}}}
+
+    def outputname(self, *args): #{{{
+        #recover options
+        options = pairoptions(*args)
+        extension = options.getfieldvalue('extension', 1)
+
+        path, name, ext = fileparts(*args)
+        if extension:
+            output = '{}{}'.format(name, ext)
+        else:
+            output = name
+
+        return output
+    #}}}
+
+    def plot(self, *args): #{{{
+        #add option
+        for i in range(len(self.boundaries)):
+            self.boundaries[i].plot('proj', self.proj, *args)
+    #}}}
+
+    def plot3d(self, *args): #{{{
+        #add option
+        for i in range(len(self.boundaries)):
+            self.boundaries[i].plot3d(*args)
+    #}}}
+
+    def contour(self, *args): #{{{
+        #recover options
+        options = pairoptions(*args)
+        x = []
+        y = []
+
+        #go through boundaries, recover edges and project them in the basin epsg reference frame
+        for i in range(len(self.boundaries)):
+            boundary = self.boundaries[i]
+            contour = boundary.edges()
+            contour.x, contour.y = gdaltransform(contour.x, contour.y, boundary.proj, self.proj)
+            x = [[x], [contour.x]]
+            y = [[y], [contour.y]]
+
+        #close the contour
+        if x[-1] != x[0] or y[-1] != y[0]:
+            x.append(x[0])
+            y.append(y[0])
+
+        setattr(out, 'x', x)
+        setattr(out, 'y', y)
+        setattr(out, 'nods', len(x))
+
+        return out
+    #}}}
+
+    def checkconsistency(self, *args): #{{{
+        #recover options
+        options = pairoptions(*args)
+
+        #figure out if two boundaries are identical
+        for i in range(len(self.boundaries)):
+            namei = self.boundaries[i].shpfilename
+            for j in range(1, len(self.boundaries)):
+                namej = self.boundaries[j].shpfilename
+                if namei == namej:
+                    raise RuntimeError('Basin {} has two identical boundaries named {}'.format(self.name, namei))
+
+        #go through boundaries and check their consistency
+        for i in range(len(self.boundaries)):
+            boundary == self.boundaries[i]
+            boundary.checkconsistency()
+    #}}}
+
+    def contourplot(self, *args): #{{{
+        contour = self.contour()
+        plot(contour.x, contour.y, 'r*-')
+    #}}}
+
+    def shapefilecrop(self, *args): #{{{
+        #recover options
+        options = pairoptions(*args)
+        threshold = options.getfieldvalue('threshold', .65) #.65 degrees lat, long
+        inshapefile = options.getfieldvalue('shapefile')
+        outshapefile = options.getfieldvalue('outshapefile', '')
+
+        if options.exist('epsgshapefile') and options.exist('projshapefile'):
+            raise RuntimeError('Basin shapefilecrop error messgae: cannot specify epsgshapefile and projshapefile at the same time')
+
+        if options.exist('epsgshapefile'):
+            projshapefile = epsg2proj(options.getfieldvalue('epsgshapefile'))
+        elif options.exist('projshapefile'):
+            projshapefile = options.getfieldvalue('projshapefile')
+        else:
+            raise RuntimeError('Basin shapefilecrop error message: epsgshapefile or projshapefile should be specified')
+
+        #create list of contours that have critical length > threshold (in lat, long)
+        contours = shpread(inshapefile)
+        llist = []
+        for i in range(len(contours)):
+            contour = contours[i]
+            carea = polyarea(contour.x, contour.y)
+            clength = math.sqrt(carea)
+            if clength < threshold:
+                llist = np.concatenate(llist, i)
+        setattr(contours, llist, [])
+
+        #project onto reference frame
+        for i in range(len(contours)):
+            h = contours[i]
+            h.x, h.y = gdaltransform(h.x, h.y, projshapefile, self.proj)
+            contours[i].x = h.x
+            contours[i].y = h.y
+
+        #only keep the contours that are inside the basin (take centroids)
+        ba = self.contour()
+        flags = np.zeros(len(contours))
+        for i in range(len(contours))
+            h = contours[i]
+            isin = inpolygon(h.x, h.y, ba.x, ba.y)
+            if len(np.where(isin == 0, isin, isin)):
+                flags[i] = 1
+        pos = flags.nonzero()
+        contours[pos] = []
+
+        #Two options
+        if outputshapefile == '':
+            output = contours
+        else:
+            shpwrite(contours, outputshapefile)
+    #}}}
+#}}}
Index: ../trunk-jpl/src/m/classes/boundary.m
===================================================================
--- ../trunk-jpl/src/m/classes/boundary.m       (revision 25064)
+++ ../trunk-jpl/src/m/classes/boundary.m       (revision 25065)
@@ -19,35 +19,32 @@
        end
        methods
                function self = boundary(varargin) % {{{
-                       switch nargin
-                               case 0
-                                       self=setdefaultparameters(self);
-                               otherwise
-                                       self=setdefaultparameters(self);
-                                       options=pairoptions(varargin{:}); 
-                       
-                                       %recover field values: 
-                                       self.shppath=getfieldvalue(options,'shppath','');
-                                       self.shpfilename=getfieldvalue(options,'shpfilename','');
-                                       self.orientation=getfieldvalue(options,'orientation','normal');
+                       self=setdefaultparameters(self);
 
-                                       %figure out projection string: 
+                       if nargin > 0,
+                               options=pairoptions(varargin{:});
+
+                               %recover field values: 
+                               self.shppath=getfieldvalue(options,'shppath','');
+                               self.shpfilename=getfieldvalue(options,'shpfilename','');
+                               self.orientation=getfieldvalue(options,'orientation','normal');
+
+                               %figure out projection string: 
+                               if exist(options,'epsg'),
+                                       if exist(options,'proj'),
+                                               error(['Error in constructor for boundary ' self.shppath '. Cannot supply epsg and proj at the same time!']);
+                                       end
+                                       epsg=getfieldvalue(options,'epsg');
+                                       proj=epsg2proj(epsg); % convert to PROJ.4 format
+                               elseif exist(options,'proj'),
                                        if exist(options,'epsg'),
-                                               if exist(options,'proj'),
-                                                       error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
-                                               end
-                                               epsg=getfieldvalue(options,'epsg');
-                                               %convert into proj4: 
-                                               proj=epsg2proj(epsg);
-                                       elseif exist(options,'proj'),
-                                               if exist(options,'epsg'),
-                                                       error(['Error in constructor for  boundary ' self.shppath '.Cannot supply epsg and proj at the same time!']);
-                                               end
-                                               proj=getfieldvalue(options,'proj');
-                                       else
-                                               proj= '+proj=longlat +datum=WGS84 +no_defs';
+                                               error(['Error in constructor for boundary ' self.shppath '. Cannot supply proj and epsg at the same time!']);
                                        end
-                                       self.proj=proj;
+                                       proj=getfieldvalue(options,'proj');
+                               else
+                                       proj= '+proj=longlat +datum=WGS84 +no_defs';
+                               end
+                               self.proj=proj;
                        end
                end % }}}
                function self = setdefaultparameters(self) % {{{
@@ -54,7 +51,7 @@
                self.shppath='';
                self.shpfilename='';
                self.orientation='normal';
-               self.proj='+proj=longlat +datum=WGS84 +no_defs '; %EPSG 4326
+               self.proj='+proj=longlat +datum=WGS84 +no_defs'; %EPSG 4326
                end % }}}
                function disp(self) % {{{
                        disp(sprintf('   boundary parameters:'));
@@ -69,19 +66,19 @@
                        output=self.shpfilename;
                end % }}}
                function output=edges(self) % {{{
-               
-               %read domain:
-               [path,name,ext]=fileparts(self.shpfilename);
-               if ~strcmpi(ext,'shp'),
-                       ext='shp';
-               end
-               output=shpread([self.shppath '/' name '.' ext]);
 
-               %do we reverse? 
-               if strcmpi(self.orientation,'reverse'),
-                       output.x=flipud(output.x);
-                       output.y=flipud(output.y);
-               end
+                       %read domain:
+                       [path,name,ext]=fileparts(self.shpfilename);
+                       if ~strcmpi(ext,'shp'),
+                               ext='shp';
+                       end
+                       output=shpread([self.shppath '/' name '.' ext]);
+
+                       %do we reverse? 
+                       if strcmpi(self.orientation,'reverse'),
+                               output.x=flipud(output.x);
+                               output.y=flipud(output.y);
+                       end
                end % }}}
                function plot(self,varargin) % {{{
                        %recover options
@@ -100,10 +97,9 @@
                        fontsize=getfieldvalue(options,'fontsize',10);
                        label=getfieldvalue(options,'label','on');
 
-                       if (exist(options,'epsg') &  exist(options,'proj')), 
-                               error('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection');
-                       end
-                       if exist(options,'epsg'),
+                       if (exist(options,'epsg'),
+                               if exist(options,'proj')), 
+                                       error('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection');
                                proj=epsg2proj(getfieldvalue(options,'epsg'));
                        elseif exist(options,'proj'),
                                proj=getfieldvalue(options,'proj');
@@ -118,7 +114,7 @@
                        end
                        domain=shpread([self.shppath '/' name '.' ext]);
 
-                       %convert boundary to another reference frame:  {{{
+                       %convert boundary to another reference frame: {{{
                        for i=1:length(domain),
                                try, 
                                        [x,y] = gdaltransform(domain(i).x,domain(i).y,self.proj,proj);
@@ -168,7 +164,7 @@
                                dmax=max(distance);
                                isidentical=find(distance<dmax*tolerance);
                                if ~isempty(isidentical),
-                                       error(['boundary ' shpfilename ' has two  vertices extremely close to one another']);
+                                       error(['boundary ' shpfilename ' has two vertices extremely close to one another']);
                                end
                        end
                end % }}}
@@ -189,6 +185,16 @@
                        fontsize=getfieldvalue(options,'fontsize',10);
                        label=getfieldvalue(options,'label','on');
 
+                       if (exist(options,'epsg'),
+                               if exist(options,'proj')), 
+                                       error('Boundary plot error message: cannot specify epsg and proj at the same time for plot projection');
+                               proj=epsg2proj(getfieldvalue(options,'epsg'));
+                       elseif exist(options,'proj'),
+                               proj=getfieldvalue(options,'proj');
+                       else
+                               proj=epsg2proj(4326);
+                       end
+
                        %read domain:
                        [path,name,ext]=fileparts(self.shpfilename);
                        if ~strcmpi(ext,'shp'),
Index: ../trunk-jpl/src/m/plot/plotmodel.py
===================================================================
--- ../trunk-jpl/src/m/plot/plotmodel.py        (revision 25064)
+++ ../trunk-jpl/src/m/plot/plotmodel.py        (revision 25065)
@@ -1,6 +1,3 @@
-import numpy as np
-from plotoptions import plotoptions
-from plot_manager import plot_manager
 from math import ceil, sqrt
 
 try:
@@ -8,23 +5,37 @@
     from mpl_toolkits.axes_grid1 import ImageGrid
 except ImportError:
     print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
+import numpy as np
 
+from plot_manager import plot_manager
+from plotoptions import plotoptions
 
+
 def plotmodel(md, *args):
-    ''' at command prompt, type 'plotdoc()' for additional documentation
     '''
+    PLOTMODEL - At command prompt, type 'plotdoc()' for additional 
+    documentation.
 
+    TODO:
+        - Fix 'plotdoc()', as it is not currently working.
+    '''
+
     #First process options
     options = plotoptions(*args)
-    #get number of subplots
-    subplotwidth = ceil(sqrt(options.numberofplots))
+
     #Get figure number and number of plots
     figurenumber = options.figurenumber
     numberofplots = options.numberofplots
 
-    #get hold
-    hold = options.list[0].getfieldvalue('hold', False)
+    #get number of subplots
+    subplotwidth = ceil(sqrt(numberofplots))
 
+    # TODO: Check that commenting this out is correct; we do not need a hold 
+    # under matplotlib, right?
+    #
+    # #get hold
+    # hold = options.list[0].getfieldvalue('hold', False)
+
     #if nrows and ncols specified, then bypass
     if options.list[0].exist('nrows'):
         nrows = options.list[0].getfieldvalue('nrows')
@@ -43,14 +54,18 @@
     nrows = int(nrows)
 
     #check that nrows and ncols were given at the same time!
-    if not nr == nc:
-        raise Exception('error: nrows and ncols need to be specified together, or not at all')
+    if nr != nc:
+        raise Exception('plotmodel error message: nrows and ncols need to be specified together, or not at all')
 
-    #Go through plots
+    # Go through plots
+    #
+    # NOTE: The following is where Python + matplolib differs substantially in 
+    #       implementation and inteface from MATLAB.
+    #
+    # Sources:
+    # - https://matplotlib.org/api/_as_gen/mpl_toolkits.axes_grid1.axes_grid.ImageGrid.html
+    #
     if numberofplots:
-        #if plt.fignum_exists(figurenumber):
-        #       plt.cla()
-
         #if figsize specified
         if options.list[0].exist('figsize'):
             figsize = options.list[0].getfieldvalue('figsize')
@@ -62,36 +77,64 @@
         backgroundcolor = options.list[0].getfieldvalue('backgroundcolor', (0.7, 0.7, 0.7))
         fig.set_facecolor(backgroundcolor)
 
-        translator = {'on': 'each',
-                      'off': 'None',
-                      'one': 'single'}
-    # options needed to define plot grid
+        # options needed to define plot grid
         plotnum = options.numberofplots
         if plotnum == 1:
             plotnum = None
-        direction = options.list[0].getfieldvalue('direction', 'row')  # row, column
-        axes_pad = options.list[0].getfieldvalue('axes_pad', 0.25)
-        add_all = options.list[0].getfieldvalue('add_all', True)  # True, False
-        share_all = options.list[0].getfieldvalue('share_all', True)  # True, False
-        label_mode = options.list[0].getfieldvalue('label_mode', 'L')  # 1, L, all
+
+        # NOTE: The inline comments for each of the following parameters are 
+        #       taken from https://matplotlib.org/api/_as_gen/mpl_toolkits.axes_grid1.axes_grid.ImageGrid.html
+        #
+        direction = options.list[0].getfieldvalue('direction', 'row') # {"row", "column"}, default: "row"
+        axes_pad = options.list[0].getfieldvalue('axes_pad', 0.25) # float or (float, float), default : 0.02; Padding or (horizonal padding, vertical padding) between axes, in inches
+        add_all = options.list[0].getfieldvalue('add_all', True) # bool, default: True
+        share_all = options.list[0].getfieldvalue('share_all', True) # bool, default: False
+        label_mode = options.list[0].getfieldvalue('label_mode', 'L') # {"L", "1", "all"}, default: "L"; Determines which axes will get tick labels: "L": All axes on the left column get vertical tick labels; all axes on the bottom row get horizontal tick labels;. "1": Only the bottom left axes is labelled. "all": all axes are labelled.
+
+        # Translate MATLAB colorbar mode to matplotlib
+        #
+        # TODO:
+        # - Add 'edge' option (research if there is a corresponding option in 
+        #   MATLAB)?
+        #
         colorbar = options.list[0].getfieldvalue('colorbar', 'on')  # on, off (single)
-        cbar_mode = translator[colorbar]
-        cbar_location = options.list[0].getfieldvalue('colorbarpos', 'right')  # right, top
-        cbar_size = options.list[0].getfieldvalue('colorbarsize', '5%')
-        cbar_pad = options.list[0].getfieldvalue('colorbarpad', 0.025)  # None or %
 
-        axgrid = ImageGrid(fig, 111,
-                           nrows_ncols=(nrows, ncols),
-                           direction=direction,
-                           axes_pad=axes_pad,
-                           add_all=add_all,
-                           share_all=share_all,
-                           label_mode=label_mode,
-                           cbar_mode=cbar_mode,
-                           cbar_location=cbar_location,
-                           cbar_size=cbar_size,
-                           cbar_pad=cbar_pad)
+        if colorbar == 'on':
+            colorbar = 'each'
+        elif colorbar == 'one':
+            colorbar = 'single'
+        elif colorbar == 'off':
+            colorbar = 'None'
+        else:
+            raise RuntimeError('plotmodel error: colorbar mode \'{}\' is not a valid option'.format(colorbar))
 
+        cbar_mode = colorbar # {"each", "single", "edge", None }, default: None
+        cbar_location = options.list[0].getfieldvalue('colorbarpos', 'right') # {"left", "right", "bottom", "top"}, default: "right"
+        cbar_pad = options.list[0].getfieldvalue('colorbarpad', 0.025) # float, default: None
+        cbar_size = options.list[0].getfieldvalue('colorbarsize', '5%') # size specification (see Size.from_any), default: "5%"
+
+        # NOTE: Second parameter is:
+        #
+        #   rect(float, float, float, float) or int
+        #
+        # The axes position, as a (left, bottom, width, height) tuple or as a 
+        # three-digit subplot position code (e.g., "121").
+        #
+        axgrid = ImageGrid(
+            fig,
+            111,
+            nrows_ncols=(nrows, ncols),
+            direction=direction,
+            axes_pad=axes_pad,
+            add_all=add_all,
+            share_all=share_all,
+            label_mode=label_mode,
+            cbar_mode=cbar_mode,
+            cbar_location=cbar_location,
+            cbar_size=cbar_size,
+            cbar_pad=cbar_pad
+        )
+
         if cbar_mode == 'None':
             for ax in axgrid.cbar_axes:
                 fig._axstack.remove(ax)
Index: ../trunk-jpl/src/m/plot/plot_unit.py
===================================================================
--- ../trunk-jpl/src/m/plot/plot_unit.py        (revision 25064)
+++ ../trunk-jpl/src/m/plot/plot_unit.py        (revision 25065)
@@ -1,6 +1,4 @@
 from cmaptools import getcolormap, truncate_colormap
-from plot_quiver import plot_quiver
-import numpy as np
 try:
     import matplotlib as mpl
     import matplotlib.pyplot as plt
@@ -9,17 +7,20 @@
     from mpl_toolkits.mplot3d.art3d import Poly3DCollection
 except ImportError:
     print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
+import numpy as np
 
+from plot_quiver import plot_quiver
 
+
 def plot_unit(x, y, z, elements, data, is2d, isplanet, datatype, options, fig, axgrid, gridindex):
-    """
+    '''
     PLOT_UNIT - unit plot, display data
 
     Usage:
-    plot_unit(x, y, z, elements, data, is2d, isplanet, datatype, options)
+        plot_unit(x, y, z, elements, data, is2d, isplanet, datatype, options)
 
     See also: PLOTMODEL, PLOT_MANAGER
-    """
+    '''
     #if we are plotting 3d replace the current axis
     if not is2d:
         axgrid[gridindex].axis('off')
@@ -227,21 +228,16 @@
         else:
             raise ValueError('plot_unit error: 3D node plot not supported yet')
         return
-
     # }}}
     # {{{ plotting P1 Patch (TODO)
-
     elif datatype == 4:
         print('plot_unit message: P1 patch plot not implemented yet')
         return
-
     # }}}
     # {{{ plotting P0 Patch (TODO)
-
     elif datatype == 5:
         print('plot_unit message: P0 patch plot not implemented yet')
         return
-
     # }}}
     else:
         raise ValueError('datatype = %d not supported' % datatype)
Index: ../trunk-jpl/src/m/plot/applyoptions.py
===================================================================
--- ../trunk-jpl/src/m/plot/applyoptions.py     (revision 25064)
+++ ../trunk-jpl/src/m/plot/applyoptions.py     (revision 25065)
@@ -1,17 +1,17 @@
-import numpy as np
 from cmaptools import getcolormap
-from plot_contour import plot_contour
-from plot_streamlines import plot_streamlines
-from expdisp import expdisp
-
 try:
-    from matplotlib.ticker import MaxNLocator
     import matplotlib as mpl
     import matplotlib.pyplot as plt
+    from matplotlib.ticker import MaxNLocator
 except ImportError:
     print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
+import numpy as np
 
+from expdisp import expdisp
+from plot_contour import plot_contour
+from plot_streamlines import plot_streamlines
 
+
 def applyoptions(md, data, options, fig, axgrid, gridindex):
     '''
     APPLYOPTIONS - apply options to current plot
@@ -19,10 +19,10 @@
     'plotobj' is the object returned by the specific plot call used to
     render the data.  This object is used for adding a colorbar.
 
-        Usage:
-            applyoptions(md, data, options)
+    Usage:
+        applyoptions(md, data, options)
 
-        See also: PLOTMODEL, PARSE_OPTIONS
+    See also: PLOTMODEL, PARSE_OPTIONS
     '''
 
     # get handle to current figure and axes instance
@@ -33,9 +33,11 @@
     fontsize = options.getfieldvalue('fontsize', 8)
     fontweight = options.getfieldvalue('fontweight', 'normal')
     fontfamily = options.getfieldvalue('fontfamily', 'sans - serif')
-    font = {'fontsize': fontsize,
-            'fontweight': fontweight,
-            'family': fontfamily}
+    font = {
+        'fontsize': fontsize,
+        'fontweight': fontweight,
+        'family': fontfamily
+    }
     # }}}
     # {{{ title
     if options.exist('title'):