source: anuga_work/development/boxingday08/project.py @ 5402

from anuga.utilities.numerical_tools import ensure_numeric
from Scientific.IO.NetCDF import NetCDFFile
from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
     compress,zeros,fabs
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
    x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
    y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
    
    #find indices of the southernmost gauges
    #if more than one find gauges furthest left
    #will fail if domain crosses 180 degrees
    #i.e. if left most most point is not west i.e x=1 and other points 170<x<180
    #or even in utm just a diffenet zone
    south=min(y)
    south_indices=compress(y-south<1e-6,arange(len(y)))
    west=x[south_indices[0]]
    west_index=south_indices[0]
    if south_indices>1:
         #find most western gauge
         for i,south_index in enumerate(south_indices[1:]):
              if x[south_index]<west:
                   west_index=south_index
    
    n=len(x)
    boundary=zeros((n,2),Float)
    indices=arange(n).tolist()
    boundary[0][0]=x[west_index]
    boundary[0][1]=y[west_index]
    indices.remove(west_index)
    for i in range(n-1):
         min_dist=1e12
         for j,k in enumerate(indices):
              dist=sqrt((boundary[i][0]-x[k])**2+(boundary[i][1]-y[k])**2)
              if dist<min_dist:
                   min_dist=dist
                   index=k
              elif fabs(dist-min_dist)<1e-6:
                   msg='Two or more points are equal distance apart. Please resolve'
                   raise msg
         indices.remove(index)
         boundary[i+1][0]=x[index]
         boundary[i+1][1]=y[index]
    #ensure boundary is a polygon i.e no intesections
    #e.g.
    #msg='Boundary is not a regular polygon. Intersections found'
    #assert boundary.is_polygon, msg
              
    #So far i cannot tell if the wrong point is closer.
    return boundary

#read in polyline boundary
from anuga.shallow_water.data_manager import urs2sts
base_name='polyline'
import os
if os.path.exists(base_name+'.sts'):
        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()

#Read in large domain boundary
dir='mesh_polygons'
extent = 'extent.csv'
bp = read_polygon(dir+sep+extent)

indices=inside_polygon(urs_boundary, bp, closed=True, verbose=False)
bounding_polygon=[]
for i in indices[:-2]:
    #disregard last two entries because they are not on boundary
    bounding_polygon.append(urs_boundary[i])

bounding_polygon.extend(bp[2:len(bp)])

###############################
# Domain definitions
###############################

# Bathymetry and topography filenames
dir='data'

time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
poor_combined_dir_name = dir+sep+'poor_phuket_bathymetry'
good_combined_dir_name = dir+sep+'good_phuket_bathymetry'
meshname = 'phuket_mesh'
mesh_elevname = 'phuket_mesh_elev'
boundary_most_in = dir+sep+'out'
boundary_most_out = dir+sep+'most_boundary_condition'

dir='mesh_polygons'
extent = 'extent.csv'
patong_bay_polygon = read_polygon(dir+sep+'patong_bay_polygon.csv')

###############################
# Interior region definitions
###############################

contour20m = read_polygon(dir+sep+'20m.csv')
contour20m[0][1]+=3000
contour20m[-2][0]+=3000
contour20m[-2][1]-=13000
contour20m[-1][1]+=9000
patong = read_polygon(dir+sep+'patong.csv')
bay = read_polygon(dir+sep+'bay.csv')

plot=False
if plot:
    from pylab import plot,show,hold
    bounding_polygon=ensure_numeric(bounding_polygon)#for plotting only
    plot(bounding_polygon[:,0],bounding_polygon[:,1],'o')
    plot(bounding_polygon[:,0],bounding_polygon[:,1])
    contour20m = ensure_numeric(contour20m)
    plot(contour20m[:,0],contour20m[:,1])
    patong=ensure_numeric(patong)
    plot(patong[:,0],patong[:,1])
    show()

east  = 97.7
west  = 96.7
north = 9.0
south = 5.9

check_ferret_extent=False
if check_ferret_extent:
    east  = east
    west  = west
    north = north
    south = south
    
    from anuga.coordinate_transforms.redfearn import redfearn
    zone0, easting0, northing0 = redfearn(north,east)
    print zone0,easting0,northing0
    zone1, easting1, northing1 = redfearn(south,east)
    print zone1,easting1,northing1
    zone2, easting2, northing2 = redfearn(south,west)
    print zone2,easting2,northing2
    zone3, easting3, northing3 = redfearn(north,west)
    print zone3,easting3,northing3
    
    assert zone0==zone1==zone2==zone3
    
    ferret_polygon=[[easting0, northing0],[easting1, northing1],[easting2, northing2],[easting3, northing3]]
    points_outside=[]
    from anuga.utilities.polygon import is_inside_polygon
    for i,point in enumerate(bounding_polygon):
        if not is_inside_polygon(point, ferret_polygon, closed=True, verbose=False):
            print 'point %d = (%f,%f) not inside boundary'%(i,point[0],point[1])
            points_outside.append(point)

    from pylab import plot,show
    ferret_polygon=ensure_numeric(ferret_polygon)
    points_outside=ensure_numeric(points_outside)
    plot(ferret_polygon[:,0],ferret_polygon[:,1])
    plot(points_outside[:,0],points_outside[:,1],'o')
    show()
    print ferret_polygon[0]
    print ferret_polygon[1]
    print ferret_polygon[2]
    print ferret_polygon[3]
    print points_outside