Index: /issm/trunk-jpl/src/m/classes/basin.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/basin.m	(revision 25064)
+++ /issm/trunk-jpl/src/m/classes/basin.m	(revision 25065)
@@ -19,34 +19,31 @@
 	methods
 		function self = basin(varargin) % {{{
-			switch nargin
-				case 0
-					self=setdefaultparameters(self);
-				otherwise
-
-					self=setdefaultparameters(self);
-					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: 
+			self=setdefaultparameters(self);
+
+			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: 
+				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 % }}}
@@ -54,5 +51,5 @@
 			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={};
 
@@ -93,9 +90,9 @@
 			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 % }}}
@@ -174,5 +171,5 @@
 			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
@@ -186,5 +183,5 @@
 
 
-			%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=[];
@@ -212,6 +209,6 @@
 			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
Index: /issm/trunk-jpl/src/m/classes/basin.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/basin.py	(revision 25065)
+++ /issm/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: /issm/trunk-jpl/src/m/classes/boundary.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/boundary.m	(revision 25064)
+++ /issm/trunk-jpl/src/m/classes/boundary.m	(revision 25065)
@@ -20,33 +20,30 @@
 	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');
-
-					%figure out projection string: 
+			self=setdefaultparameters(self);
+
+			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';
-					end
-					self.proj=proj;
+						error(['Error in constructor for boundary ' self.shppath '. Cannot supply proj and epsg at the same time!']);
+					end
+					proj=getfieldvalue(options,'proj');
+				else
+					proj= '+proj=longlat +datum=WGS84 +no_defs';
+				end
+				self.proj=proj;
 			end
 		end % }}}
@@ -55,5 +52,5 @@
 		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) % {{{
@@ -70,17 +67,17 @@
 		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) % {{{
@@ -101,8 +98,7 @@
 			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'),
@@ -119,5 +115,5 @@
 			domain=shpread([self.shppath '/' name '.' ext]);
 
-			%convert boundary to another reference frame:  {{{
+			%convert boundary to another reference frame: {{{
 			for i=1:length(domain),
 				try, 
@@ -169,5 +165,5 @@
 				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
@@ -190,4 +186,14 @@
 			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);
Index: /issm/trunk-jpl/src/m/classes/boundary.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/boundary.py	(revision 25065)
+++ /issm/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: /issm/trunk-jpl/src/m/classes/plotoptions.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/plotoptions.py	(revision 25064)
+++ /issm/trunk-jpl/src/m/classes/plotoptions.py	(revision 25065)
@@ -1,3 +1,4 @@
 from collections import OrderedDict
+
 import pairoptions
 
@@ -7,6 +8,6 @@
     PLOTOPTIONS class definition
 
-        Usage:
-            plotoptions = plotoptions(*arg)
+    Usage:
+        plotoptions = plotoptions(*arg)
     '''
 
@@ -36,5 +37,5 @@
         #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)
@@ -103,5 +104,5 @@
                             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: /issm/trunk-jpl/src/m/classes/sealevelmodel.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/sealevelmodel.m	(revision 25064)
+++ /issm/trunk-jpl/src/m/classes/sealevelmodel.m	(revision 25065)
@@ -163,8 +163,8 @@
 		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: /issm/trunk-jpl/src/m/classes/sealevelmodel.py
===================================================================
--- /issm/trunk-jpl/src/m/classes/sealevelmodel.py	(revision 25064)
+++ /issm/trunk-jpl/src/m/classes/sealevelmodel.py	(revision 25065)
@@ -11,19 +11,20 @@
 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
+        )
     '''
 
@@ -59,9 +60,9 @@
 
     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'))
     #}}}
 
@@ -85,5 +86,5 @@
         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:
@@ -93,10 +94,10 @@
         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
@@ -104,5 +105,5 @@
             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))
     #}}}
 
@@ -138,5 +139,5 @@
     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))
     #}}}
 
@@ -185,5 +186,5 @@
             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))
Index: /issm/trunk-jpl/src/m/coordsystems/gdaltransform.m
===================================================================
--- /issm/trunk-jpl/src/m/coordsystems/gdaltransform.m	(revision 25064)
+++ /issm/trunk-jpl/src/m/coordsystems/gdaltransform.m	(revision 25065)
@@ -20,7 +20,7 @@
 %			+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
Index: /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py
===================================================================
--- /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py	(revision 25064)
+++ /issm/trunk-jpl/src/m/coordsystems/gdaltransform.py	(revision 25065)
@@ -6,29 +6,33 @@
 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
     '''
 
@@ -42,5 +46,5 @@
     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
 
@@ -53,3 +57,5 @@
     os.remove(filename_in)
     os.remove(filename_out)
+
+    return [xout, yout]
 #}}}
Index: /issm/trunk-jpl/src/m/coordsystems/laea.py
===================================================================
--- /issm/trunk-jpl/src/m/coordsystems/laea.py	(revision 25065)
+++ /issm/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: /issm/trunk-jpl/src/m/geometry/AboveGround.py
===================================================================
--- /issm/trunk-jpl/src/m/geometry/AboveGround.py	(revision 25065)
+++ /issm/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: /issm/trunk-jpl/src/m/geometry/polyarea.py
===================================================================
--- /issm/trunk-jpl/src/m/geometry/polyarea.py	(revision 25065)
+++ /issm/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: /issm/trunk-jpl/src/m/mech/analyticaldamage.py
===================================================================
--- /issm/trunk-jpl/src/m/mech/analyticaldamage.py	(revision 25064)
+++ /issm/trunk-jpl/src/m/mech/analyticaldamage.py	(revision 25065)
@@ -1,5 +1,5 @@
 import numpy as np
+
 from averaging import averaging
-#from plotmodel import plotmodel
 from thomasparams import thomasparams
 
Index: /issm/trunk-jpl/src/m/mesh/meshintersect3d.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/meshintersect3d.py	(revision 25064)
+++ /issm/trunk-jpl/src/m/mesh/meshintersect3d.py	(revision 25065)
@@ -7,7 +7,9 @@
 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).
     '''
 
Index: /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py
===================================================================
--- /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py	(revision 25064)
+++ /issm/trunk-jpl/src/m/mesh/planet/gmsh/gmshplanet.py	(revision 25065)
@@ -1,3 +1,2 @@
-import os
 import subprocess
 
@@ -10,5 +9,5 @@
 
 
-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
@@ -34,5 +33,5 @@
 
     #process options
-    options = pairoptions(*varargin)
+    options = pairoptions(*args)
     #options = deleteduplicates(options, 1)
 
Index: /issm/trunk-jpl/src/m/plot/applyoptions.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/applyoptions.m	(revision 25064)
+++ /issm/trunk-jpl/src/m/plot/applyoptions.m	(revision 25065)
@@ -307,5 +307,4 @@
 	end
 end
-
 
 %text (default value is empty, not NaN...)
@@ -352,5 +351,5 @@
 end
 
-%streamliness
+%streamlines
 if exist(options,'streamlines'),
 	plot_streamlines(md,options);
Index: /issm/trunk-jpl/src/m/plot/applyoptions.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/applyoptions.py	(revision 25064)
+++ /issm/trunk-jpl/src/m/plot/applyoptions.py	(revision 25065)
@@ -1,14 +1,14 @@
+from cmaptools import getcolormap
+try:
+    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 cmaptools import getcolormap
+
+from expdisp import expdisp
 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
-except ImportError:
-    print("could not import pylab, matplotlib has not been installed, no plotting capabilities enabled")
 
 
@@ -20,8 +20,8 @@
     render the data.  This object is used for adding a colorbar.
 
-        Usage:
-            applyoptions(md, data, options)
-
-        See also: PLOTMODEL, PARSE_OPTIONS
+    Usage:
+        applyoptions(md, data, options)
+
+    See also: PLOTMODEL, PARSE_OPTIONS
     '''
 
@@ -34,7 +34,9 @@
     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
Index: /issm/trunk-jpl/src/m/plot/plot_manager.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_manager.m	(revision 25064)
+++ /issm/trunk-jpl/src/m/plot/plot_manager.m	(revision 25065)
@@ -15,7 +15,5 @@
 %figure out if this is a special plot
 if ischar(data),
-
 	switch data,
-
 		case 'boundaries',
 			plot_boundaries(md,options,subplotwidth,i);
@@ -33,21 +31,16 @@
 			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);
@@ -121,8 +114,6 @@
 			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',...
@@ -138,11 +129,8 @@
 		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 ';']);
@@ -159,5 +147,5 @@
 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);
@@ -165,5 +153,5 @@
 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);
@@ -171,5 +159,5 @@
 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);
Index: /issm/trunk-jpl/src/m/plot/plot_manager.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_manager.py	(revision 25064)
+++ /issm/trunk-jpl/src/m/plot/plot_manager.py	(revision 25065)
@@ -1,13 +1,2 @@
-
-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
@@ -17,6 +6,15 @@
     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
 
 
@@ -27,6 +25,9 @@
     '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:
@@ -37,7 +38,12 @@
     #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')
 
@@ -58,6 +64,5 @@
         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':
Index: /issm/trunk-jpl/src/m/plot/plot_overlay.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_overlay.py	(revision 25064)
+++ /issm/trunk-jpl/src/m/plot/plot_overlay.py	(revision 25065)
@@ -1,13 +1,11 @@
-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
@@ -15,9 +13,18 @@
     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
     '''
 
@@ -81,5 +88,5 @@
         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
Index: /issm/trunk-jpl/src/m/plot/plot_unit.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_unit.py	(revision 25064)
+++ /issm/trunk-jpl/src/m/plot/plot_unit.py	(revision 25065)
@@ -1,5 +1,3 @@
 from cmaptools import getcolormap, truncate_colormap
-from plot_quiver import plot_quiver
-import numpy as np
 try:
     import matplotlib as mpl
@@ -10,15 +8,18 @@
 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:
@@ -228,19 +229,14 @@
             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:
Index: /issm/trunk-jpl/src/m/plot/plotmodel.m
===================================================================
--- /issm/trunk-jpl/src/m/plot/plotmodel.m	(revision 25064)
+++ /issm/trunk-jpl/src/m/plot/plotmodel.m	(revision 25065)
@@ -10,5 +10,5 @@
 
 %get number of subplots
-subplotwidth=ceil(sqrt(options.numberofplots));
+subplotwidth=ceil(sqrt(numberofplots));
 
 %if nlines and ncols specified, then bypass.
@@ -104,5 +104,5 @@
 	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'));
 
@@ -125,4 +125,4 @@
 	end
 else
-	error('plotmodel error message: no output data found. ');
+	error('plotmodel error message: no output data found.');
 end
Index: /issm/trunk-jpl/src/m/plot/plotmodel.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plotmodel.py	(revision 25064)
+++ /issm/trunk-jpl/src/m/plot/plotmodel.py	(revision 25065)
@@ -1,5 +1,2 @@
-import numpy as np
-from plotoptions import plotoptions
-from plot_manager import plot_manager
 from math import ceil, sqrt
 
@@ -9,20 +6,34 @@
 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
@@ -44,12 +55,16 @@
 
     #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'):
@@ -63,33 +78,61 @@
         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':
Index: /issm/trunk-jpl/src/m/qmu/helpers.py
===================================================================
--- /issm/trunk-jpl/src/m/qmu/helpers.py	(revision 25064)
+++ /issm/trunk-jpl/src/m/qmu/helpers.py	(revision 25065)
@@ -1,5 +1,6 @@
-import numpy as np
 from collections import OrderedDict
 from copy import deepcopy
+
+import numpy as np
 
 
@@ -11,38 +12,39 @@
 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']
-
-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
-
-Example uses:
-
-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')
-
-Credit: https://github.com/Vectorized/Python-Attribute-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']
+
+    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
+
+    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
+
+    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')
+
+    Sources:
+    -https://github.com/Vectorized/Python-Attribute-List
+    '''
 
     def __new__(self, *args, **kwargs):
@@ -63,34 +65,32 @@
 class OrderedStruct(object):
     '''
-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
+    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
     also easier to work with fixed valued keys in-code.
 
-Eg:
-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
-
-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)
-
-    # note below that the output fields as iterables are always
-    #    in the same order as the inputs
+    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
+
+    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)
+
+    note below that the output fields as iterables are always in the same 
+    order as the inputs
 
     x.keys() -> ['y', 'z']
@@ -111,14 +111,16 @@
     ('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'''
-
-    # keys and values
+        '''
+        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
         self._k = []
         self._v = []
@@ -185,7 +187,9 @@
 
     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()):
@@ -194,9 +198,11 @@
 
     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()):
@@ -212,5 +218,5 @@
         k = self._k.pop(i)
         v = self._v.pop(i)
-    #exec('del self.%s')%(key)
+        #exec('del self.%s')%(key)
         return (k, v)
 
@@ -227,6 +233,8 @@
 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]:
@@ -262,6 +270,9 @@
 
 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())
 
@@ -273,6 +284,9 @@
 
 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)
 
@@ -281,5 +295,6 @@
     '''
     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
@@ -297,7 +312,7 @@
 def fullfile(*args):
     '''
-    use:
-
-    fullfile(path, path, ... , file_name + ext)
+    usage:
+        fullfile(path, path, ... , file_name + ext)
+
     returns: "path/path/.../file_name.ext"
 
@@ -305,7 +320,7 @@
 
     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:])):
@@ -321,6 +336,9 @@
     '''
     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:
@@ -331,18 +349,19 @@
 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
-
-    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
-'''
+    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
+
+    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)
Index: /issm/trunk-jpl/src/m/shp/shpread.m
===================================================================
--- /issm/trunk-jpl/src/m/shp/shpread.m	(revision 25064)
+++ /issm/trunk-jpl/src/m/shp/shpread.m	(revision 25065)
@@ -1,21 +1,21 @@
 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
@@ -27,7 +27,4 @@
 end
 
-%initialize number of profile
-count=0;
-
 %read shapefile
 shp=shaperead(filename);
@@ -36,5 +33,24 @@
 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)=[];
+
+		Struct(end+1).x=x;
+		Struct(end).y=y;
+		Struct(end).density=1;
+		if isfield(shp,'id'),
+			Struct(end).name=num2str(shp(i).id);
+		else
+			Struct(end).name='';
+		end
+		for j=1:length(fields),
+			field=fields{j};
+			if ~(strcmpi(field,'X') | strcmpi(field,'Y') | strcmpi(field,'id')),
+				Struct(end).(field)=shp(i).(field);
+			end
+		end
+	elseif strcmpi(shp(i).Geometry,'Line'),
 		x=shp(i).X'; y=shp(i).Y';
 		ids=find(isnan(x));
@@ -57,7 +73,5 @@
 			end
 		end
-	end
-	
-	if strcmpi(shp(i).Geometry,'Line'),
+	elseif strcmpi(shp(i).Geometry,'Polygon'),
 		x=shp(i).X'; y=shp(i).Y';
 		ids=find(isnan(x));
@@ -69,26 +83,4 @@
 		Struct(end).density=1;
 		Struct(end).closed=1;
-		if isfield(shp,'id'),
-			Struct(end).name=num2str(shp(i).id);
-		else
-			Struct(end).name='';
-		end
-		for j=1:length(fields),
-			field=fields{j};
-			if ~(strcmpi(field,'X') | strcmpi(field,'Y') | strcmpi(field,'id')),
-				Struct(end).(field)=shp(i).(field);
-			end
-		end
-	end
-
-
-	if strcmpi(shp(i).Geometry,'Point'),
-		x=shp(i).X'; y=shp(i).Y';
-		ids=find(isnan(x));
-		x(ids)=[]; y(ids)=[];
-
-		Struct(end+1).x=x;
-		Struct(end).y=y;
-		Struct(end).density=1;
 		if isfield(shp,'id'),
 			Struct(end).name=num2str(shp(i).id);
Index: /issm/trunk-jpl/src/m/shp/shpread.py
===================================================================
--- /issm/trunk-jpl/src/m/shp/shpread.py	(revision 25065)
+++ /issm/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: /issm/trunk-jpl/src/m/shp/shpwrite.m
===================================================================
--- /issm/trunk-jpl/src/m/shp/shpwrite.m	(revision 25064)
+++ /issm/trunk-jpl/src/m/shp/shpwrite.m	(revision 25065)
@@ -10,8 +10,4 @@
 %   See also SHPREAD
 
-
-%initialize number of profile
-count=0;
-
 contours=struct([]);
 for i=1:length(shp),
@@ -22,4 +18,6 @@
 	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;
@@ -27,5 +25,5 @@
 	contours(i).Y=shp(i).y;
 end
-	
+
 shapewrite(contours,filename);
 end
Index: /issm/trunk-jpl/src/m/shp/shpwrite.py
===================================================================
--- /issm/trunk-jpl/src/m/shp/shpwrite.py	(revision 25065)
+++ /issm/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()
+#}}}