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

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,allclose,ones
from anuga.utilities.polygon import inside_polygon,read_polygon
from os import sep
from time import localtime, strftime, gmtime
from anuga.shallow_water.data_manager import urs2sts,create_sts_boundary
import os

def create_sts_boundary_order_file(urs_filename,order_filename,sts_filename,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.
    """
    from anuga.shallow_water.data_manager import read_mux2_py
    times, y, x, elevation, quantity, starttime = read_mux2_py([urs_filename],weights=ones(1,Float))
    
    #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)
    order=[]
    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]
         order.append(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.
    
    boundary=ensure_numeric(boundary)
    from anuga.coordinate_transforms import redfearn,convert_from_latlon_to_utm
    boundary_utm,zone=convert_from_latlon_to_utm(longitudes=list(boundary[:,0]),latitudes=list(boundary[:,1]))

    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,boundary_utm

######################
# 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'
boundary_most_out = dir+sep+'tide_most_boundary_condition'

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

######################################
# Create urs boundary from mux2files #
######################################
urs_filename='polyline'
base_name='tide_polyline'
order_filename='boxingday_boundary_order.txt'

# Create ordering file
# Usually this would be provided before hand
boundary,boundary_utm=create_sts_boundary_order_file(urs_filename+'-z-mux2',order_filename,base_name,verbose=False)

# Create ordered sts file
#if not os.path.exists(base_name+'.sts'):
print 'creating sts file'
urs2sts(urs_filename,basename_out=base_name,
        ordering_filename=order_filename,
        mean_stage=tide,
        verbose=False)
#else:
#    print 'sts file exists'

# Read in boundary from ordered sts file
urs_boundary=create_sts_boundary(base_name)

# Checks that the boundary read from the sts file is the same as the one
# defined by the function create_sts_boundary_order_file
assert allclose(boundary_utm,urs_boundary)

#############################################################
# Add urs boundary segment to user defined boundary polygon #
#############################################################

# Read in specific user defined boundary
dir='mesh_polygons'
extent = 'extent.csv'
user_boundary_polygon = read_polygon(dir+sep+extent)

# The following is specific to this scenario
indices=inside_polygon(urs_boundary, user_boundary_polygon, 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])

# Append clipped sts boundary to user defined domain boundary
bounding_polygon.extend(user_boundary_polygon[2:len(user_boundary_polygon)])

###############################
# 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 urs boundary and internal regions #
##########################################

plot=False
if plot:
    from pylab import plot,show,hold
    bounding_polygon=ensure_numeric(bounding_polygon)#for plotting only
    #boundary=ensure_numeric(boundary)
    #plot(boundary[:,0],boundary[:,1],'o')
    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()

    
#######################
# For MOST SIMULATION #
#######################
east  = 97.7
west  = 96.7
north = 9.0
south = 5.9


bay_west = 417348.0
bay_east = 425656.0
bay_south = 870565.0
bay_north = 877008.0

bay_grid_poly=[[bay_west,bay_south],[bay_west,bay_north],[bay_east,bay_north],[bay_east,bay_south]]

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