def create_sts_boundary(filename,verbose=False):
    from Scientific.IO.NetCDF import NetCDFFile
    from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
        compress,zeros,fabs
    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.

    #if verbose:
         #from pylab import plot,show
         #plot(x,y,'o')
         #plot(boundary[:,0],boundary[:,1])
         #show()
    return 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)
bounding_polygon=boundary.tolist()
del bounding_polygon[0:2] #delete first two elements
n=len(bounding_polygon)
del bounding_polygon[n-8:n]

import project
bp=project.bounding_polygon
bounding_polygon.insert(0,bp[0])
bounding_polygon.extend(bp[len(bp)-5:len(bp)])
from anuga.utilities.numerical_tools import ensure_numeric
bounding_polygon=ensure_numeric(bounding_polygon)
#from pylab import plot,show
#plot(bounding_polygon[:,0],bounding_polygon[:,1])
#bp=ensure_numeric(bp)
#plot(bp[:,0],bp[:,1])
#show()

check_ferret_extent=False
if check_ferret_extent:
    east  = project.east
    west  = project.west
    north = project.north
    south = project.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