source: issm/trunk/src/m/exp/expcoarsen.py@ 21341

import os.path
import numpy as  np
from collections import OrderedDict
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=np.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=np.sqrt((x[j]-x[j+1])**2+(y[j]-y[j+1])**2)
                        if distance<resolution and j<numpoints-2:   #do not remove last point
                                x=np.delete(x,j+1,0)
                                y=np.delete(y,j+1,0)
                                numpoints=numpoints-1
                        else:
                                division=int(np.floor(distance/resolution)+1)
                                if division>=2:
                                        xi=np.linspace(x[j],x[j+1],division)
                                        yi=np.linspace(y[j],y[j+1],division)
                                        
                                        x=np.hstack((x[0:j+1],xi[1:-1],x[j+1:]))
                                        y=np.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 np.size(x)>1:
                        #keep the (x,y) contour arond
                        newcontour=OrderedDict()
                        newcontour['nods']=np.size(x)
                        newcontour['density']=contour['density']
                        newcontour['x']=x
                        newcontour['y']=y
                        newcontours.append(newcontour)

        #write output
        expwrite(newcontours,newfile)