Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 5371)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 5372)
@@ -15,5 +15,5 @@
 
 from anuga.utilities.numerical_tools import ensure_numeric
-from Numeric import arange, choose, zeros, Float, array
+from Numeric import arange, choose, zeros, Float, array, allclose, take, compress
     
 from anuga.geospatial_data.geospatial_data import ensure_absolute
@@ -33,5 +33,6 @@
                   time_thinning=1,
                   verbose=False,
-                  use_cache=False):
+                  use_cache=False,
+                  boundary_polygon=None):
     """Read time history of spatial data from NetCDF file and return
     a callable object.
@@ -114,5 +115,6 @@
               'domain_starttime': domain_starttime,
               'time_thinning': time_thinning,                   
-              'verbose': verbose}
+              'verbose': verbose,
+              'boundary_polygon':boundary_polygon}
 
 
@@ -130,5 +132,6 @@
                              dependencies=[filename],
                              compression=False,                  
-                             verbose=verbose)
+                             verbose=verbose,
+                             boundary_polygon=boundary_polygon)
 
     else:
@@ -166,5 +169,6 @@
                    domain_starttime=None,
                    time_thinning=1,
-                   verbose=False):
+                   verbose=False,
+                   boundary_polygon=None):
     """Internal function
     
@@ -193,5 +197,6 @@
                                         domain_starttime,
                                         time_thinning=time_thinning,
-                                        verbose=verbose)
+                                        verbose=verbose,
+                                        boundary_polygon=boundary_polygon)
     else:
         raise 'Must be a NetCDF File'
@@ -204,5 +209,6 @@
                              domain_starttime=None,                            
                              time_thinning=1,                             
-                             verbose=False):
+                             verbose=False,
+                             boundary_polygon=None):
     """Read time history of spatial data from NetCDF sww file and
     return a callable object f(t,x,y)
@@ -282,4 +288,9 @@
         raise msg
 
+    if filename[-3:] == 'sts' and boundary_polygon is None:
+        #What if mux file only contains one point
+        msg = 'Files of type sts require boundary polygon'        
+        raise msg
+
     # Get first timestep
     try:
@@ -308,4 +319,44 @@
         y = reshape(y, (len(y),1))
         vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
+
+        if boundary_polygon is not None:
+            #removes sts points that do not lie on boundary
+            boundary_polygon=ensure_numeric(boundary_polygon)
+            boundary_polygon[:,0] -= xllcorner
+            boundary_polygon[:,1] -= yllcorner
+            temp=[]
+            boundary_id=[]
+            gauge_id=[]
+            for i in range(len(boundary_polygon)):
+                for j in range(len(x)):
+                    if allclose(vertex_coordinates[j],boundary_polygon[i],1e-4):
+                        #FIX ME:
+                        #currently gauges lat and long is stored as float and
+                        #then cast to double. This cuases slight repositioning
+                        #of vertex_coordinates.
+                        temp.append(boundary_polygon[i])
+                        gauge_id.append(j)
+                        boundary_id.append(i)
+                        break
+            gauge_neighbour_id=[]
+            for i in range(len(boundary_id)-1):
+                if boundary_id[i]+1==boundary_id[i+1]:
+                    gauge_neighbour_id.append(i+1)
+                else:
+                    gauge_neighbour_id.append(-1)
+            if boundary_id[len(boundary_id)-1]==len(boundary_polygon)-1 and boundary_id[0]==0:
+                gauge_neighbour_id.append(0)
+            else:
+                gauge_neighbour_id.append(-1)
+            gauge_neighbour_id=ensure_numeric(gauge_neighbour_id)
+            if len(compress(gauge_neighbour_id>=0,gauge_neighbour_id))!=len(temp)-1:
+                msg='incorrect number of segments'
+                raise msg
+            vertex_coordinates=ensure_numeric(temp)
+            if len(vertex_coordinates)==0:
+                msg = 'None of the sts gauges fall on the boundary'
+                raise msg
+        else:
+            gauge_neighbour_id=None
 
         if interpolation_points is not None:
@@ -342,6 +393,8 @@
     for i, name in enumerate(quantity_names):
         quantities[name] = fid.variables[name][:]
+        if boundary_polygon is not None:
+            #removes sts points that do not lie on boundary
+            quantities[name] = take(quantities[name],gauge_id,1)
     fid.close()
-    
 
     from anuga.fit_interpolate.interpolate import Interpolation_function
@@ -362,5 +415,5 @@
                                   interpolation_points,
                                   time_thinning=time_thinning,
-                                  verbose=verbose), starttime
+                                  verbose=verbose,gauge_neighbour_id=gauge_neighbour_id), starttime
 
     # NOTE (Ole): Caching Interpolation function is too slow as
Index: /anuga_core/source/anuga/fit_interpolate/interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 5371)
+++ /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 5372)
@@ -93,6 +93,6 @@
                              f,
                              vertex_coordinates,
+                             gauge_neighbour_id,
                              point_coordinates=None,
-                             start_blocking_len=500000,
                              verbose=False):
  
@@ -112,46 +112,30 @@
         assert f.shape[0]>0,msg
 
-        #print point_on_line_py(point_coordinates[3],[vertex_coordinates[0],vertex_coordinates[1]])
-        #print vertex_coordinates[0],vertex_coordinates[1]
-
-        #print point_coordinates
-        #print vertex_coordinates
         n=f.shape[0]
         if n==1:
             z=f*z
         elif n>1:
-            index=0#index of segment on which last point was found
-            k=0#number of points on line
             for i in range(len(point_coordinates)):
                 found = False
-                #find the segment the cetnroid lies on
-                #Start with the segment the last point was found on
-                #For n points there will be n-1 segments
-                for j in range(n-1):
-                    #print 'searcing segment', index+j
-                    #reached last segment look at first segment
-                    if (index+j)==n-2:
-                        index=0
-                    if point_on_line_py(point_coordinates[i],[vertex_coordinates[(j)+index],vertex_coordinates[(j+1)+index]]):
-                        found=True
-                        x0=vertex_coordinates[j][0];y0=vertex_coordinates[j][1]
-                        x1=vertex_coordinates[j+1][0];y1=vertex_coordinates[j+1][1]
-                        x2=point_coordinates[i][0];y2=point_coordinates[i][1]
-                        
-                        segment_len=sqrt((x1-x0)**2+(y1-y0)**2)
-                        dist=sqrt((x2-x0)**2+(y2-y0)**2)
-                        z[i]=(f[j+1]-f[j])/segment_len*dist+f[j]
-                        #print 'element found on segment',j+index
-                        index=j
-                        break
-                    
-                
+                for j in range(n):
+                    if gauge_neighbour_id[j]>=0:
+                        if point_on_line_py(point_coordinates[i],[vertex_coordinates[j],vertex_coordinates[gauge_neighbour_id[j]]]):
+                            found=True
+                            x0=vertex_coordinates[j][0]
+                            y0=vertex_coordinates[j][1]
+                            x1=vertex_coordinates[gauge_neighbour_id[j]][0]
+                            y1=vertex_coordinates[gauge_neighbour_id[j]][1]
+                            x2=point_coordinates[i][0]
+                            y2=point_coordinates[i][1]
+                            
+                            segment_len=sqrt((x1-x0)**2+(y1-y0)**2)
+                            dist=sqrt((x2-x0)**2+(y2-y0)**2)
+                            z[i]=(f[gauge_neighbour_id[j]]-f[j])/segment_len*dist+f[j]
+                            #print 'element found on segment'
+                            break
+                                  
                 if not found:
                     z[i]=0.0
                     #print 'point not on urs boundary'
-                else:
-                    k+=1
-        #print 'number of midpoints on urs boundary',k
-        #print z
         return z
 
@@ -550,5 +534,6 @@
                  interpolation_points=None,
                  time_thinning=1,
-                 verbose=False):
+                 verbose=False,
+                 gauge_neighbour_id=None):
         """Initialise object and build spatial interpolation if required
 
@@ -749,5 +734,5 @@
                                                       verbose=False) # Don't clutter
                     elif triangles is None and vertex_coordinates is not None:
-                        result=interpol.interpolate_polyline(Q,vertex_coordinates,point_coordinates=self.interpolation_points)
+                        result=interpol.interpolate_polyline(Q,vertex_coordinates,gauge_neighbour_id,point_coordinates=self.interpolation_points)
 
                     #assert len(result), len(interpolation_points)
Index: /anuga_core/source/anuga/shallow_water/test_data_manager.py
===================================================================
--- /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 5371)
+++ /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 5372)
@@ -6441,5 +6441,5 @@
         domain_fbound = pmesh_to_domain_instance(meshname, Domain)
         domain_fbound.set_quantity('stage', tide)
-        Bf = File_boundary(sts_file+'.sts', domain_fbound)
+        Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=bounding_polygon)
         Br = Reflective_boundary(domain_fbound)
         Bd=Dirichlet_boundary([2.0,220,-220])
@@ -6466,5 +6466,85 @@
             i+=1
 
-        assert temp_fbound-temp_drchlt<epsilon
+        assert allclose(temp_fbound,temp_drchlt)
+        os.remove(sts_file+'.sts')
+        os.remove(meshname)
+
+    def test_file_boundary_sts2(self):
+        # mux2 file has points not included in boundary
+        # mux2 gauges are not stored with the same order as they are 
+        # found in bounding_polygon
+        from anuga.shallow_water import Domain
+        from anuga.shallow_water import Reflective_boundary
+        from anuga.shallow_water import Dirichlet_boundary
+        from anuga.shallow_water import File_boundary
+        from anuga.pmesh.mesh_interface import create_mesh_from_regions
+        from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
+        bounding_polygon=[[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02],[6.0,97.0]]
+        tide=0.
+        time_step_count = 5
+        time_step = 2
+        lat_long_points=bounding_polygon[0:2]
+        lat_long_points.insert(0,bounding_polygon[len(bounding_polygon)-1])
+        lat_long_points.insert(0,[6.0,97.01])
+        n=len(lat_long_points)
+        first_tstep=ones(n,Int)
+        last_tstep=(time_step_count)*ones(n,Int)
+        gauge_depth=20*ones(n,Float)
+        ha=2*ones((n,time_step_count),Float)
+        ua=10*ones((n,time_step_count),Float)
+        va=-10*ones((n,time_step_count),Float)
+        base_name, files = self.write_mux2(lat_long_points,
+                                   time_step_count, time_step,
+                                   first_tstep, last_tstep,
+                                   depth=gauge_depth,
+                                   ha=ha,
+                                   ua=ua,
+                                   va=va)
+
+        sts_file=base_name
+        urs2sts(base_name,sts_file,mean_stage=tide,verbose=False)
+        self.delete_mux(files)
+
+        #print 'start create mesh from regions'
+        for i in range(len(bounding_polygon)):
+            zone,bounding_polygon[i][0],bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],bounding_polygon[i][1])
+        extent_res=1000000
+        meshname = 'urs_test_mesh' + '.tsh'
+        interior_regions=None
+        boundary_tags={'ocean': [0,4], 'otherocean': [1,2,3]}
+        #have to change boundary tags from last example because now bounding
+        #polygon starts in different place.
+        create_mesh_from_regions(bounding_polygon,boundary_tags=boundary_tags,
+                         maximum_triangle_area=extent_res,filename=meshname,
+                         interior_regions=interior_regions,verbose=False)
+        
+        domain_fbound = pmesh_to_domain_instance(meshname, Domain)
+        domain_fbound.set_quantity('stage', tide)
+        Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=bounding_polygon)
+        Br = Reflective_boundary(domain_fbound)
+        Bd=Dirichlet_boundary([2.0,220,-220])
+        domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br})
+        finaltime=time_step*(time_step_count-1)
+        yieldstep=time_step
+        temp_fbound=zeros(int(finaltime/yieldstep)+1,Float)
+        i=0
+        for t in domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime, 
+                                      skip_initial_step = False):
+            temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2]
+            i+=1
+        
+        domain_drchlt = pmesh_to_domain_instance(meshname, Domain)
+        domain_drchlt.set_quantity('stage', tide)
+        Br = Reflective_boundary(domain_drchlt)
+        Bd=Dirichlet_boundary([2.0,220,-220])
+        domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br})
+        temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float)
+        i=0
+        for t in domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime, 
+                                      skip_initial_step = False):
+            temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2]
+            i+=1
+
+        assert allclose(temp_fbound,temp_drchlt)
         os.remove(sts_file+'.sts')
         os.remove(meshname)
Index: /anuga_core/source/anuga/utilities/polygon.py
===================================================================
--- /anuga_core/source/anuga/utilities/polygon.py	(revision 5371)
+++ /anuga_core/source/anuga/utilities/polygon.py	(revision 5372)
@@ -16,4 +16,46 @@
 from anuga.utilities.numerical_tools import ensure_numeric
 from anuga.geospatial_data.geospatial_data import ensure_absolute
+
+def point_on_line_py(point, line):
+    from Numeric import fabs
+    point = ensure_numeric(point)
+    line = ensure_numeric(line)
+
+
+    x=point[0];y=point[1]
+    x0=line[0,0];y0=line[0,1]
+    x1=line[1,0];y1=line[1,1]
+    #from pylab import plot,show
+    #plot(line[:,0],line[:,1])
+    #plot([x],[y],'o')
+    #show()
+    a0 = x - x0;
+    a1 = y - y0;
+
+    a_normal0 = a1;
+    a_normal1 = -a0;
+
+    b0 = x1 - x0;
+    b1 = y1 - y0;
+
+    #if ( a_normal0*b0 + a_normal1*b1 == 0.0 ):
+    eps=200
+    #print 'normal',a_normal0*b0 + a_normal1*b1
+    if ( fabs(a_normal0*b0 + a_normal1*b1) < eps ):
+        #for some reason (perhaps catastrophic cancellation) urs_boundary.py
+        #example only works for eps=2
+        #point is somewhere on the infinite extension of the line
+        # FIXME (Ole): Perhaps add a tolerance here instead of 0.0 
+
+        len_a = sqrt(a0*a0 + a1*a1);
+        len_b = sqrt(b0*b0 + b1*b1);
+
+        if (a0*b0 + a1*b1 >= 0 and len_a <= len_b):
+            return 1
+        else:
+            return 0
+    else:
+        return 0
+
 
 def point_on_line(point, line):
