source: anuga_core/source/obsolete_code/combine_pts.py @ 7585

Last change on this file since 7585 was 4235, checked in by duncan, 17 years ago

obsolete code

File size: 6.5 KB
Line 
1"""Functionality for aggregation of points files
2
3   Ole Nielsen, Stephen Roberts, Duncan Gray
4   Geoscience Australia, 2005.   
5"""
6
7from anuga.utilities.polygon import outside_polygon, inside_polygon
8from Numeric import take, concatenate
9import time
10
11import exceptions
12class AttributeError(exceptions.Exception): pass
13
14def combine_rectangular_points_files(fine_points_file,
15                                     coarse_points_file,
16                                     output_points_file,
17                                     verbose = False):
18    """Add two point files by the union of the fine and coarse points.
19   
20       Course points that are in the extent of the fine points are removed.
21
22       The extent of the fine points file is assumed to be a rectangle, parallel
23       to the x and y axis.
24    """
25   
26    from load_mesh.loadASCII import import_points_file, extent, point_atts2array, export_points_file, add_point_dictionaries, take_points
27
28   
29    # load fine points file
30    if verbose: print "loading Point files"
31    fine = import_points_file(fine_points_file)
32   
33    # load course points file
34    coarse = import_points_file(coarse_points_file)
35       
36    if verbose: print "doing other stuff"
37    # convert to Numeric   
38    fine = point_atts2array(fine)
39    coarse = point_atts2array(coarse)
40   
41    #check that the attribute info is the same
42    #FIXME is not sorting this a problem?
43    if not fine['attributelist'].keys() == coarse['attributelist'].keys():
44        raise AttributeError
45   
46    # set points to the same geo ref -  say the fine points geo ref
47    if fine.has_key('geo_reference')and not fine['geo_reference'] is None:
48        if coarse.has_key('geo_reference'):
49            coarse['pointlist'] = \
50                fine['geo_reference'].change_points_geo_ref(coarse['pointlist'],
51                                                            points_geo_ref=coarse['geo_reference']) 
52    # find extent of course points
53    extent = extent(fine['pointlist'])   
54
55    # find out_points - coarse points outside of fine points extent
56
57    if verbose:
58        print "starting outside_polygon"
59        t0 = time.time()
60    outside_coarse_indices = outside_polygon(coarse['pointlist'],
61                                             extent, closed=True)
62    if verbose:
63        print "Points outside determined"
64        print 'That took %.2f seconds' %(time.time()-t0)
65
66    # Remove points from the coarse data
67    coarse = take_points(coarse, outside_coarse_indices)
68       
69    # add fine points and out_points
70    if verbose: print "Adding points"
71    combined = add_point_dictionaries(fine, coarse)
72
73    # save
74    if verbose: print "writing points"
75    export_points_file(output_points_file, combined)
76
77def add_points_files(fine_points_file,
78                     coarse_points_file,
79                     output_points_file,
80                     verbose = False):
81    """
82    """
83   
84    from load_mesh.loadASCII import import_points_file, extent, point_atts2array, export_points_file, add_point_dictionaries,take_points
85
86   
87    # load fine points file
88    if verbose:print "loading Point files"
89    try:
90        fine = import_points_file(fine_points_file,delimiter = ',')
91    except SyntaxError,e:
92        fine = import_points_file(fine_points_file,delimiter = ' ')
93   
94   
95    # load course points file
96    try:
97        coarse = import_points_file(coarse_points_file,
98                                              delimiter = ',')
99    except SyntaxError,e:
100        coarse = import_points_file(coarse_points_file,
101                                              delimiter = ' ')
102       
103    if verbose:print "doing other stuff"
104    # convert to Numeric   
105    fine = point_atts2array(fine)
106    coarse = point_atts2array(coarse)
107   
108    #check that the attribute info is the same
109    #FIXME is not sorting this a problem?
110    if not fine['attributelist'].keys() == coarse['attributelist'].keys():
111        raise AttributeError
112   
113    # set points to the same geo ref -  say the fine points geo ref
114    if fine.has_key('geo_reference')and not fine['geo_reference'] is None:
115        if coarse.has_key('geo_reference'):
116            coarse['pointlist'] = \
117             fine['geo_reference'].change_points_geo_ref(coarse['pointlist'],
118                                                         points_geo_ref=coarse['geo_reference']) 
119       
120    # add fine points and out_points
121    if verbose: print "Adding points"
122    combined = add_point_dictionaries(fine, coarse)
123
124    # save
125    if verbose: print "writing points"
126    export_points_file(output_points_file, combined)
127
128def reduce_points_to_mesh_extent(points_file,
129                                 mesh_file,
130                                 output_points_file,
131                                 verbose = False):
132    """
133    """   
134   
135    from load_mesh.loadASCII import import_points_file, extent, point_atts2array, export_points_file, add_point_dictionaries,take_points, import_mesh_file
136
137   
138    # load  points file
139    try:
140        points = import_points_file(points_file,delimiter = ',')
141    except SyntaxError,e:
142        points = import_points_file(points_file,delimiter = ' ')
143   
144    # load  mesh file
145    mesh = import_mesh_file(mesh_file)
146
147    # get extent of mesh
148    extent = extent(mesh['vertices'])
149    #print "extent",extent
150    # set points to the same geo ref - the points geo ref
151    if points.has_key('geo_reference')and not points['geo_reference'] is None:
152        if mesh.has_key('geo_reference'):
153            extent = \
154             points['geo_reference'].change_points_geo_ref(extent,
155                                                           points_geo_ref=mesh['geo_reference']) 
156    #print "extent",extent
157
158    inside_indices = inside_polygon(points['pointlist'],
159                                    extent, closed=True)
160    #create a list of in points
161    points = take_points(points, inside_indices)
162
163    # save the list, along with geo-ref
164    export_points_file(output_points_file,points)
165   
166    #-------------------------------------------------------------
167if __name__ == "__main__":
168    """
169    Load in two data points files.
170    Save a new points file.
171    """
172    import os, sys
173    usage = "usage: %s fine_input.pts coarse_input.pts output.pts " \
174            %os.path.basename(sys.argv[0])
175
176    if len(sys.argv) < 4:
177        print usage
178    else:
179        fine_file  = sys.argv[1]
180        coarse_file = sys.argv[2]
181        output_file = sys.argv[3]
182       
183        combine_rectangular_points_files(fine_file,
184                                         coarse_file,
185                                         output_file)
Note: See TracBrowser for help on using the repository browser.