"""Add two point files by the union of the fine and coarse points. Course points that are in the extent of the fine points are removed. The extent of the fine points file is assumed to be a rectangle, parallel to the x and y axis. Ole Nielsen, Stephen Roberts, Duncan Gray Geoscience Australia, 2005. """ from util import outside_polygon, inside_polygon from Numeric import take, concatenate import time import exceptions class AttributeError(exceptions.Exception): pass def combine_rectangular_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']) # find extent of course points extent = extent(fine['pointlist']) # find out_points - coarse points outside of fine points extent 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 coarse = take_points(coarse, outside_coarse_indices) # 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__": """ Load in two data points files. Save a new points file. """ import os, sys usage = "usage: %s fine_input.pts coarse_input.pts output.pts " \ %os.path.basename(sys.argv[0]) if len(sys.argv) < 4: print usage else: fine_file = sys.argv[1] coarse_file = sys.argv[2] output_file = sys.argv[3] combine_rectangular_points_files(fine_file, coarse_file, output_file)