Changeset 2570

Timestamp:

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

Author:

nick

Message:

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

Location:

inundation/geospatial_data

Files:

• geospatial_data.py (modified) (7 diffs)
• test_geospatial_data.py (modified) (4 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):            
+        

inundation/geospatial_data/test_geospatial_data.py

@@ -6 +6 @@
 from Numeric import zeros, array, allclose, concatenate
 from math import sqrt, pi
+import tempfile
 
 from geospatial_data import *
@@ -231 +232 @@
         assert allclose(G.get_data_points(), [[2.0, 4.1], [4.0, 7.3], [5.1, 9.1], [6.1, 6.3]])                             
 
-    def xtest_create_from_file(self):
-        """Check that object can be created from a points file
+
+
+
+    def test_create_from_xya_file(self):
+        """Check that object can be created from a points file (.pts and .xya)
         """
 
@@ -245 +249 @@
                       'attributelist': {'att1': attributes,
                                         'att2': array(attributes) + 1}} 
-        
+     
         # Create points as an xya file 
         FN = 'test_points.xya'
@@ -252 +256 @@
 
         #Create object from file
-        G = Geospatial_data(filename = FN)
-
+        G = Geospatial_data(file_name = FN)
+        
         assert allclose(G.get_data_points(), points)
         assert allclose(G.get_attributes('att1'), attributes)
         assert allclose(G.get_attributes('att2'), array(attributes) + 1)
         
-
-        
         os.remove(FN)
         
+    def test_loadxya(self):
+        """
+        comma delimited
+        """
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("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()
+        #print fileName
+        dict = import_points_file(fileName,delimiter = ',')
+        os.remove(fileName)
+        assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
+        assert allclose(dict['attributelist']['speed'], [0.0, 10.0, 40.0])
+
+    def test_loadxya2(self):
+        """
+        space delimited
+        """
+        import os
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  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()
+        #print fileName
+#        dict = import_points_file(fileName,delimiter = ' ')
+        dict = import_points_file(fileName,delimiter = ' ')
+        os.remove(fileName)
+        assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
+        assert allclose(dict['attributelist']['speed'], [0.0, 10.0, 40.0])
+     
+    def test_loadxya3(self):
+        """
+        space delimited
+        """
+        import os
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  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\
+#geocrap\n\
+56\n\
+56.6\n\
+3\n")
+        file.close()
+        #print fileName
+        dict = import_points_file(fileName,delimiter = ' ')
+        os.remove(fileName)
+        assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
+        assert allclose(dict['attributelist']['speed'], [0.0, 10.0, 40.0])
+
+    def test_loadxy_bad_no_file_xya(self):
+        import os
+       
+        fileName = tempfile.mktemp(".xya")
+        #print fileName
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        '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 Xtest_read_write_points_file_bad2(self):
+        dict = {}
+        att_dict = {}
+        dict['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)
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad points file extension did not raise error!')
+                   
+    def test_loadxy_bad(self):
+        import os
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  elevation   \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()
+        #print fileName
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad xya file did not raise error!')
+        os.remove(fileName)
+       
+    def test_loadxy_bad2(self):
+        import os
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("elevation\n\
+1.0 0.0 10.0 \n\
+0.0 1.0\n\
+1.0 \n")
+        file.close()
+        #print fileName
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad xya file did not raise error!')    
+        os.remove(fileName)
+   
+    def test_loadxy_bad3(self):
+        """
+        specifying wrong delimiter
+        """
+        import os
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  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()
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad xya file did not raise error!')    
+        os.remove(fileName)
+     
+    def test_loadxy_bad4(self):
+        """
+        specifying wrong delimiter
+        """
+        import os
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  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\
+yeah")
+        file.close()
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad xya file did not raise error!')    
+        os.remove(fileName)
+  
+    def test_loadxy_bad4(self):
+        """
+        specifying wrong delimiter
+        """
+        import os
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  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\
+#geocrap")
+        file.close()
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad xya file did not raise error!')    
+        os.remove(fileName)
+
+    def test_loadxy_bad5(self):
+        """
+        specifying wrong delimiter
+        """
+        import os
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  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\
+#geocrap\n\
+crap")
+        file.close()
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad xya file did not raise error!')    
+        os.remove(fileName)              
+
+    def test_loadxy_bad_no_file_xya(self):
+        import os
+       
+        fileName = tempfile.mktemp(".xya")
+        #print fileName
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'imaginary file did not raise error!')
+                        
+  ###################### .XYA ##############################
+        
+    def test_export_xya_file(self):
+        dict = {}
+        att_dict = {}
+        dict['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)
+        dict2 = import_points_file(fileName)
+        #print "fileName",fileName 
+        os.remove(fileName)
+        #print "dict2",dict2 
+        
+        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])
+        answer = [10.0, 0.0, 10.4]
+        assert allclose(dict2['attributelist']['brightness'], answer)
+        #print "dict2['geo_reference']",dict2['geo_reference'] 
+        self.failUnless(dict['geo_reference'] == dict2['geo_reference'],
+                         'test_writepts failed. Test geo_reference')
+
+    def test_export_xya_file2(self):
+        dict = {}
+        att_dict = {}
+        dict['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)
+        dict2 = import_points_file(fileName)
+        #print "fileName",fileName 
+        os.remove(fileName)
+        #print "dict2",dict2 
+        
+        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])
+        answer = [10.0, 0.0, 10.4]
+        assert allclose(dict2['attributelist']['brightness'], answer)
+
+
+    def test_loadpts(self):
+        
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        fileName = tempfile.mktemp(".pts")
+        # NetCDF file definition
+        outfile = NetCDFFile(fileName, 'w')
+        
+        # dimension definitions
+        outfile.createDimension('number_of_points', 3)    
+        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
+    
+        # variable definitions
+        outfile.createVariable('points', Float, ('number_of_points',
+                                                 'number_of_dimensions'))
+        outfile.createVariable('elevation', Float, ('number_of_points',))
+    
+        # Get handles to the variables
+        points = outfile.variables['points']
+        elevation = outfile.variables['elevation']
+ 
+        points[0, :] = [1.0,0.0]
+        elevation[0] = 10.0 
+        points[1, :] = [0.0,1.0]
+        elevation[1] = 0.0  
+        points[2, :] = [1.0,0.0]
+        elevation[2] = 10.4    
+
+        outfile.close()
+        
+        dict = import_points_file(fileName)
+        os.remove(fileName)
+        answer =  [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
+        #print "dict['pointlist']",dict['pointlist']
+        #print "answer",answer
+        assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
+        
+    def test_writepts(self):
+        dict = {}
+        att_dict = {}
+        dict['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(".pts")
+        export_points_file(fileName, dict)
+        dict2 = import_points_file(fileName)
+        #print "fileName",fileName 
+        os.remove(fileName)
+        #print "dict2",dict2 
+        
+        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])
+        answer = [10.0, 0.0, 10.4]
+        assert allclose(dict2['attributelist']['brightness'], answer)
+
+        #print "dict['geo_reference'] ",dict['geo_reference'] 
+        #print "dict2['geo_reference']",dict2['geo_reference']
+        
+        self.failUnless(dict['geo_reference'] == dict2['geo_reference'],
+                         'test_writepts failed. Test geo_reference')
+        
+ ########################## BAD .PTS ##########################          
+
+    def test_load_bad_no_file_pts(self):
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".pts")
+        #print fileName
+        try:
+            dict = import_points_file(fileName)
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'imaginary file did not raise error!')
+
+
+    def test_create_from_pts_file(self):
+        
+        from Scientific.IO.NetCDF import NetCDFFile
+
+
+#        fileName = tempfile.mktemp(".pts")
+        FN = 'test_points.pts'
+        # NetCDF file definition
+        outfile = NetCDFFile(FN, 'w')
+        
+        # dimension definitions
+        outfile.createDimension('number_of_points', 3)    
+        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
+    
+        # variable definitions
+        outfile.createVariable('points', Float, ('number_of_points',
+                                                 'number_of_dimensions'))
+        outfile.createVariable('elevation', Float, ('number_of_points',))
+    
+        # Get handles to the variables
+        points = outfile.variables['points']
+        elevation = outfile.variables['elevation']
+ 
+        points[0, :] = [1.0,0.0]
+        elevation[0] = 10.0 
+        points[1, :] = [0.0,1.0]
+        elevation[1] = 0.0  
+        points[2, :] = [1.0,0.0]
+        elevation[2] = 10.4    
+
+        outfile.close()
+
+        G = Geospatial_data(file_name = FN)
+
+        assert allclose(G.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
+        os.remove(FN)
 
 if __name__ == "__main__":