source: inundation/pyvolution/combine_pts.py @ 1740

Last change on this file since 1740 was 1423, checked in by duncan, 19 years ago

I've removed the function mesh_file_to_mesh_dictionary. Please use the function import_mesh_file instead.

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