source: inundation/geospatial_data/geospatial_data.py @ 2570

"""Class Geospatial_data - Manipulation of locations on the planet and 
associated attributes.

"""


from utilities.numerical_tools import ensure_numeric

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 = None,
                 attributes = None,
                 geo_reference = None,
                 default_attribute_name = None,
                 file_name = None):
         

        """Create instance from data points and associated attributes


        data_points: x,y coordinates in meters. Type must be either a
        sequence of 2-tuples or an Mx2 Numeric array of floats.

        attributes: Associated values for each data point. The type
        must be either a list or an array of length M or a dictionary
        of lists (or arrays) of length M. In the latter case the keys
        in the dictionary represent the attribute names, in the former
        the attribute will get the default name 'attribute'.
        
        geo_reference: Object representing the origin of the data
        points. It contains UTM zone, easting and northing and data
        points are assumed to be relative to this origin.
        If geo_reference is None, the default geo ref object is used

        default_attribute_name: Name of default attribute to be used with
        get_attribute_values. The idea is that the dataset can be
        equipped with information about which attribute to return.
        If None, the default is the 'first'
        
        file_name: Name of input file.....
        
        """

        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):
        """Check and assign attributes dictionary
        """
        
        from types import DictType
       
        if attributes is None:
            self.attributes = None
            return
        
        if not isinstance(attributes, DictType):
            #Convert single attribute into dictionary
            attributes = {'attribute': attributes}

        #Check input attributes    
        for key in attributes.keys():
            try:
                attributes[key] = ensure_numeric(attributes[key])
            except:
                msg = 'Attribute %s could not be converted' %key
                msg += 'to a numeric vector'
                raise msg

        self.attributes = attributes    


    def set_geo_reference(self, geo_reference):

        from coordinate_transforms.geo_reference import Geo_reference


        if geo_reference is None:
            self.geo_reference = Geo_reference() # Use default
        elif isinstance(geo_reference, Geo_reference):
            self.geo_reference = geo_reference
        else:
            msg = 'Argument geo_reference must be a valid Geo_reference \n'
            msg += 'object or None.'
            raise msg


    def set_default_attribute_name(self, default_attribute_name):
        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):
        return self.geo_reference
       
    def get_data_points(self, absolute = True):
        if absolute is True:
            xll = self.geo_reference.xllcorner
            yll = self.geo_reference.yllcorner
            return self.data_points + [xll, yll]
            
        else: 
            return self.data_points
    
    def get_attributes(self, attribute_name = None):
        """Return values for one named attribute.

        If attribute_name is None, default_attribute_name is used
        """

        if attribute_name is None:
            if self.default_attribute_name is not None:
                attribute_name = self.default_attribute_name
            else:
                attribute_name = self.attributes.keys()[0] 
                # above line takes the first one from keys
                

        msg = 'Attribute name %s does not exist in data set' %attribute_name
        assert self.attributes.has_key(attribute_name), msg

        return self.attributes[attribute_name]

    def __add__(self,other):
        """ 
        returns the addition of 2 geospatical objects,
        objects are concatencated to the end of each other
            
        NOTE: doesn't add if objects contain different 
        attributes  
        """

        # sets xll and yll as the smallest from self and other
        if self.geo_reference.xllcorner <= other.geo_reference.xllcorner:
            xll = self.geo_reference.xllcorner
        else:
            xll = other.geo_reference.xllcorner

        if self.geo_reference.yllcorner <= other.geo_reference.yllcorner:
            yll = self.geo_reference.yllcorner
        else:
            yll = other.geo_reference.yllcorner

        # find objects zone and checks if the same
        geo_ref1 = self.get_geo_reference()
        zone1 = geo_ref1.get_zone()
        
        geo_ref2 = other.get_geo_reference()
        zone2 = geo_ref2.get_zone()
        
        geo_ref = Geo_reference(zone1, xll, yll)
        
        if zone1 == zone2:
            a_rel_points = self.get_data_points()
            b_rel_points = other.get_data_points()
            
            # subtracts xll and yll from self's and other's
            # reletive data point and concatenates
            c_points = concatenate((a_rel_points-[xll, yll],
                                    b_rel_points-[xll, yll]), axis = 0)

            new_dictionary = {}
            for x in self.attributes.keys():
                if other.attributes.has_key(x):

                    a_attrib = self.attributes[x]
                    b_attrib = other.attributes[x]
                    new_dictionary[x] = concatenate((a_attrib, b_attrib))

                else:
                    msg = 'Both geospatial_data objects must have the same \n'
                    msg += 'attributes to allow addition.'
                    raise msg
 
            return Geospatial_data(c_points, new_dictionary, geo_ref)
            
        else:
            msg = 'Both geospatial_data objects must be in the same \
            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

   


def geospatial_data2points_dictionary(geospatial_data):
    """Convert geospatial data to points_dictionary
    """

    points_dictionary = {}
    points_dictionary['pointlist'] = geospatial_data.data_points

    points_dictionary['attributelist'] = {}

    for attribute_name in geospatial_data.attributes.keys():
        val = geospatial_data.attributes[attribute_name]
        points_dictionary['attributelist'][attribute_name] = val

    points_dictionary['geo_reference'] = geospatial_data.geo_reference

    return points_dictionary

    
def points_dictionary2geospatial_data(points_dictionary):
    """Convert points_dictionary to geospatial data object
    """

    msg = 'Points dictionary must have key pointlist' 
    assert points_dictionary.has_key('pointlist'), msg

    msg = 'Points dictionary must have key attributelist'     
    assert points_dictionary.has_key('attributelist'), msg        

    if points_dictionary.has_key('geo_reference'):
        geo = points_dictionary['geo_reference']
    else:
        geo = None
    
    return Geospatial_data(points_dictionary['pointlist'],
                           points_dictionary['attributelist'],
                           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):