Changeset 4175

Timestamp:

Jan 15, 2007, 12:10:20 PM (18 years ago)

Author:

duncan

Message:

added blocking on netcdf files & a fix in using verbose

Location:

anuga_core/source/anuga

Files:

• fit_interpolate/fit.py (modified) (2 diffs)
• fit_interpolate/test_fit.py (modified) (2 diffs)
• geospatial_data/geospatial_data.py (modified) (8 diffs)
• geospatial_data/test_geospatial_data.py (modified) (3 diffs)

anuga_core/source/anuga/fit_interpolate/fit.py

@@ -320 +320 @@
           
         """
-
         # use blocking to load in the point info
         if type(point_coordinates_or_filename) == types.StringType:
@@ -328 +327 @@
             for i,geo_block in  enumerate(Geospatial_data(filename,
                                               max_read_lines=max_read_lines,
-                                              load_file_now=False)):
+                                              load_file_now=False,
+                                              verbose=verbose)):
                 if verbose is True and 0 == i%200: # round every 5 minutes
                     print 'Block %i' %i

anuga_core/source/anuga/fit_interpolate/test_fit.py

@@ -241 +241 @@
         os.remove(fileName)
 
-    def test_fit_to_mesh(self):
-
-        a = [-1.0, 0.0]
-        b = [3.0, 4.0]
-        c = [4.0,1.0]
-        d = [-3.0, 2.0] #3
-        e = [-1.0,-2.0]
-        f = [1.0, -2.0] #5
-
-        vertices = [a, b, c, d,e,f]
-        triangles = [[0,1,3], [1,0,2], [0,4,5], [0,5,2]] #abd bac aef afc
-
-
-        fileName = tempfile.mktemp(".ddd")
-        file = open(fileName,"w")
-        file.write(" x,y, elevation \n\
--2.0, 2.0, 0.\n\
--1.0, 1.0, 0.\n\
-0.0, 2.0 , 2.\n\
-1.0, 1.0 , 2.\n\
- 2.0,  1.0 ,3. \n\
- 0.0,  0.0 , 0.\n\
- 1.0,  0.0 , 1.\n\
- 0.0,  -1.0, -1.\n\
- -0.2, -0.5, -0.7\n\
- -0.9, -1.5, -2.4\n\
- 0.5,  -1.9, -1.4\n\
- 3.0,  1.0 , 4.\n")
-        file.close()
-        
-        f = fit_to_mesh(vertices, triangles,fileName,
-                                alpha=0.0, max_read_lines=2)
-                        #use_cache=True, verbose=True)
-        answer = linear_function(vertices)
-        #print "f\n",f
-        #print "answer\n",answer
-        assert allclose(f, answer)
-        
-        os.remove(fileName)
    
     def test_fit_to_mesh_pts(self):
@@ -314 +275 @@
         geo.export_points_file(fileName_pts)
         f = fit_to_mesh(vertices, triangles,fileName_pts,
-                                alpha=0.0, max_read_lines=2)
+                                alpha=0.0, max_read_lines=2) #,
                         #use_cache=True, verbose=True)
         answer = linear_function(vertices)

anuga_core/source/anuga/geospatial_data/geospatial_data.py

@@ -13 +13 @@
 #from MA import tolist
 
+from Scientific.IO.NetCDF import NetCDFFile    
+    
 from anuga.utilities.numerical_tools import ensure_numeric
 from anuga.coordinate_transforms.geo_reference import Geo_reference, TitleError
@@ -586 +588 @@
             msg += 'Text file format is moving to comma seperated .txt files.'
             warn(msg, DeprecationWarning) 
-            error(msg, DeprecationWarning) 
 
         if (file_name[-4:] == ".xya"):
@@ -699 +700 @@
 
     def __iter__(self):
-        # read in the header and save the file pointer position
-
+        """
+        read in the header and save the file pointer position
+        """
+
+        from Scientific.IO.NetCDF import NetCDFFile
+        
         #FIXME - what to do if the file isn't there
 
-        #FIXME - give warning if the file format is .xya
-        if self.file_name[-4:]== ".xya" or self.file_name[-4:]== ".pts":
+        if self.file_name[-4:]== ".xya":
             #let's just read it all
             pass
+        elif self.file_name[-4:]== ".pts":
+            
+            # see if the file is there.  Throw a QUIET IO error if it isn't
+            fd = open(self.file_name,'r')
+            fd.close()
+    
+            #throws prints to screen if file not present
+            self.fid = NetCDFFile(self.file_name, 'r')
+            
+            self.blocking_georef, self.blocking_keys, self.last_row = \
+                     _read_pts_file_header(self.fid, self.verbose)
+            self.start_row=0
         else:
             file_pointer = open(self.file_name)
@@ -712 +728 @@
                          _read_csv_file_header(file_pointer)
 
-            if self.max_read_lines is None:
-                self.max_read_lines = MAX_READ_LINES
+        if self.max_read_lines is None:
+            self.max_read_lines = MAX_READ_LINES
         return self
     
     def next(self):
         # read a block, instanciate a new geospatial and return it
-        if self.file_name[-4:]== ".xya" or self.file_name[-4:]== ".pts":
+        if self.file_name[-4:]== ".xya" :
             if not hasattr(self,'finished_reading') or \
                    self.finished_reading is False:
@@ -728 +744 @@
                 self.finished_reading = False
                 
+        elif self.file_name[-4:]== ".pts":
+            if self.start_row == self.last_row:
+                # read the end of the file last iteration
+                # FIXME clean up, remove blocking atts eg
+                #self.max_read_lines
+                #self.blocking_georef,self.last_row
+                #self.start_row
+                ### self.blocking_keys
+                self.fid.close()
+                raise StopIteration
+            fin_row = self.start_row + self.max_read_lines
+            if fin_row > self.last_row:
+                fin_row = self.last_row
+
+            #call stuff
+            pointlist, att_dict, = \
+                   _read_pts_file_blocking( self.fid,
+                                            self.start_row,
+                                            fin_row,
+                                            self.blocking_keys
+                                            ) 
+            geo = Geospatial_data(pointlist, att_dict, self.blocking_georef)
+            self.start_row = fin_row
+            
         else:
             try:
                 pointlist, att_dict, self.file_pointer = \
                    _read_csv_file_blocking( self.file_pointer,
-                                         self.header[:],
-                                         max_read_lines=self.max_read_lines) 
+                                            self.header[:],
+                                            max_read_lines=self.max_read_lines,
+                                            verbose=self.verbose)
                 geo = Geospatial_data(pointlist, att_dict)
             except StopIteration:
@@ -761 +802 @@
     fid = NetCDFFile(file_name, 'r') 
     
-#    point_atts = {}  
-        # Get the variables
-#    point_atts['pointlist'] = array(fid.variables['points'])
     pointlist = array(fid.variables['points'])
     keys = fid.variables.keys()
@@ -780 +818 @@
         attributes[key] = array(fid.variables[key])
     
-#    point_atts['attributelist'] = attributes
     
     try:
         geo_reference = Geo_reference(NetCDFObject=fid)
-#        point_atts['geo_reference'] = geo_reference
     except AttributeError, e:
-        #geo_ref not compulsory 
-#        point_atts['geo_reference'] = None
         geo_reference = None
     
@@ -895 +929 @@
     return pointlist, att_dict,file_pointer 
 
+def _read_pts_file_header(fid, verbose=False):
+
+    """
+    Read the geo_reference of a .pts file
+    """
+    
+    keys = fid.variables.keys()
+    try:
+        keys.remove('points')
+    except IOError, e:       
+        fid.close()    
+        msg = 'Expected keyword "points" but could not find it'
+        raise IOError, msg
+    if verbose: print 'Got %d variables: %s' %(len(keys), keys)
+    
+    try:
+        geo_reference = Geo_reference(NetCDFObject=fid)
+    except AttributeError, e:
+        geo_reference = None
+
+    return geo_reference, keys, fid.dimensions['number_of_points']
+
+def _read_pts_file_blocking(fid, start_row, fin_row, keys):
+                            #verbose=False):
+    
+
+    """
+    Read the body of a .csv file.
+    header: The list header of the csv file, with the x and y labels.
+    """
+    
+    pointlist = array(fid.variables['points'][start_row:fin_row])
+    
+    attributes = {}
+    for key in keys:
+        attributes[key] = array(fid.variables[key][start_row:fin_row])
+
+    return pointlist, attributes
+    
+    
 def _read_xya_file(fd, delimiter):
     points = []

anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

@@ -16 +16 @@
 warnings.filterwarnings(action = 'ignore',
                         message='.xya format is deprecated.  Please use .txt.',
+                        category=DeprecationWarning)
+
+warnings.filterwarnings(action = 'ignore',
+                        message='The text file values must be ab',
                         category=DeprecationWarning)
 
@@ -974 +978 @@
                         [10.4])
            
-        os.remove(fileName)             
+        os.remove(fileName)         
         
     def test_load_csv_bad(self):
@@ -1004 +1008 @@
             raise msg
         os.remove(fileName)
-        
-    def depreciated_test_export_xya_file(self):
+
+        
+    def test_load_pts_blocking(self):
+        """ test_load_csv(self):
+        space delimited
+        """
+        import os
+       
+        fileName = tempfile.mktemp(".txt")
+        file = open(fileName,"w")
+        file.write(" x,y, elevation ,  speed \n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0 ,10.4, 40.0\n")
+        file.close()
+
+        pts_file = tempfile.mktemp(".pts")
+        
+        convert = Geospatial_data(fileName)
+        convert.export_points_file(pts_file)
+        results = Geospatial_data(pts_file, max_read_lines=2)
+
+        assert allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],
+                                                    [1.0, 0.0]])
+        assert allclose(results.get_attributes(attribute_name='elevation'),
+                        [10.0, 0.0, 10.4])
+        assert allclose(results.get_attributes(attribute_name='speed'),
+                        [0.0, 10.0, 40.0])
+
+        # Blocking
+        geo_list = []
+        for i in results:
+            geo_list.append(i) 
+        assert allclose(geo_list[0].get_data_points(),
+                        [[1.0, 0.0],[0.0, 1.0]])
+        assert allclose(geo_list[0].get_attributes(attribute_name='elevation'),
+                        [10.0, 0.0])
+        assert allclose(geo_list[1].get_data_points(),
+                        [[1.0, 0.0]])        
+        assert allclose(geo_list[1].get_attributes(attribute_name='elevation'),
+                        [10.4])
+           
+        os.remove(fileName)  
+        os.remove(pts_file)               
+
+        
+    def test_export_xya_file(self):
+        # depreciated since it's testing using geo refs in a text file
 #        dict = {}
         att_dict = {}