Changeset 2584


Ignore:
Timestamp:
Mar 23, 2006, 5:07:33 PM (18 years ago)
Author:
nick
Message:

add_points_files

Location:
inundation/geospatial_data
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • inundation/geospatial_data/geospatial_data.py

    r2570 r2584  
    158158        if file_name is None:
    159159            self.file_name = None
    160 #        else:
    161 #            if access(self.file_name, F_OK) == 0 :
    162 #                msg = 'Geospatial_data object needs have either data_points \n'
    163 #                msg += ' or a file with data points'
    164 #                raise msg
     160
    165161        else:   
    166162            self.file_name = file_name
    167163            self = self.import_points_file(file_name)
    168 #            print 'data points from set self file name', self.data_points
    169 #            print 'attributes from set self file name', self.attributes
    170 #            print 'geo_ref from set self file name', self.geo_reference
    171164        return self
    172165           
     
    230223        zone2 = geo_ref2.get_zone()
    231224       
    232         geo_ref = Geo_reference(zone1, xll, yll)
    233        
    234225        if zone1 == zone2:
    235226            a_rel_points = self.get_data_points()
     
    241232                                    b_rel_points-[xll, yll]), axis = 0)
    242233
    243             new_dictionary = {}
     234            new_attributes = {}
    244235            for x in self.attributes.keys():
    245236                if other.attributes.has_key(x):
     
    247238                    a_attrib = self.attributes[x]
    248239                    b_attrib = other.attributes[x]
    249                     new_dictionary[x] = concatenate((a_attrib, b_attrib))
     240                    new_attributes[x] = concatenate((a_attrib, b_attrib))
    250241
    251242                else:
     
    253244                    msg += 'attributes to allow addition.'
    254245                    raise msg
    255  
    256             return Geospatial_data(c_points, new_dictionary, geo_ref)
     246               
     247            geo_ref = Geo_reference(zone1, xll, yll)
     248            return Geospatial_data(c_points, new_attributes, geo_ref)
    257249           
    258250        else:
     
    263255     
    264256   
    265     ###
    266     #  IMPORT/EXPORT POINTS FILES
    267     ###
    268     '''
    269     def export_points_file(ofile, point_dict):
    270         """
    271         write a points file, ofile, as a binary (.xya) or text (.pts) file
    272 
    273         ofile is the file name, including the extension
    274 
    275         The point_dict is defined at the top of this file.
    276         """
    277         #this was done for all keys in the mesh file.
    278         #if not mesh_dict.has_key('points'):
    279         #    mesh_dict['points'] = []
    280         if (ofile[-4:] == ".xya"):
    281             _write_xya_file(ofile, point_dict)
    282         elif (ofile[-4:] == ".pts"):
    283             _write_pts_file(ofile, point_dict)
    284         else:
    285             msg = 'Unknown file type %s ' %ofile
    286         raise IOError, msg
    287 
    288     '''       
     257###
     258#  IMPORT/EXPORT POINTS FILES
     259###
     260
    289261def import_points_file( ofile, delimiter = None, verbose = False):
    290262    """ load an .xya or .pts file
     
    312284            msg = 'Could not open file %s ' %ofile
    313285            raise IOError, msg
    314         except IOError:
     286        except IOError, e:
    315287            # Catch this to add an error message
    316             msg = 'Could not open file %s ' %ofile
     288            msg = 'Could not open file or incorrect file format %s:%s' %(ofile, e)
    317289            raise IOError, msg
    318290           
     
    416388    line = fd.readline()
    417389    numbers = clean_line(line,delimiter)
     390    print 'read xya numbers', numbers
    418391    while len(numbers) > 1:
    419392        if numbers != []:
     
    608581    return numbers
    609582           
    610 #def add_points_files(file):           
    611        
     583def add_points_files(add_file1, add_file2, results_file):
     584    ''' adds the points and attruibutes of 2 pts or xya files and
     585    writes it to a pts file
     586   
     587    NOTE will add the function to check and remove points from one set
     588    that are shared. This will require some work and maybe __subtract__ function
     589    '''
     590#    print 'load %s' %add_file1
     591#    if access(add_file1,F_OK) == 1 :
     592#        print'file1', add_file1
     593
     594#    print'file1 again', add_file1
     595#    print'file2 again', add_file2
     596    print 'create G1'
     597    G1 = Geospatial_data(file_name = add_file1)
     598   
     599#    print 'G1 dict'
     600 #   G1_points_dict = geospatial_data2points_dictionary(G1)
     601
     602#    print 'G1 xya file'
     603#    export_points_file("v:\inundation\onslow_tsunami_scenario_2006\\topographies\onshore.xya", G1_points_dict)
     604
     605    print 'create G2'
     606#    print 'load %s' %add_file2
     607#    if access(add_file2,F_OK) == 1 :
     608#        print'file2', add_file2
     609    G2 = Geospatial_data(file_name = add_file2)
     610    print 'G2 dict'
     611#    G2_points_dict = geospatial_data2points_dictionary(G2)
     612
     613    new_add_file2 = add_file2[:-4] + '.pts'
     614    print 'G2 pts file', new_add_file2
     615   
     616#    export_points_file(new_add_file2, G2_points_dict)
     617    print 'G1+g2'   
     618    G = G1 + G2
     619    #FIXME remove dependance on points to dict in export only!
     620    G_points_dict = geospatial_data2points_dictionary(G)
     621   
     622    export_points_file(results_file, G_points_dict)
     623   
     624   
     625     
     626       
  • inundation/geospatial_data/test_geospatial_data.py

    r2570 r2584  
    190190        G2 = Geospatial_data(points, attributes1)
    191191       
    192 #        print 'G1 depth=', G1.get_attributes('depth')
    193 #        print 'G1 attrib keys=', G1.attributes.keys()
    194         g3 = geospatial_data2points_dictionary(G1)
     192#        g3 = geospatial_data2points_dictionary(G1)
    195193#        print 'g3=', g3
    196194       
    197195        G = G1 + G2
    198196
    199 #        g = geospatial_data2points_dictionary(G)
    200 #        print 'G=', g
    201 #        print 'G points =', G.get_data_points()
    202 #        print 'G attrib keys =', G.attributes.keys()
    203 #        print 'G test =', G.get_attributes('elevation')
    204197        assert G.attributes.has_key('depth')
    205198        assert G.attributes.has_key('elevation')
     
    239232        """
    240233
    241         from load_mesh.loadASCII import export_points_file
    242        
    243 
    244234        points = [[1.0, 2.1], [3.0, 5.3], [5.0, 6.1], [6.0, 3.3]]
    245235        attributes = [2, 4, 5, 76]
     
    352342                        'bad points file extension did not raise error!')
    353343
    354     def Xtest_read_write_points_file_bad2(self):
     344    def test_read_write_points_file_bad2(self):
    355345        dict = {}
    356346        att_dict = {}
     
    682672        assert allclose(G.get_attributes(), [10.0, 0.0, 10.4])
    683673        os.remove(FN)
     674       
     675    def test_add_points_files(self):
     676        '''adds an xya and pts files and checks resulting file is correct
     677        '''
     678       
     679        # create files
     680        dict1 = {}
     681        att_dict1 = {}
     682        dict1['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
     683        att_dict1['elevation'] = array([-10.0, 0.0, 10.4])
     684        att_dict1['brightness'] = array([10.0, 0.0, 10.4])
     685        dict1['attributelist'] = att_dict1
     686        dict1['geo_reference'] = Geo_reference(56,1.9,1.9)
     687       
     688        dict2 = {}
     689        att_dict2 = {}
     690        dict2['pointlist'] = array([[2.0, 1.0],[1.0, 2.0],[2.0, 1.0]])
     691        att_dict2['elevation'] = array([1.0, 15.0, 1.4])
     692        att_dict2['brightness'] = array([14.0, 1.0, -12.4])
     693        dict2['attributelist'] = att_dict2
     694        dict2['geo_reference'] = Geo_reference(56,1.9,1.9)
     695       
     696        fileName1 = tempfile.mktemp(".xya")
     697        fileName2 = tempfile.mktemp(".pts")
     698
     699        export_points_file(fileName1, dict1)
     700        export_points_file(fileName2, dict2)
     701#        export_points_file('hello.xya', dict1)
     702#        export_points_file('hello.pts', dict2)
     703       
     704        results_file = 'resulting_points.pts'
     705       
     706        # add files
     707        add_points_files(fileName1, fileName2, results_file )
     708#        add_points_files('hello.xya', 'hello.pts', results_file )
     709       
     710        #read results
     711        results_dict = import_points_file(results_file)
     712       
     713        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]])
     714        assert allclose(results_dict['attributelist']['elevation'], [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4])
     715        answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4]
     716        assert allclose(results_dict['attributelist']['brightness'], answer)
     717        #print "dict2['geo_reference']",dict2['geo_reference']
     718        self.failUnless(results_dict['geo_reference'] == dict2['geo_reference'],
     719                         'test_writepts failed. Test geo_reference')
     720                         
     721        os.remove(fileName1)
     722        os.remove(fileName2)
     723        os.remove(results_file)
    684724
    685725if __name__ == "__main__":
Note: See TracChangeset for help on using the changeset viewer.