Index: /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
===================================================================
--- /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 5360)
+++ /anuga_core/source/anuga/abstract_2d_finite_volumes/util.py	(revision 5361)
@@ -19,5 +19,5 @@
 from anuga.geospatial_data.geospatial_data import ensure_absolute
 from math import sqrt, atan, degrees
-from exceptions import IOError
+
 
 
@@ -110,5 +110,4 @@
     # Instead we pass in those attributes that are needed (and return them
     # if modified)
-
     kwargs = {'quantities': quantities,
               'interpolation_points': interpolation_points,
@@ -179,8 +178,8 @@
     try:
         fid = open(filename)
-    except IOError, e:
+    except Exception, e:
         msg = 'File "%s" could not be opened: Error="%s"'\
                   %(filename, e)
-        raise IOError, msg # So IOErrors can be caught
+        raise msg
 
     line = fid.readline()
@@ -244,5 +243,5 @@
         raise Exception, msg
 
-
+ 
     if interpolation_points is not None:
         interpolation_points = ensure_absolute(interpolation_points)
@@ -250,6 +249,4 @@
         assert interpolation_points.shape[1] == 2, msg
 
-    if verbose:
-        print 'File_function using quantities %s from file %s' %(str(quantity_names), filename)
 
     # Now assert that requested quantitites (and the independent ones)
@@ -273,10 +270,15 @@
 
     if filename[-3:] == 'tms' and spatial is True:
-        msg = 'Files of type tms must not contain spatial information'
+        msg = 'Files of type tms must not contain spatial  information'
         raise msg
 
     if filename[-3:] == 'sww' and spatial is False:
         msg = 'Files of type sww must contain spatial information'        
-        raise msg        
+        raise msg
+
+    if filename[-3:] == 'sts' and spatial is False:
+        #What if mux file only contains one point
+        msg = 'Files of type sts must contain spatial information'        
+        raise msg
 
     # Get first timestep
@@ -300,5 +302,6 @@
         x = fid.variables['x'][:]
         y = fid.variables['y'][:]
-        triangles = fid.variables['volumes'][:]
+        if filename[-3:] == 'sww':
+            triangles = fid.variables['volumes'][:]
 
         x = reshape(x, (len(x),1))
@@ -308,11 +311,6 @@
         if interpolation_points is not None:
             # Adjust for georef
-
-            # FIXME (Ole): Use geo_reference.get_relative
             interpolation_points[:,0] -= xllcorner
             interpolation_points[:,1] -= yllcorner        
-        
-
-
 
     if domain_starttime is not None:
@@ -350,6 +348,8 @@
 
     if not spatial:
-        vertex_coordinates = triangles = interpolation_points = None         
-
+        vertex_coordinates = triangles = interpolation_points = None
+    if filename[-3:] == 'sts':#added
+        triangles = None
+        #vertex coordinates is position of urs gauges
 
     # Return Interpolation_function instance as well as
@@ -782,5 +782,4 @@
         #print 'label', label
         leg_label.append(label)
-
 
         
@@ -2145,8 +2144,8 @@
         
         callable_sww = file_function(sww_file,
-                                     quantities=core_quantities,
-                                     interpolation_points=points_array,
-                                     verbose=verbose,
-                                     use_cache=use_cache)
+                                 quantities=core_quantities,
+                                 interpolation_points=points_array,
+                                 verbose=verbose,
+                                 use_cache=use_cache)
 
     gauge_file='gauge_'
Index: /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py	(revision 5360)
+++ /anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py	(revision 5361)
@@ -92,29 +92,33 @@
             max_vertices_per_cell = MAX_VERTICES_PER_CELL
         if mesh is None:
-            # Fixme (DSG) Throw errors if triangles or vertex_coordinates
-            # are None
+            if vertex_coordinates is not None and  triangles is not None:
+                # Fixme (DSG) Throw errors if triangles or vertex_coordinates
+                # are None
             
-            #Convert input to Numeric arrays
-            triangles = ensure_numeric(triangles, Int)
-            vertex_coordinates = ensure_absolute(vertex_coordinates,
+                #Convert input to Numeric arrays
+                triangles = ensure_numeric(triangles, Int)
+                vertex_coordinates = ensure_absolute(vertex_coordinates,
                                                  geo_reference = mesh_origin)
 
-            if verbose: print 'FitInterpolate: Building mesh'        
-            self.mesh = Mesh(vertex_coordinates, triangles)
-            #self.mesh.check_integrity() # Time consuming
+                if verbose: print 'FitInterpolate: Building mesh'        
+                self.mesh = Mesh(vertex_coordinates, triangles)
+                #self.mesh.check_integrity() # Time consuming
+            else:
+                self.mesh = None
         else:
             self.mesh = mesh
+
+        if self.mesh is not None:
+            if verbose: print 'FitInterpolate: Building quad tree'
+            #This stores indices of vertices
+            t0 = time.time()
+            #print "self.mesh.get_extent(absolute=True)", \
+            #self.mesh.get_extent(absolute=True)
+            self.root = build_quadtree(self.mesh,
+                                       max_points_per_cell = max_vertices_per_cell)
+            #print "self.root",self.root.show()
         
-        if verbose: print 'FitInterpolate: Building quad tree'
-        # This stores indices of vertices
-        t0 = time.time()
-        #print "self.mesh.get_extent(absolute=True)", \
-              #self.mesh.get_extent(absolute=True)
-        self.root = build_quadtree(self.mesh,
-                                   max_points_per_cell = max_vertices_per_cell)
-        #print "self.root",self.root.show()
-        
-        build_quadtree_time =  time.time()-t0
-        set_last_triangle()
+            build_quadtree_time =  time.time()-t0
+            set_last_triangle()
         
     def __repr__(self):
Index: /anuga_core/source/anuga/fit_interpolate/interpolate.py
===================================================================
--- /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 5360)
+++ /anuga_core/source/anuga/fit_interpolate/interpolate.py	(revision 5361)
@@ -89,4 +89,70 @@
                                 verbose=verbose,
                                 max_vertices_per_cell=max_vertices_per_cell)
+
+    def interpolate_polyline(self,
+                             f,
+                             vertex_coordinates,
+                             point_coordinates=None,
+                             start_blocking_len=500000,
+                             verbose=False):
+ 
+        if isinstance(point_coordinates, Geospatial_data):
+            point_coordinates = point_coordinates.get_data_points( \
+                absolute = True)
+ 
+        from utilities.polygon import point_on_line,point_on_line_py
+        from Numeric import ones
+        z=ones(len(point_coordinates),Float)
+
+        msg='point coordinates are not given (interpolate.py)'
+        assert point_coordinates is not None, msg
+        msg='function value must be specified at every interpolation node'
+        assert f.shape[0]==vertex_coordinates.shape[0],msg
+        msg='Must define function value at one or more nodes'
+        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
+                    
+                
+                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
 
     # FIXME: What is a good start_blocking_len value?
@@ -529,7 +595,7 @@
             vertex_coordinates = ensure_absolute(vertex_coordinates)
 
-            assert triangles is not None, 'Triangles array must be specified'
-            triangles = ensure_numeric(triangles)
-            self.spatial = True            
+            if triangles is not None:
+                triangles = ensure_numeric(triangles)
+            self.spatial = True          
 
         # Thin timesteps if needed
@@ -557,7 +623,8 @@
         # Precomputed spatial interpolation if requested
         if interpolation_points is not None:
-            if self.spatial is False:
-                raise 'Triangles and vertex_coordinates must be specified'
-            
+            #no longer true. sts files have spatial = True but
+            #if self.spatial is False:
+            #    raise 'Triangles and vertex_coordinates must be specified'
+            # 
             try:
 	        self.interpolation_points = interpolation_points = ensure_numeric(interpolation_points)
@@ -569,49 +636,55 @@
                 raise msg
 
-
-            # Check that all interpolation points fall within
-            # mesh boundary as defined by triangles and vertex_coordinates.
-            from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
-            from anuga.utilities.polygon import outside_polygon            
-
-            # Create temporary mesh object from mesh info passed
-            # into this function. 
-            mesh = Mesh(vertex_coordinates, triangles)
-            mesh_boundary_polygon = mesh.get_boundary_polygon()
-
-            
-            indices = outside_polygon(interpolation_points,
-                                      mesh_boundary_polygon)
-
-            # Record result
-            #self.mesh_boundary_polygon = mesh_boundary_polygon
-            self.indices_outside_mesh = indices
-
-            # Report
-            if len(indices) > 0:
-                msg = 'Interpolation points in Interpolation function fall ' 
-                msg += 'outside specified mesh. '
-                msg += 'Offending points:\n'
-                out_interp_pts = []
-                for i in indices:
-                    msg += '%d: %s\n' %(i, interpolation_points[i])
-                    out_interp_pts.append(ensure_numeric(interpolation_points[i]))
-
-                if verbose is True:
-                    import sys
-                    if sys.platform == 'win32':
-                        from anuga.utilities.polygon import plot_polygons
-                        #out_interp_pts = take(interpolation_points,[indices])
-                        title = 'Interpolation points fall outside specified mesh'
-                        plot_polygons([mesh_boundary_polygon,interpolation_points,out_interp_pts],
-                                      ['line','point','outside'],figname='points_boundary_out',label=title,verbose=verbose)
-
-                # Joaquim Luis suggested this as an Exception, so
-                # that the user can now what the problem is rather than
-                # looking for NaN's. However, NANs are handy as they can
-                # be ignored leaving good points for continued processing.
-                if verbose:
-                    print msg
-                #raise Exception(msg)
+            if triangles is not None and vertex_coordinates is not None:
+                # Check that all interpolation points fall within
+                # mesh boundary as defined by triangles and vertex_coordinates.
+                from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh
+                from anuga.utilities.polygon import outside_polygon            
+
+                # Create temporary mesh object from mesh info passed
+                # into this function. 
+                mesh = Mesh(vertex_coordinates, triangles)
+                mesh_boundary_polygon = mesh.get_boundary_polygon()
+
+            
+                indices = outside_polygon(interpolation_points,
+                                          mesh_boundary_polygon)
+
+                # Record result
+                #self.mesh_boundary_polygon = mesh_boundary_polygon
+                self.indices_outside_mesh = indices
+
+                # Report
+                if len(indices) > 0:
+                    msg = 'Interpolation points in Interpolation function fall ' 
+                    msg += 'outside specified mesh. '
+                    msg += 'Offending points:\n'
+                    out_interp_pts = []
+                    for i in indices:
+                        msg += '%d: %s\n' %(i, interpolation_points[i])
+                        out_interp_pts.append(ensure_numeric(interpolation_points[i]))
+
+                    if verbose is True:
+                        import sys
+                        if sys.platform == 'win32':
+                            from anuga.utilities.polygon import plot_polygons
+                            #out_interp_pts = take(interpolation_points,[indices])
+                            title = 'Interpolation points fall outside specified mesh'
+                            plot_polygons([mesh_boundary_polygon,interpolation_points,out_interp_pts],
+                                          ['line','point','outside'],figname='points_boundary_out',label=title,verbose=verbose)
+
+                    # Joaquim Luis suggested this as an Exception, so
+                    # that the user can now what the problem is rather than
+                    # looking for NaN's. However, NANs are handy as they can
+                    # be ignored leaving good points for continued processing.
+                    if verbose:
+                        print msg
+                    #raise Exception(msg)
+            elif triangles is None and vertex_coordinates is not None:#jj
+                #Dealing with sts file
+                pass
+            else:
+                msg = 'Sww file function requires both triangles and vertex_coordinates. sts file file function requires the later.'
+                raise Exception(msg)
 
             # Plot boundary and interpolation points
@@ -632,7 +705,11 @@
             # Build interpolator
             if verbose:
-                msg = 'Building interpolation matrix from source mesh '
-                msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
-                                                       triangles.shape[0])
+                if triangles is not None and vertex_coordinates is not None:
+                    msg = 'Building interpolation matrix from source mesh '
+                    msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0],
+                                                           triangles.shape[0])
+                elif triangles is None and vertex_coordinates is not None:
+                    msg = 'Building interpolation matrix from source points'
+                
                 print msg
 
@@ -665,9 +742,12 @@
                         print '    quantity %s, size=%d' %(name, len(Q))
                         
-                    # Interpolate    
-                    result = interpol.interpolate(Q,
-                                                  point_coordinates=\
-                                                  self.interpolation_points,
-                                                  verbose=False) # Don't clutter
+                    # Interpolate 
+                    if triangles is not None and vertex_coordinates is not None:   
+                        result = interpol.interpolate(Q,
+                                                      point_coordinates=\
+                                                      self.interpolation_points,
+                                                      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)
 
                     #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 5360)
+++ /anuga_core/source/anuga/shallow_water/test_data_manager.py	(revision 5361)
@@ -6289,15 +6289,4 @@
 
     def test_urs2sts(self):
-        """tide = 1
-        time_step_count = 3
-        time_step = 2
-        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
-        depth=20
-        ha=2
-        ua=5
-        va=-10 #-ve added to take into account mux file format where south
-               # is positive.
-        """
-        #tide = 1
         tide=0
         time_step_count = 3
@@ -6407,4 +6396,77 @@
         self.delete_mux(files)
         os.remove(sts_file)
+
+    def test_file_boundary_sts(self):
+        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.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
+        tide=0.
+        time_step_count = 5
+        time_step = 2
+        lat_long_points =bounding_polygon[0:3]
+        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,1], 'otherocean': [2,3,4]}
+        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)
+        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 temp_fbound-temp_drchlt<epsilon
+        os.remove(sts_file+'.sts')
+        os.remove(meshname)
 
     def test_lon_lat2grid(self):
