Changeset 5419


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

adjustedproject file to read in order of boundary nodes from file

File:
1 edited

Legend:

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

    r5404 r5419  
    22from Scientific.IO.NetCDF import NetCDFFile
    33from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
    4      compress,zeros,fabs
     4    compress,zeros,fabs,allclose
    55from anuga.utilities.polygon import inside_polygon,read_polygon
    66from os import sep
    77from time import localtime, strftime, gmtime
    8 def create_sts_boundary(filename,verbose=False):
    9     fid = NetCDFFile(filename+'.sts', 'r')    #Open existing file for read
     8
     9def create_sts_boundary_order_file(order_filename,stsfilename,verbose=False):
     10    """
     11    Note This is not a robust algorithm. Be Careful when applying that
     12    polygon does not intersect itself and that resulting boundary actually
     13    follows the path it is supposed too.
     14    """
     15    fid = NetCDFFile(stsfilename+'.sts', 'r')    #Open existing file for read
    1016    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
    1117    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
     
    3339    boundary[0][1]=y[west_index]
    3440    indices.remove(west_index)
     41    order=[]
    3542    for i in range(n-1):
    3643         min_dist=1e12
     
    4653         boundary[i+1][0]=x[index]
    4754         boundary[i+1][1]=y[index]
     55         order.append(index)
    4856    #ensure boundary is a polygon i.e no intesections
    4957    #e.g.
     
    5260             
    5361    #So far i cannot tell if the wrong point is closer.
     62    d=","
     63    fid=open(order_filename,'w')
     64    header="index,longitude,latitude\n"
     65    fid.write(header)
     66    line=str(west_index)+d+str(x[west_index])+d+\
     67        str(y[west_index])+"\n"
     68    fid.write(line)
     69    for i in order:
     70            line=str(i)+d+str(x[i])+d+\
     71                str(y[i])+"\n"
     72            fid.write(line)
     73    fid.close()
    5474    return boundary
     75           
    5576
    5677#read in polyline boundary
    57 from anuga.shallow_water.data_manager import urs2sts
     78from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
    5879base_name='polyline'
    5980import os
    6081if os.path.exists(base_name+'.sts'):
    61         print 'sts boundary file already created.'
     82        print 'sts boundary file already created.'
    6283else:
    63     urs2sts(base_name,mean_stage=0.0,verbose=False)
    64 boundary=create_sts_boundary(base_name,verbose=True)
    65 urs_boundary=boundary.tolist()
     84        urs2sts(base_name,mean_stage=0.0,verbose=False)
     85
     86order_filename='sts_order.txt'
     87boundary=create_sts_boundary_order_file(order_filename,base_name,verbose=False)
     88urs_boundary=create_sts_boundary(order_filename,base_name,lat_long=False)
     89assert allclose(boundary.tolist(),urs_boundary)
    6690
    6791#Read in large domain boundary
     
    7094bp = read_polygon(dir+sep+extent)
    7195
     96#Now order of boundary is correct clip sts boundary to small region for
     97#this particular scenario
    7298indices=inside_polygon(urs_boundary, bp, closed=True, verbose=False)
    7399bounding_polygon=[]
     
    76102    bounding_polygon.append(urs_boundary[i])
    77103
     104#Append clipped sts boundary to large domain boundary
    78105bounding_polygon.extend(bp[2:len(bp)])
    79106
Note: See TracChangeset for help on using the changeset viewer.