Changeset 1101

Timestamp:

Mar 17, 2005, 2:38:40 PM (20 years ago)

Author:

duncan

Message:

bug fix

Location:

inundation/ga/storm_surge/pyvolution

Files:

• combine_pts.py (modified) (7 diffs)
• test_combine_pts.py (modified) (3 diffs)

inundation/ga/storm_surge/pyvolution/combine_pts.py

@@ -8 +8 @@
    Geoscience Australia, 2005.   
 """
-from util import outside_polygon
+from util import outside_polygon, inside_polygon
 from Numeric import take, concatenate
+import time
 
 import exceptions
@@ -16 +17 @@
 def combine_rectangular_points_files(fine_points_file,
                                     coarse_points_file,
-                                    output_points_file):
+                                    output_points_file,
+                                     verbose = False):
     
     from load_mesh.loadASCII import load_points_file, extent, point_atts2array, write_pts, add_point_dictionaries,take_points
@@ -22 +24 @@
     
     # load fine points file
+    if verbose:print "loading Point files"
     try:
         fine = load_points_file(fine_points_file,delimiter = ',')
@@ -35 +38 @@
         coarse = load_points_file(coarse_points_file,
                                               delimiter = ' ')
+        
+    if verbose:print "doing other stuff"
     # convert to Numeric    
     fine = point_atts2array(fine)
@@ -43 +48 @@
     if not fine['attributelist'].keys() == coarse['attributelist'].keys():
         raise AttributeError
+    
     # set points to the same geo ref -  say the fine points geo ref
     if fine.has_key('geo_reference')and not fine['geo_reference'] is None:
@@ -54 +60 @@
     # find out_points - coarse points outside of fine points extent
 
-    #FIXME the closed doesn't seem to work here,
-    #though it works in the test case
+    if verbose:
+        print "starting outside_polygon"
+        t0 = time.time()
     outside_coarse_indices = outside_polygon(coarse['pointlist'],
                                              extent,closed=True)
+    if verbose:
+        print "Points outside determined"
+        print 'That took %.2f seconds' %(time.time()-t0)
 
     # Remove points from the coarse data
@@ -63 +73 @@
         
     # add fine points and out_points
+    if verbose:print "Adding points"
     combined = add_point_dictionaries(fine, coarse)
 
     # save
+    if verbose:print "writing points"
     write_pts(output_points_file, combined)
 
+def add_points_files(fine_points_file,
+                                    coarse_points_file,
+                                    output_points_file,
+                                     verbose = False):
+    
+    from load_mesh.loadASCII import load_points_file, extent, point_atts2array, write_pts, add_point_dictionaries,take_points
+
+    
+    # load fine points file
+    if verbose:print "loading Point files"
+    try:
+        fine = load_points_file(fine_points_file,delimiter = ',')
+    except SyntaxError,e:
+        fine = load_points_file(fine_points_file,delimiter = ' ')
+    
+    
+    # load course points file
+    try:
+        coarse = load_points_file(coarse_points_file,
+                                              delimiter = ',')
+    except SyntaxError,e:
+        coarse = load_points_file(coarse_points_file,
+                                              delimiter = ' ')
+        
+    if verbose:print "doing other stuff"
+    # convert to Numeric    
+    fine = point_atts2array(fine)
+    coarse = point_atts2array(coarse)
+    
+    #check that the attribute info is the same
+    #FIXME is not sorting this a problem?
+    if not fine['attributelist'].keys() == coarse['attributelist'].keys():
+        raise AttributeError
+    
+    # set points to the same geo ref -  say the fine points geo ref
+    if fine.has_key('geo_reference')and not fine['geo_reference'] is None:
+        if coarse.has_key('geo_reference'):
+            coarse['pointlist'] = \
+             fine['geo_reference'].change_points_geo_ref(coarse['pointlist'],
+                                                         points_geo_ref=coarse['geo_reference'])  
+        
+    # add fine points and out_points
+    if verbose:print "Adding points"
+    combined = add_point_dictionaries(fine, coarse)
+
+    # save
+    if verbose:print "writing points"
+    write_pts(output_points_file, combined)
+
+def reduce_points_to_mesh_extent(points_file,
+                                    mesh_file,
+                                    output_points_file,
+                                     verbose = False):
+    
+    from load_mesh.loadASCII import load_points_file, extent, point_atts2array, write_pts, add_point_dictionaries,take_points, mesh_file_to_mesh_dictionary
+
+    
+    # load  points file
+    try:
+        points = load_points_file(points_file,delimiter = ',')
+    except SyntaxError,e:
+        points = load_points_file(points_file,delimiter = ' ')
+    
+    # load  mesh file
+    mesh = mesh_file_to_mesh_dictionary(mesh_file)
+
+    # get extent of mesh
+    extent = extent(mesh['vertices'])
+    #print "extent",extent 
+    # set points to the same geo ref - the points geo ref
+    if points.has_key('geo_reference')and not points['geo_reference'] is None:
+        if mesh.has_key('geo_reference'):
+            extent = \
+             points['geo_reference'].change_points_geo_ref(extent,
+                                                         points_geo_ref=mesh['geo_reference'])  
+    #print "extent",extent 
+
+    #get a list of in point indexes
+    #FIXME closed doesn't seems to work here.
+    inside_indices = inside_polygon(points['pointlist'],
+                                             extent,closed=True)
+    #create a list of in points
+    points = take_points(points, inside_indices)
+
+    # save the list, along with geo-ref
+    write_pts(output_points_file,points)
+    
     #-------------------------------------------------------------
 if __name__ == "__main__":

inundation/ga/storm_surge/pyvolution/test_combine_pts.py

@@ -27 +27 @@
         import tempfile
         import os
+
 
         # create a fine .pts file
@@ -63 +64 @@
         results = load_points_file(out_file,
                                   delimiter = ',')
+        answer = [[2.0, 0.0],
+                  [0.0, 1.0],
+                  [0.0, 4.0],
+                  [4.0, 4.0],
+                  [4.0, 1.0]]
         #print "results",results 
-        assert allclose(results['pointlist'], [[4.0, 3.0],
-                                            [2.0, 0.0],
-                                            [0.0, 1.0],
-                                            [0.0, 4.0],
-                                            [4.0, 4.0],
-                                            [4.0, 1.0]])
-        assert allclose(results['attributelist']['sonic'], [ 3., -2.,
+        #print "answer",answer
+        
+        self.failUnless(len(results['pointlist']) == len(answer),
+                         'final number of points wrong. failed.')
+        
+        assert allclose(results['pointlist'], answer)
+        assert allclose(results['attributelist']['sonic'], [ -2.,
                                                              1.,  4.,
                                                              8.,  5.])
-        assert allclose(results['attributelist']['elevation'],[30., -20.,
+        assert allclose(results['attributelist']['elevation'],[ -20.,
                                                                10.,  40.,
                                                                80.,  50.])
+        
+        self.failUnless(results['geo_reference'] == fine_dict['geo_reference'],
+                         ' failed.')
+        #clean up
+        os.remove(out_file)
+
+    def test_combine_rectangular_points_filesII(self):
+        from load_mesh.loadASCII import write_pts
+        import tempfile
+        import os
+
+        # create a fine .pts file
+        fine_dict = {}
+        fine_dict['pointlist']=[[0,1],[0,4],[4,4],[4,1],[3,1],[2,2],[1,3],[3,4]]
+        att_dict = {}
+        fine_dict['attributelist'] = att_dict
+        fine_dict['geo_reference'] = Geo_reference(56,2.0,1.0)
+        
+        fine_file = tempfile.mktemp(".pts")
+        write_pts(fine_file, fine_dict)
+
+
+        # create a coarse .pts file
+        coarse_dict = {}
+        coarse_dict['pointlist']=[[0,1],[0,0],[0.5,0.5],[1,1],
+                                  [1.5,1.5],[2,1],[0,-2],[100,10],
+                                  [-20,4],[-50,5],[60,70]]
+        att_dict = {}
+        coarse_dict['attributelist'] = att_dict
+        coarse_dict['geo_reference'] = Geo_reference(56,4.0,3.0)
+        
+        coarse_file = tempfile.mktemp(".pts")
+        write_pts(coarse_file, coarse_dict)
+        
+        out_file = tempfile.mktemp(".pts")
+
+        combine_rectangular_points_files(fine_file,coarse_file,out_file)
+
+        #clean up
+        #os.remove(fine_file)
+        #os.remove(coarse_file)
+
+        results = load_points_file(out_file,
+                                  delimiter = ',')
+        answer = [[2.0, 0.0],
+                  [102.,12.],
+                  [-18.,6.],
+                  [-48.,7.],
+                  [62.,72.],
+                  [0.0, 1.0],
+                  [0.0, 4.0],
+                  [4.0, 4.0],
+                  [4.0, 1.0],
+                  [3.,1.],
+                  [2.,2.],
+                  [1.,3.],
+                  [3.,4.]]
+        #print "results",results['pointlist']
+        #print "answer",answer
+        #print "len(results['pointlist']",len(results['pointlist'])
+        #print "len(answer)",len(answer)
+        
+        self.failUnless(len(results['pointlist']) == len(answer),
+                         'final number of points wrong. failed.')
+        assert allclose(results['pointlist'], answer)
         
         self.failUnless(results['geo_reference'] == fine_dict['geo_reference'],
@@ -123 +194 @@
         os.remove(fine_file)
         os.remove(coarse_file)
+      
+    def test_reduce_points_to_mesh_extent(self):
+        from load_mesh.loadASCII import write_pts, export_mesh_file
+        import tempfile
+        import os
+        x_origin = -45435345.
+        y_origin = 433432432.
+        # create a fine .pts file
+        fine_dict = {}
+        fine_dict['pointlist']=[[1.,1.],
+                                [1.,7.],
+                                [7.,1.],
+                                [11.,11.],
+                                [7.,8.],
+                                [8.,8.],
+                                [9.,8.]]
+        att_dict = {}
+        att_dict['elevation'] = [10,40,80,50,78,78,45]
+        att_dict['sonic'] = [1,4,8,5,56,34,213]
+        fine_dict['attributelist'] = att_dict
+        fine_dict['geo_reference'] = Geo_reference(56,x_origin,y_origin)
+        
+        points_file = tempfile.mktemp(".pts")
+        write_pts(points_file, fine_dict)
+
+
+        # create a coarse .pts file
+        mesh = {}
+        mesh['vertices']=[[0,0],
+                          [0,3],
+                          [3,3],
+                          [3,0],
+                          [1,2]
+                          ]
+        mesh['vertex_attributes']=[[],
+                          [],
+                          [],
+                          [],
+                          []
+                          ]
+        mesh['geo_reference'] = Geo_reference(56,x_origin+7.,y_origin+7.)
+        
+        mesh_file = tempfile.mktemp(".tsh")
+        export_mesh_file(mesh_file, mesh)
+        
+        out_file = tempfile.mktemp(".pts")
+
+        reduce_points_to_mesh_extent(points_file,mesh_file,out_file)
+
+        #clean up
+        #os.remove(fine_file)
+        #os.remove(coarse_file)
+
+        results = load_points_file(out_file,
+                                  delimiter = ',')
+        answer = [
+            [7.,8.], 
+            [8.,8.],
+            [9.,8.]]
+        #print "results",results['pointlist']
+        #print "answer",answer
+        
+        self.failUnless(len(results['pointlist']) == len(answer),
+                         'final number of points wrong. failed.')
+        assert allclose(results['pointlist'],answer)
+
+        answer = [78., 78.,45.]
+        
+        self.failUnless(len(results['attributelist']['elevation']) == len(answer),
+                         'final number of points wrong. failed.')
+        assert allclose(results['attributelist']['elevation'], answer)
+        
+        #clean up
+        os.remove(out_file)
+
+        #FIXME do test for add points files
         
 #-------------------------------------------------------------
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_combine_pts,'test')
-
-    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_msh_netcdf_fileII')
-    #suite = unittest.makeSuite(Test_Least_Squares,'test_fit_to_mesh_fileII')
+    #suite = unittest.makeSuite(Test_combine_pts,'test_reduce_points_to_mesh_extent')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)