Changeset 2624

Timestamp:

Mar 29, 2006, 8:42:34 AM (18 years ago)

Author:

nick

Message:

fixed geospatial_data to not use point_dict in export and updated onslow run to use new geospatial_data methods

Files:

• inundation/geospatial_data/geospatial_data.py (modified) (16 diffs)
• inundation/geospatial_data/test_geospatial_data.py (modified) (16 diffs)
• production/onslow_2006/project.py (modified) (1 diff)
• production/onslow_2006/run_onslow.py (modified) (8 diffs)

inundation/geospatial_data/geospatial_data.py

@@ -20 +20 @@
                  geo_reference = None,
                  default_attribute_name = None,
-                 file_name = None):
-         
+                 file_name = None,
+                 absolute = True):
 
         """Create instance from data points and associated attributes
@@ -47 +47 @@
         file_name: Name of input file.....
         
+        absolute: is the data point absolute or relative to a xll and yll 
+        in a zone
+        
         """
 
@@ -56 +59 @@
             self.set_geo_reference(geo_reference)
             self.set_default_attribute_name(default_attribute_name)
+            self.set_absolute(absolute)
             
         
@@ -152 +156 @@
                 self.file_name = file_name
                 print 'file name from set', self.file_name
-
+    
+    def set_absolute(self, absolute = True):
+        if absolute is False:
+            self.absolute = False
+        else:
+            self.absolute = True
+    
+    '''
     def set_self_from_file(self, file_name = None):
         # check if file_name is a string and also that the file exists
@@ -162 +173 @@
             self = self.import_points_file(file_name)
         return self
-            
+        '''    
 
     def get_geo_reference(self):
@@ -202 +213 @@
         return self.attributes[attribute_name]
 
+    def get_all_attributes(self):
+        """
+        Return values for all attributes.
+        """
+
+        return self.attributes
 
     def __add__(self, other):
-        """ 
+        '''
         returns the addition of 2 geospatical objects,
         objects are concatencated to the end of each other
@@ -210 +227 @@
         NOTE: doesn't add if objects contain different 
         attributes  
-        """
+        
+        Always return relative points!
+        '''
 
 
@@ -237 +256 @@
 
         xll = yll = 0. 
-        relative_points1 = self.get_data_points(absolute=False)
-        relative_points2 = other.get_data_points(absolute=False)
-            
-        #print 'R1', relative_points1
-        #print 'R2', relative_points2        
+        
+        relative_points1 = self.get_data_points(absolute = False)
+        relative_points2 = other.get_data_points(absolute = False)
 
         
         new_relative_points1 = new_geo_ref.change_points_geo_ref(relative_points1, geo_ref1)
         new_relative_points2 = new_geo_ref.change_points_geo_ref(relative_points2, geo_ref2)
-
-        #print 'NR1', new_relative_points1
-        #print 'NR2', new_relative_points2        
         
         #Now both point sets are relative to new_geo_ref and zones have been reconciled
@@ -256 +270 @@
         new_points = concatenate((new_relative_points1,
                                   new_relative_points2),
-                                 axis = 0)
+                                  axis = 0)
+ #       print 'new points:', new_points
 
       
@@ -277 +292 @@
         return Geospatial_data(new_points,
                                new_attributes,
-                               new_geo_ref)
-
-
-
+                               new_geo_ref,
+                               absolute = False)
+
+
+    '''
     def xxx__add__(self, other):
         """ 
@@ -339 +355 @@
             raise msg
 
-      
-    
-###
-#  IMPORT/EXPORT POINTS FILES
-###
+    '''   
+    
+    ###
+    #  IMPORT/EXPORT POINTS FILES
+    ###
+    def new_import_points_file(self, ofile, delimiter = None, verbose = False):
+        """ load an .xya or .pts file
+        Note: will throw an IOError if it can't load the file.
+        Catch these!
+        """
+        
+        all_data = {}
+        if ofile[-4:]== ".xya":
+            try:
+                if delimiter == None:
+                    try:
+                        fd = open(ofile)
+                        all_data = _read_xya_file(fd, ',')
+                    except SyntaxError:
+                        fd.close()
+                        fd = open(ofile)
+                        all_data = _read_xya_file(fd, ' ')
+                else:
+                    fd = open(ofile)
+                    all_data = _read_xya_file(fd, delimiter)
+                fd.close()
+            except (IndexError,ValueError,SyntaxError):
+                fd.close()    
+                msg = 'Could not open file %s ' %ofile
+                raise IOError, msg
+            except IOError, e:
+                # Catch this to add an error message
+                msg = 'Could not open file or incorrect file format %s:%s' %(ofile, e)
+                raise IOError, msg
+                
+        elif ofile[-4:]== ".pts":
+            try:
+    #            print 'hi from import_points_file'
+                all_data = _read_pts_file(ofile, verbose)
+    #            print 'hi1 from import_points_file', all_data
+            except IOError, e:    
+                msg = 'Could not open file %s ' %ofile
+                raise IOError, msg        
+        else:      
+            msg = 'Extension %s is unknown' %ofile[-4:]
+            raise IOError, msg
+    
+        return all_data
+    
+    def new_export_points_file(self, ofile):
+        """
+        write a points file, ofile, as a text (.xya) or binary (.pts) file
+    
+        ofile is the file name, including the extension
+    
+        The point_dict is defined at the top of this file.
+        """
+        #this was done for all keys in the mesh file.
+        #if not mesh_dict.has_key('points'):
+        #    mesh_dict['points'] = []
+        if (ofile[-4:] == ".xya"):
+            _new_write_xya_file(ofile, self.get_data_points(self.absolute), 
+                                   self.get_all_attributes(),
+                                   self.get_geo_reference())
+        elif (ofile[-4:] == ".pts"):
+            _new_write_pts_file(ofile, self.get_data_points(self.absolute), 
+                                   self.get_all_attributes(),
+                                   self.get_geo_reference())
+        else:
+            msg = 'Unknown file type %s ' %ofile
+            raise IOError, msg 
+    
 
 def import_points_file( ofile, delimiter = None, verbose = False):
@@ -388 +471 @@
 
     return all_data
-
+'''
 def export_points_file(ofile, point_dict):
     """
@@ -407 +490 @@
         msg = 'Unknown file type %s ' %ofile
         raise IOError, msg 
-
+'''
 def _read_pts_file(file_name, verbose = False):
     """Read .pts NetCDF file
@@ -520 +603 @@
 
 
-
+def _new_write_pts_file(file_name, write_data_points,
+                                   write_attributes, 
+                                   write_geo_reference = None):
+    """Write .pts NetCDF file   
+
+    WARNING: This function mangles the point_atts data structure
+    """
+    #FIXME: (DSG)This format has issues.
+    # There can't be an attribute called points 
+    # consider format change
+    # method changed by NB not sure if above statement is correct
+
+    ''' should create new test for this
+    legal_keys = ['pointlist', 'attributelist', 'geo_reference']
+    for key in point_atts.keys():
+        msg = 'Key %s is illegal. Valid keys are %s' %(key, legal_keys) 
+        assert key in legal_keys, msg
+    '''    
+    from Scientific.IO.NetCDF import NetCDFFile
+    # NetCDF file definition
+    outfile = NetCDFFile(file_name, 'w')
+    
+    #Create new file
+    outfile.institution = 'Geoscience Australia'
+    outfile.description = 'NetCDF format for compact and portable storage ' +\
+                      'of spatial point data'
+    
+    # dimension definitions
+    shape = write_data_points.shape[0]
+    outfile.createDimension('number_of_points', shape)  
+    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
+    
+    # variable definition
+    outfile.createVariable('points', Float, ('number_of_points',
+                                             'number_of_dimensions'))
+
+    #create variables  
+    outfile.variables['points'][:] = write_data_points #.astype(Float32)
+
+    for key in write_attributes.keys():
+        outfile.createVariable(key, Float, ('number_of_points',))
+        outfile.variables[key][:] = write_attributes[key] #.astype(Float32)
+        
+    if write_geo_reference is not None:
+        write_geo_reference.write_NetCDF(outfile)
+        
+    outfile.close() 
+  
+
+
+def _new_write_xya_file( file_name, write_data_points,
+                                    write_attributes, 
+                                    write_geo_reference = None, 
+                                    delimiter = ','):
+    """
+    export a file, ofile, with the xya format
+    
+    """
+    points = write_data_points 
+    pointattributes = write_attributes
+    
+    fd = open(file_name,'w')
+    titlelist = ""
+    for title in pointattributes.keys():
+        titlelist = titlelist + title + delimiter
+    titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
+    fd.write(titlelist+"\n")
+    #<vertex #> <x> <y> [attributes]
+    for i,vert in enumerate( points):
+        
+        attlist = ","
+        for att in pointattributes.keys():
+            attlist = attlist + str(pointattributes[att][i])+ delimiter
+        attlist = attlist[0:-len(delimiter)] # remove the last delimiter
+        attlist.strip()
+        fd.write( str(vert[0]) + delimiter
+                  + str(vert[1])
+                  + attlist + "\n")
+    '''    
+    # geo_reference info
+    if xya_dict.has_key('geo_reference') and \
+           not xya_dict['geo_reference'] is None:
+        xya_dict['geo_reference'].write_ASCII(fd)
+    '''
+    if  write_geo_reference is not None:
+        write_geo_reference.write_ASCII(fd)
+    fd.close()
+'''
 def _write_pts_file(file_name, point_atts):
     """Write .pts NetCDF file   
@@ -600 +770 @@
         xya_dict['geo_reference'].write_ASCII(fd)
     fd.close()
-    
+    '''
 def _point_atts2array(point_atts):
     point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
@@ -682 +852 @@
     
     #FIXME remove dependance on points to dict in export only!
-    G_points_dict = geospatial_data2points_dictionary(G)
-    export_points_file(results_file, G_points_dict)
+#    G_points_dict = geospatial_data2points_dictionary(G)
+#    export_points_file(results_file, G_points_dict)
+
+#    G_points_dict = geospatial_data2points_dictionary(G)
+
+    G.new_export_points_file(results_file)
     
     

inundation/geospatial_data/test_geospatial_data.py

@@ -301 +301 @@
         points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
         attributes = [2, 4, 5, 76]
-
+        '''
         # Use old pointsdict format
         pointsdict = {'pointlist': points,
                       'attributelist': {'att1': attributes,
                                         'att2': array(attributes) + 1}} 
-     
+        '''
+        att_dict = {'att1': attributes,
+                    'att2': array(attributes) +1}
+                    
         # Create points as an xya file 
         FN = 'test_points.xya'
-        export_points_file(FN, pointsdict)
-
+        G1 = Geospatial_data(points, att_dict)
+        G1.new_export_points_file(FN)
+#        G1.new_export_points_file(ofile)
 
         #Create object from file
@@ -397 +401 @@
                         'imaginary file did not raise error!')
   
-    def Xtest_read_write_points_file_bad(self):
-        #not used yet as i'm uncertain if "tri_dict" is nesscessary to geospatial_data
-        dict = self.tri_dict.copy()
-        fileName = tempfile.mktemp(".xxx")
-        try:
-            export_points_file(fileName,dict)
-        except IOError:
-            pass
-        else:
-            self.failUnless(0 ==1,
-                        'bad points file extension did not raise error!')
-
     def test_read_write_points_file_bad2(self):
-        dict = {}
+#        dict = {}
         att_dict = {}
-        dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
         att_dict['elevation'] = array([10.0, 0.0, 10.4])
         att_dict['brightness'] = array([10.0, 0.0, 10.4])
-        dict['attributelist'] = att_dict
-        dict['geo_reference'] = Geo_reference(56,1.9,1.9)
-        try:
-            export_points_file("_???/yeah.xya",dict)
+#        dict['attributelist'] = att_dict
+        geo_reference = Geo_reference(56,1.9,1.9)
+        
+        G = Geospatial_data(pointlist, att_dict, geo_reference)
+        
+        try:
+#            export_points_file("_???/yeah.xya",dict)
+            G.new_export_points_file("_???/yeah.xya")
+            
         except IOError:
             pass
         else:
-            self.failUnless(0 ==1,
+            self.failUnless(0 == 1,
                         'bad points file extension did not raise error!')
                    
@@ -573 +570 @@
         
     def test_export_xya_file(self):
-        dict = {}
+#        dict = {}
         att_dict = {}
-        dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
         att_dict['elevation'] = array([10.0, 0.0, 10.4])
         att_dict['brightness'] = array([10.0, 0.0, 10.4])
-        dict['attributelist'] = att_dict
-        dict['geo_reference'] = Geo_reference(56,1.9,1.9)
-        
-        
-        fileName = tempfile.mktemp(".xya")
-        export_points_file(fileName, dict)
+#        dict['attributelist'] = att_dict
+        geo_reference = Geo_reference(56,1.9,1.9)
+        
+        
+        fileName = tempfile.mktemp(".xya")
+        G = Geospatial_data( pointlist, att_dict, geo_reference, absolute = False)
+        G.new_export_points_file(fileName)
         dict2 = import_points_file(fileName)
         #print "fileName",fileName 
@@ -594 +592 @@
         assert allclose(dict2['attributelist']['brightness'], answer)
         #print "dict2['geo_reference']",dict2['geo_reference'] 
-        self.failUnless(dict['geo_reference'] == dict2['geo_reference'],
+        self.failUnless(geo_reference == dict2['geo_reference'],
                          'test_writepts failed. Test geo_reference')
 
     def test_export_xya_file2(self):
-        dict = {}
+#        dict = {}
         att_dict = {}
-        dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
         att_dict['elevation'] = array([10.0, 0.0, 10.4])
         att_dict['brightness'] = array([10.0, 0.0, 10.4])
-        dict['attributelist'] = att_dict
-        
-        
-        fileName = tempfile.mktemp(".xya")
-        export_points_file(fileName, dict)
+#        dict['attributelist'] = att_dict
+        
+        fileName = tempfile.mktemp(".xya")
+        G = Geospatial_data(pointlist, att_dict)
+        G.new_export_points_file(fileName)
         dict2 = import_points_file(fileName)
         #print "fileName",fileName 
@@ -616 +614 @@
         assert allclose(dict2['attributelist']['elevation'], [10.0, 0.0, 10.4])
         answer = [10.0, 0.0, 10.4]
+        assert allclose(dict2['attributelist']['brightness'], answer)
+
+    def test_new_export_xya_file2(self):
+        dict = {}
+        att_dict = {}
+        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        att_dict['elevation'] = array([10.1, 0.0, 10.4])
+        att_dict['brightness'] = array([10.0, 1.0, 10.4])
+        
+        fileName = tempfile.mktemp(".xya")
+        
+        G = Geospatial_data(pointlist, att_dict)
+        
+        G.new_export_points_file(fileName)
+
+        dict2 = import_points_file(fileName)
+
+        os.remove(fileName)
+        
+        assert allclose(dict2['pointlist'],[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        assert allclose(dict2['attributelist']['elevation'], [10.1, 0.0, 10.4])
+        answer = [10.0, 1.0, 10.4]
+        assert allclose(dict2['attributelist']['brightness'], answer)
+
+    def test_new_export_pts_file(self):
+        att_dict = {}
+        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        att_dict['elevation'] = array([10.1, 0.0, 10.4])
+        att_dict['brightness'] = array([10.0, 1.0, 10.4])
+        
+        fileName = tempfile.mktemp(".pts")
+        
+        G = Geospatial_data(pointlist, att_dict)
+        
+        G.new_export_points_file(fileName)
+
+        dict2 = import_points_file(fileName)
+
+        os.remove(fileName)
+        
+        assert allclose(dict2['pointlist'],[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        assert allclose(dict2['attributelist']['elevation'], [10.1, 0.0, 10.4])
+        answer = [10.0, 1.0, 10.4]
         assert allclose(dict2['attributelist']['brightness'], answer)
 
@@ -658 +699 @@
         
     def test_writepts(self):
-        dict = {}
+#        dict = {}
         att_dict = {}
-        dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        pointlist = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
         att_dict['elevation'] = array([10.0, 0.0, 10.4])
         att_dict['brightness'] = array([10.0, 0.0, 10.4])
-        dict['attributelist'] = att_dict
-        dict['geo_reference'] = Geo_reference(56,1.9,1.9)
-        
+#        dict['attributelist'] = att_dict
+        geo_reference = Geo_reference(56,1.9,1.9)
         
         fileName = tempfile.mktemp(".pts")
-        export_points_file(fileName, dict)
+        
+        G = Geospatial_data(pointlist, att_dict, geo_reference, absolute = False)
+        
+        G.new_export_points_file(fileName)
+        
         dict2 = import_points_file(fileName)
-        #print "fileName",fileName 
-        os.remove(fileName)
-        #print "dict2",dict2 
-        
+
+        os.remove(fileName)
+
         assert allclose(dict2['pointlist'],[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
         assert allclose(dict2['attributelist']['elevation'], [10.0, 0.0, 10.4])
@@ -682 +725 @@
         #print "dict2['geo_reference']",dict2['geo_reference']
         
-        self.failUnless(dict['geo_reference'] == dict2['geo_reference'],
+        self.failUnless(geo_reference == dict2['geo_reference'],
                          'test_writepts failed. Test geo_reference')
         
@@ -705 +748 @@
         
         from Scientific.IO.NetCDF import NetCDFFile
-
 
 #        fileName = tempfile.mktemp(".pts")
@@ -738 +780 @@
         assert allclose(G.get_geo_reference().get_xllcorner(), 0.0)
         assert allclose(G.get_geo_reference().get_yllcorner(), 0.0)
-        
 
         assert allclose(G.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
@@ -749 +790 @@
         
         from Scientific.IO.NetCDF import NetCDFFile
-
 
         FN = 'test_points.pts'
@@ -761 +801 @@
         geo_reference.write_NetCDF(outfile)
 
-        
         # dimension definitions
         outfile.createDimension('number_of_points', 3)    
@@ -775 +814 @@
         elevation = outfile.variables['elevation']
 
- 
         points[0, :] = [1.0,0.0]
         elevation[0] = 10.0 
@@ -787 +825 @@
         G = Geospatial_data(file_name = FN)
 
-
-
         assert allclose(G.get_geo_reference().get_xllcorner(), xll)
         assert allclose(G.get_geo_reference().get_yllcorner(), yll)
-        
 
         assert allclose(G.get_data_points(), [[1.0+xll, 0.0+yll],
@@ -801 +836 @@
 
         
-    def xtest_add_points_files(self):
+    def test_add_points_files(self):
         '''adds an xya and pts files and checks resulting file is correct
         '''
@@ -808 +843 @@
         
         # create files
-        dict1 = {}
+#        dict1 = {}
         att_dict1 = {}
-        dict1['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        pointlist1 = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
         att_dict1['elevation'] = array([-10.0, 0.0, 10.4])
         att_dict1['brightness'] = array([10.0, 0.0, 10.4])
-        dict1['attributelist'] = att_dict1
-        dict1['geo_reference'] = Geo_reference(56, 1.9, 1.9)
-        
-        dict2 = {}
+#        dict1['attributelist'] = att_dict1
+        geo_reference1 = Geo_reference(56, 2.0, 1.0)
+        
+#        dict2 = {}
         att_dict2 = {}
-        dict2['pointlist'] = array([[2.0, 1.0],[1.0, 2.0],[2.0, 1.0]])
+        pointlist2 = array([[2.0, 1.0],[1.0, 2.0],[2.0, 1.0]])
         att_dict2['elevation'] = array([1.0, 15.0, 1.4])
         att_dict2['brightness'] = array([14.0, 1.0, -12.4])
-        dict2['attributelist'] = att_dict2
-        dict2['geo_reference'] = Geo_reference(56, 1.0, 1.0) #FIXME (Ole): I changed this, why does it still pass
+#        dict2['attributelist'] = att_dict2
+        geo_reference2 = Geo_reference(56, 1.0, 2.0) #FIXME (Ole): I changed this, why does it still pass
         #OK - it fails now, after the fix revealed by test_create_from_pts_file_with_geo')
+
+        G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1)
+        G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2)
         
         fileName1 = tempfile.mktemp(".xya")
         fileName2 = tempfile.mktemp(".pts")
 
-        export_points_file(fileName1, dict1)
-        export_points_file(fileName2, dict2)
+#        G1.new_export_points_file(fileName1)
+#        G2.new_export_points_file(fileName2)
         
         results_file = 'resulting_points.pts'
         
         # add files
-        add_points_files(fileName1, fileName2, results_file )
+        
+        G = G1 + G2
+        
+        G.new_export_points_file(results_file)
+#        add_points_files(fileName1, fileName2, results_file )
         
         #read results
         results_dict = import_points_file(results_file)
+#        print 'results points:', results_dict['pointlist']
+#        print 'results geo_ref:', results_dict['geo_reference']
         
         assert allclose(results_dict['pointlist'],
-                        [[1.0, 0.0],[0.0, 1.0],
-                         [1.0, 0.0],[2.0, 1.0],
-                         [1.0, 2.0],[2.0, 1.0]])
+                        [[2.0, 0.0],[1.0, 1.0],
+                         [2.0, 0.0],[2.0, 2.0],
+                         [1.0, 3.0],[2.0, 2.0]])
         assert allclose(results_dict['attributelist']['elevation'],
                         [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4])
@@ -849 +893 @@
         assert allclose(results_dict['attributelist']['brightness'], answer)
         
-        self.failUnless(results_dict['geo_reference'] == dict2['geo_reference'],
+        self.failUnless(results_dict['geo_reference'] == geo_reference1,
                          'test_writepts failed. Test geo_reference')
                          
-        os.remove(fileName1)
-        os.remove(fileName2)
+#        os.remove(fileName1)
+#        os.remove(fileName2)
         os.remove(results_file)
 

production/onslow_2006/project.py

@@ -47 +47 @@
 polygondir = home+sep+scenario_dir_name+sep+'polygons'+sep
 boundarydir = home+sep+scenario_dir_name+sep+'boundaries'+sep
+print'bound', boundarydir
+
+# boundary source data
+MOST_dir = 'f:'+sep+'3'+sep+'ehn'+sep+'users'+sep+'davidb'+sep+'tsunami'+sep+'WA_project'+sep+'SU-AU_90'+sep+'most_2'+sep+'detailed'+sep
 
 codedir = getcwd()+sep

production/onslow_2006/run_onslow.py

@@ -55 +55 @@
 source_dir = project.boundarydir
 
-# creates copy of code in output dir
+# creates copy of code in output dir if dir doesn't exist
 if access(project.outputdir,F_OK) == 0 :
     mkdir (project.outputdir)
 copy (project.codedirname, project.outputdir + project.codename)
 copy (project.codedir + 'run_onslow.py', project.outputdir + 'run_onslow.py')
-
+print'hello'
+print' most file', project.MOST_dir + project.boundary_basename+'_ha.nc'
+if access(project.MOST_dir + project.boundary_basename+'_ha.nc',F_OK) == 1 :
+    print' most file', project.MOST_dir + project.boundary_basename
 
 '''
@@ -79 +82 @@
 
         '''
+'''
 print'create G1'
 G1 = Geospatial_data(file_name = project.offshore_dem_name + '.xya')
@@ -91 +95 @@
 G.new_export_points_file(project.combined_dem_name + '.pts')
 
-'''
-add_points_files(
-                  project.offshore_dem_name + '.xya',
-                  project.onshore_dem_name + '.pts',
-                  project.combined_dem_name + '.pts')
 '''
 #-------------------------------------------------------------------------------                                 
@@ -106 +105 @@
 
 # original
-interior_res = 5000
+interior_res = 50000
 interior_regions = [[project.poly_onslow, interior_res],
                     [project.poly_thevenard, interior_res],
@@ -119 +118 @@
                              'left': [2], 'bottom': [3],
                              'bottomright': [4], 'topright': [5]},
-           'maximum_triangle_area': 100000,
+           'maximum_triangle_area': 1000000,
            'filename': meshname,           
            'interior_regions': interior_regions},
@@ -167 +166 @@
 domain.set_quantity('friction', 0.0) 
 print 'hi and file',project.combined_dem_name + '.pts'
+
 domain.set_quantity('elevation', 
 #                    0.
@@ -176 +176 @@
                     alpha = 0.1
                     )
+
 print 'hi1'
 
@@ -190 +191 @@
 
 cache(ferret2sww,
-      (source_dir + project.boundary_basename,
+#      (source_dir + project.boundary_basename,
+#       source_dir + project.boundary_basename), 
+      (project.MOST_dir + project.boundary_basename,
        source_dir + project.boundary_basename), 
       {'verbose': True,