Changeset 2570 for inundation/geospatial_data/geospatial_data.py

Timestamp:

Mar 21, 2006, 3:51:39 PM (17 years ago)

Author:

nick

Message:

Added functionality to geospatial_data.py including import and exporting pts and xya files

File:

• inundation/geospatial_data/geospatial_data.py (modified) (7 diffs)

inundation/geospatial_data/geospatial_data.py

@@ -7 +7 @@
 from utilities.numerical_tools import ensure_numeric
 
-from Numeric import concatenate
+from Numeric import concatenate, array, Float, shape
 
 from coordinate_transforms.geo_reference import Geo_reference
+
+from os import access, F_OK, R_OK
         
 class Geospatial_data:
 
     def __init__(self,
-                 data_points,
+                 data_points = None,
                  attributes = None,
                  geo_reference = None,
-                 default_attribute_name = None):
+                 default_attribute_name = None,
+                 file_name = None):
+         
 
         """Create instance from data points and associated attributes
@@ -41 +45 @@
         If None, the default is the 'first'
         
+        file_name: Name of input file.....
+        
         """
 
-        self.data_points = ensure_numeric(data_points)
-        self.set_attributes(attributes)
-        self.set_geo_reference(geo_reference)
-        self.set_default_attribute_name(default_attribute_name)
-        
-
-
+        if file_name is None:
+            self.file_name = None
+            self.check_data_points(data_points)
+
+            self.set_attributes(attributes)
+            self.set_geo_reference(geo_reference)
+            self.set_default_attribute_name(default_attribute_name)
+            
+        
+        else:
+#            print 'from init file name:', file_name
+            if access(file_name, F_OK) == 0 :
+                msg = 'File %s does not exist or is not accessable' %file_name
+                raise msg
+            else:
+                data = {}    
+# watch for case where file name and points, attributes etc are provided!!
+# if file name then all provided info will be removed!
+                self.file_name = file_name
+                data = import_points_file(self.file_name)
+                
+#                print 'data: point in init', data['pointlist']
+#                print 'attrib: point in init', data['attributelist']
+#                print 'geo_ref: point in init', data['geo_reference']
+                
+                data_points = data['pointlist']
+                attributes = data['attributelist']
+                get_reference = data['geo_reference']
+                
+                self.check_data_points(data_points)
+                self.set_attributes(attributes)
+                self.set_geo_reference(geo_reference)
+                self.set_default_attribute_name(default_attribute_name)
+
+        
+    def check_data_points(self, data_points):
+        """Checks data points
+        """
+        
+        if data_points is None:
+            self.data_points = None
+#           FIXME: should throw an error if bad file name and no data!
+            print 'There is no data or files to read for data!!'
+            msg = 'file name %s does not exist in data set and the ' %self.file_name
+#           some assert 
+            
+        else:
+            self.data_points = ensure_numeric(data_points)
 
     def set_attributes(self, attributes):
@@ -95 +142 @@
         self.default_attribute_name = default_attribute_name
 
+    def set_file_name(self, file_name = None):
+        if file_name is None:
+            self.file_name = None
+        else:
+            if access(file_name, F_OK) == 0 :
+                msg = 'Geospatial_data object needs have either data_points \n'
+                msg += ' or a file with data points'
+                raise msg
+            else:    
+                self.file_name = file_name
+                print 'file name from set', self.file_name
+
+    def set_self_from_file(self, file_name = None):
+        # check if file_name is a string and also that the file exists
+        if file_name is None:
+            self.file_name = None
+#        else:
+#            if access(self.file_name, F_OK) == 0 :
+#                msg = 'Geospatial_data object needs have either data_points \n'
+#                msg += ' or a file with data points'
+#                raise msg
+        else:    
+            self.file_name = file_name
+            self = self.import_points_file(file_name)
+#            print 'data points from set self file name', self.data_points
+#            print 'attributes from set self file name', self.attributes
+#            print 'geo_ref from set self file name', self.geo_reference
+        return self
+            
 
     def get_geo_reference(self):
@@ -156 +232 @@
         geo_ref = Geo_reference(zone1, xll, yll)
         
-        if zone1 == zone2:  
+        if zone1 == zone2:
             a_rel_points = self.get_data_points()
             b_rel_points = other.get_data_points()
@@ -184 +260 @@
             ZONE to allow addition.'
             raise msg
+
       
+    
+    ###
+    #  IMPORT/EXPORT POINTS FILES
+    ###
+    '''
+    def export_points_file(ofile, point_dict):
+        """
+        write a points file, ofile, as a binary (.xya) or text (.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"):
+            _write_xya_file(ofile, point_dict)
+        elif (ofile[-4:] == ".pts"):
+            _write_pts_file(ofile, point_dict)
+        else:
+            msg = 'Unknown file type %s ' %ofile
+        raise IOError, msg 
+
+    '''       
+def import_points_file( 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:
+            # Catch this to add an error message
+            msg = 'Could not open file %s ' %ofile
+            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 export_points_file(ofile, point_dict):
+    """
+    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"):
+        _write_xya_file(ofile, point_dict)
+    elif (ofile[-4:] == ".pts"):
+        _write_pts_file(ofile, point_dict)
+    else:
+        msg = 'Unknown file type %s ' %ofile
+        raise IOError, msg 
+
+def _read_pts_file(file_name, verbose = False):
+    """Read .pts NetCDF file
+    
+    Return a dic of array of points, and dic of array of attribute
+    eg
+    dic['points'] = [[1.0,2.0],[3.0,5.0]]
+    dic['attributelist']['elevation'] = [[7.0,5.0]
+    """    
+    #FIXME: (DSG) This format has issues.
+    # There can't be an attribute called points 
+    # consider format change
+
+#    print 'hi for read points file'
+    from Scientific.IO.NetCDF import NetCDFFile
+    
+    if verbose: print 'Reading ', file_name
+    
+        
+    # see if the file is there.  Throw a QUIET IO error if it isn't
+    fd = open(file_name,'r')
+    fd.close()
+    
+    #throws prints to screen if file not present
+    fid = NetCDFFile(file_name, 'r') 
+    
+    point_atts = {}  
+        # Get the variables
+    point_atts['pointlist'] = array(fid.variables['points'])
+    keys = fid.variables.keys()
+    if verbose: print 'Got %d variables: %s' %(len(keys), keys)
+    try:
+        keys.remove('points')
+    except IOError, e:       
+        fid.close()    
+        msg = 'Expected keyword "points" but could not find it'
+        raise IOError, msg
+    
+    attributes = {}
+    for key in keys:
+        if verbose: print "reading attribute '%s'" %key
+            
+        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
+    
+    fid.close()
+    return point_atts
+
+
+def _read_xya_file( fd, delimiter):
+#    print 'hello from read xya data'
+    points = []
+    pointattributes = []
+    title = fd.readline()
+    att_names = clean_line(title,delimiter)
+    
+    att_dict = {}
+    line = fd.readline()
+    numbers = clean_line(line,delimiter)
+    while len(numbers) > 1:
+        if numbers != []:
+            try:
+                x = float(numbers[0])
+                y = float(numbers[1])
+                points.append([x,y])
+                numbers.pop(0)
+                numbers.pop(0)
+                #attributes = []
+                #print "att_names",att_names
+                #print "numbers",numbers 
+                if len(att_names) != len(numbers):
+                    fd.close()
+                    # It might not be a problem with the title
+                    #raise TitleAmountError
+                    raise IOError
+                for i,num in enumerate(numbers):
+                    num.strip()
+                    if num != '\n' and num != '':
+                        #attributes.append(float(num))
+                        att_dict.setdefault(att_names[i],[]).append(float(num))
+                        
+            except ValueError:
+                raise SyntaxError
+        line = fd.readline()
+        numbers = clean_line(line,delimiter)
+    
+    if line == '':
+            # end of file
+        geo_reference = None
+    else:
+        geo_reference = Geo_reference(ASCIIFile=fd,read_title=line)
+    
+    xya_dict = {}
+    xya_dict['pointlist'] = array(points).astype(Float)
+    
+    for key in att_dict.keys():
+        att_dict[key] = array(att_dict[key]).astype(Float)
+    xya_dict['attributelist'] = att_dict
+    
+    xya_dict['geo_reference'] = geo_reference
+    #print "xya_dict",xya_dict 
+    return xya_dict
+
+
+
+def _write_pts_file(file_name, point_atts):
+    """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
+
+
+    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
+    _point_atts2array(point_atts)
+    # 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 = point_atts['pointlist'].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'][:] = point_atts['pointlist'] #.astype(Float32)
+    for key in point_atts['attributelist'].keys():
+        outfile.createVariable(key, Float, ('number_of_points',))
+        outfile.variables[key][:] = point_atts['attributelist'][key] #.astype(Float32)
+        
+    if point_atts.has_key('geo_reference') and not point_atts['geo_reference'] == None:
+        point_atts['geo_reference'].write_NetCDF(outfile)
+        
+    outfile.close() 
+  
+
+
+def _write_xya_file( file_name, xya_dict, delimiter = ','):
+    """
+    export a file, ofile, with the xya format
+    
+    """
+    points = xya_dict['pointlist'] 
+    pointattributes = xya_dict['attributelist']
+    
+    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)
+    fd.close()
+    
+def _point_atts2array(point_atts):
+    point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
+    
+    for key in point_atts['attributelist'].keys():
+        point_atts['attributelist'][key]= array(point_atts['attributelist'][key]).astype(Float)
+    return point_atts
+
+   
 
 
@@ -205 +570 @@
 
     
-    
-
-
-
-
 def points_dictionary2geospatial_data(points_dictionary):
     """Convert points_dictionary to geospatial data object
@@ -229 +589 @@
                            geo_reference = geo)
 
-
-            
-            
-        
+def clean_line(line,delimiter):      
+    """Remove whitespace
+    """
+    #print ">%s" %line
+    line = line.strip()
+    #print "stripped>%s" %line
+    numbers = line.split(delimiter)
+    i = len(numbers) - 1
+    while i >= 0:
+        if numbers[i] == '':
+            numbers.pop(i)
+        else:
+            numbers[i] = numbers[i].strip()
+        
+        i += -1
+    #for num in numbers:
+    #    print "num>%s<" %num
+    return numbers
+            
+#def add_points_files(file):            
+