Index: /issm/trunk-jpl/src/m/exp/expcoarsen.py
===================================================================
--- /issm/trunk-jpl/src/m/exp/expcoarsen.py	(revision 16098)
+++ /issm/trunk-jpl/src/m/exp/expcoarsen.py	(revision 16098)
@@ -0,0 +1,79 @@
+import os.path
+import numpy as npy
+from collections import OrderedDict
+from MatlabFuncs import *
+from expread import expread
+from expwrite import expwrite
+
+def expcoarsen(newfile,oldfile,resolution):
+	""" 
+	EXPCOARSEN - coarsen an exp contour
+
+	This routine read an Argus file and remove points with respect to
+	the resolution (in meters) given in input. 
+
+	Usage:
+	  expcoarsen(newfile,oldfile,resolution)
+
+	Example:
+	   expcoarsen('DomainOutline.exp','Antarctica.exp',4000)
+	""" 
+
+	#Some checks
+	if not os.path.exists(oldfile):
+		raise OSError("expcoarsen error message: file '%s' not found!" % oldfile)
+	if os.path.exists(newfile):
+		choice=raw_input('A file ' + newfile + ' already exists, do you want to modify it? (y/n)')
+		if choice not in 'y': 
+			print('no modification done ... exiting')
+			return 0
+
+	#Get exp oldfile
+	contours=expread(oldfile)
+	newcontours=[]
+
+	for contour in  contours:
+		
+		numpoints=npy.size(contour['x'])
+
+		j=0
+		x=contour['x']
+		y=contour['y']
+
+		#stop if we have reached end of profile (always keep the last point)
+		while j<numpoints-1:
+
+			#see whether we keep this point or not
+			distance=npy.sqrt((x[j]-x[j+1])**2+(y[j]-y[j+1])**2)
+			print j,distance,numpoints
+			if distance<resolution and j<numpoints-2:   #do not remove last point
+				x=npy.delete(x,j+1,0)
+				y=npy.delete(y,j+1,0)
+				numpoints=numpoints-1
+			else:
+				division=int(npy.floor(distance/resolution)+1)
+				if division>=2:
+					xi=npy.linspace(x[j],x[j+1],division)
+					yi=npy.linspace(y[j],y[j+1],division)
+					
+					x=npy.hstack((x[0:j+1],xi[1:-1],x[j+1:]))
+					y=npy.hstack((y[0:j+1],yi[1:-1],y[j+1:]))
+
+					#update current point
+					j=j+1+division-2
+					numpoints=numpoints+division-2
+				else:
+					#update current point
+					j=j+1
+		
+		if npy.size(x)>1:
+			#keep the (x,y) contour arond
+			newcontour=OrderedDict()
+			newcontour['nods']=npy.size(x)
+			newcontour['density']=contour['density']
+			newcontour['x']=x
+			newcontour['y']=y
+			newcontours.append(newcontour)
+
+	#write output
+	expwrite(newcontours,newfile)