# Changeset 5419

Ignore:
Timestamp:
Jun 24, 2008, 12:37:05 PM (15 years ago)
Message:

adjustedproject file to read in order of boundary nodes from file

File:
1 edited

Unmodified
Added
Removed
• ## anuga_work/development/boxingday08/project.py

 r5404 from Scientific.IO.NetCDF import NetCDFFile from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\ compress,zeros,fabs compress,zeros,fabs,allclose from anuga.utilities.polygon import inside_polygon,read_polygon from os import sep from time import localtime, strftime, gmtime def create_sts_boundary(filename,verbose=False): fid = NetCDFFile(filename+'.sts', 'r')    #Open existing file for read def create_sts_boundary_order_file(order_filename,stsfilename,verbose=False): """ Note This is not a robust algorithm. Be Careful when applying that polygon does not intersect itself and that resulting boundary actually follows the path it is supposed too. """ fid = NetCDFFile(stsfilename+'.sts', 'r')    #Open existing file for read x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices boundary[0][1]=y[west_index] indices.remove(west_index) order=[] for i in range(n-1): min_dist=1e12 boundary[i+1][0]=x[index] boundary[i+1][1]=y[index] order.append(index) #ensure boundary is a polygon i.e no intesections #e.g. #So far i cannot tell if the wrong point is closer. d="," fid=open(order_filename,'w') header="index,longitude,latitude\n" fid.write(header) line=str(west_index)+d+str(x[west_index])+d+\ str(y[west_index])+"\n" fid.write(line) for i in order: line=str(i)+d+str(x[i])+d+\ str(y[i])+"\n" fid.write(line) fid.close() return boundary #read in polyline boundary from anuga.shallow_water.data_manager import urs2sts from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary base_name='polyline' import os if os.path.exists(base_name+'.sts'): print 'sts boundary file already created.' print 'sts boundary file already created.' else: urs2sts(base_name,mean_stage=0.0,verbose=False) boundary=create_sts_boundary(base_name,verbose=True) urs_boundary=boundary.tolist() urs2sts(base_name,mean_stage=0.0,verbose=False) order_filename='sts_order.txt' boundary=create_sts_boundary_order_file(order_filename,base_name,verbose=False) urs_boundary=create_sts_boundary(order_filename,base_name,lat_long=False) assert allclose(boundary.tolist(),urs_boundary) #Read in large domain boundary bp = read_polygon(dir+sep+extent) #Now order of boundary is correct clip sts boundary to small region for #this particular scenario indices=inside_polygon(urs_boundary, bp, closed=True, verbose=False) bounding_polygon=[] bounding_polygon.append(urs_boundary[i]) #Append clipped sts boundary to large domain boundary bounding_polygon.extend(bp[2:len(bp)])
Note: See TracChangeset for help on using the changeset viewer.