Index: /trunk/anuga_core/source/anuga/file/mux.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/mux.py	(revision 7765)
+++ /trunk/anuga_core/source/anuga/file/mux.py	(revision 7766)
@@ -1,8 +1,16 @@
+"""
+    Read a mux2 file.
+"""
+
 from anuga.utilities.numerical_tools import ensure_numeric
 import numpy as num     
         
- ################################################################################
+################################################################################
 # READ MUX2 FILES line of points
 ################################################################################
+
+WAVEHEIGHT_MUX_LABEL = '-z-mux'
+EAST_VELOCITY_LABEL =  '-e-mux'
+NORTH_VELOCITY_LABEL =  '-n-mux'
 
 WAVEHEIGHT_MUX2_LABEL = '-z-mux2'
@@ -43,7 +51,7 @@
     # Convert verbose to int C flag
     if verbose:
-        verbose=1
+        verbose = 1
     else:
-        verbose=0
+        verbose = 0
 
     if weights is None:
@@ -102,5 +110,6 @@
     longitudes = num.zeros(number_of_selected_stations, num.float)
     elevation = num.zeros(number_of_selected_stations, num.float)
-    quantity = num.zeros((number_of_selected_stations, parameters_index), num.float)
+    quantity = num.zeros((number_of_selected_stations, parameters_index), \
+                                                    num.float)
 
     starttime = 1e16
@@ -116,28 +125,3 @@
 
 
-##
-# @brief ??
-# @param mux_times ??
-# @param mint ??
-# @param maxt ??
-# @return ??
-def read_time_from_mux(mux_times, mint, maxt):
-    """
-    """
 
-    if mint == None:
-        mux_times_start_i = 0
-    else:
-        mux_times_start_i = num.searchsorted(mux_times, mint)
-
-    if maxt == None:
-        mux_times_fin_i = len(mux_times)
-    else:
-        maxt += 0.5 # so if you specify a time where there is
-                    # data that time will be included
-        mux_times_fin_i = num.searchsorted(mux_times, maxt)
-
-    return mux_times_start_i, mux_times_fin_i
-
-
-       
Index: /trunk/anuga_core/source/anuga/file/sts.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/sts.py	(revision 7765)
+++ /trunk/anuga_core/source/anuga/file/sts.py	(revision 7766)
@@ -284,5 +284,5 @@
     except:
         msg = 'Cannot open %s' % stsname + '.sts'
-        raise msg
+        raise IOError(msg)
 
     xllcorner = fid.xllcorner[0]
Index: /trunk/anuga_core/source/anuga/file/test_mux.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/test_mux.py	(revision 7765)
+++ /trunk/anuga_core/source/anuga/file/test_mux.py	(revision 7766)
@@ -10,4 +10,7 @@
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 
+from mux import WAVEHEIGHT_MUX_LABEL, EAST_VELOCITY_LABEL, \
+                            NORTH_VELOCITY_LABEL
+
 from mux import WAVEHEIGHT_MUX2_LABEL, EAST_VELOCITY_MUX2_LABEL, \
                 NORTH_VELOCITY_MUX2_LABEL
@@ -15,6 +18,7 @@
 from mux import read_mux2_py
 from anuga.file_conversion.urs2sts import urs2sts
-
-class TestCase(unittest.TestCase):
+from anuga.file.urs import Read_urs
+
+class Test_Mux(unittest.TestCase):
     def setUp(self):
         pass
@@ -1487,10 +1491,45 @@
         
         self.delete_mux(files)
-        
+ 
+    def test_Urs_points(self):
+        time_step_count = 3
+        time_step = 2
+        lat_long_points =[(-21.5,114.5),(-21.5,115),(-21.,115)]
+        base_name, files = self.write_mux(lat_long_points,
+                                          time_step_count, time_step)
+        for file in files:
+            urs = Read_urs(file)
+            assert time_step_count == urs.time_step_count
+            assert time_step == urs.time_step
+
+            for lat_lon, dep in map(None, lat_long_points, urs.lonlatdep):
+                    _ , e, n = redfearn(lat_lon[0], lat_lon[1])
+                    assert num.allclose(n, dep[2])
+                        
+            count = 0
+            for slice in urs:
+                count += 1
+                #print slice
+                for lat_lon, quantity in map(None, lat_long_points, slice):
+                    _ , e, n = redfearn(lat_lon[0], lat_lon[1])
+                    #print "quantity", quantity
+                    #print "e", e
+                    #print "n", n
+                    if file[-5:] == WAVEHEIGHT_MUX_LABEL[-5:] or \
+                           file[-5:] == NORTH_VELOCITY_LABEL[-5:] :
+                        assert num.allclose(e, quantity)
+                    if file[-5:] == EAST_VELOCITY_LABEL[-5:]:
+                        assert num.allclose(n, quantity)
+            assert count == time_step_count
+                     
+        self.delete_mux(files)        
+
+
+
         
 ################################################################################
 
 if __name__ == "__main__":
-    suite = unittest.makeSuite(TestCase,'test')
+    suite = unittest.makeSuite(Test_Mux,'test')
     runner = unittest.TextTestRunner() #verbosity=2)
     runner.run(suite)
Index: /trunk/anuga_core/source/anuga/file/test_urs.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/test_urs.py	(revision 7766)
+++ /trunk/anuga_core/source/anuga/file/test_urs.py	(revision 7766)
@@ -0,0 +1,227 @@
+import unittest
+import numpy as num
+
+from urs import calculate_boundary_points, save_boundary_as_urs
+from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees
+from anuga.geospatial_data.geospatial_data import Geospatial_data
+
+class Test_Urs(unittest.TestCase):
+    def setUp(self):
+        self.verbose = True
+        pass
+
+    def tearDown(self):
+        pass
+        
+    def test_URS_points_needed(self):
+        
+        ll_lat = -21.5
+        ll_long = 114.5
+        grid_spacing = 1./60.
+        lat_amount = 30
+        long_amount = 30
+        zone = 50
+
+        boundary_polygon = [[250000,7660000],[280000,7660000],
+                             [280000,7630000],[250000,7630000]]
+        geo=calculate_boundary_points(boundary_polygon, zone, 
+                              ll_lat, ll_long, grid_spacing, 
+                              lat_amount, long_amount,
+                              verbose=self.verbose)
+        # to test this geo, can info from it be transfered onto the boundary
+        # poly?
+        #Maybe see if I can fit the data to the polygon - have to represent
+        # the poly as points though.
+        #geo.export_points_file("results.txt", as_lat_long=True)
+        results = frozenset(geo.get_data_points(as_lat_long=True))
+        #print 'results',results
+
+        # These are a set of points that have to be in results
+        points = []
+        for i in range(18):
+            lat = -21.0 - 8./60 - grid_spacing * i
+            points.append((lat,degminsec2decimal_degrees(114,35,0))) 
+            points.append((lat,degminsec2decimal_degrees(114,36,0))) 
+            points.append((lat,degminsec2decimal_degrees(114,52,0))) 
+            points.append((lat,degminsec2decimal_degrees(114,53,0)))
+        geo_answer = Geospatial_data(data_points=points,
+                                     points_are_lats_longs=True) 
+        #geo_answer.export_points_file("answer.txt", as_lat_long=True)  
+        answer = frozenset(points)
+        
+        outs = answer.difference(results)
+        #print "outs", outs
+        # This doesn't work.  Though visualising the results shows that
+        # it is correct.
+        #assert answer.issubset(results)
+        # this is why;
+        #point (-21.133333333333333, 114.58333333333333)
+        #result (-21.133333332232368, 114.58333333300342)
+        
+        for point in points:
+            found = False
+            for result in results:
+                if num.allclose(point, result):
+                    found = True
+                    break
+            if not found:
+                assert False
+        
+    
+    def dave_test_URS_points_needed(self):
+        ll_lat = -21.51667
+        ll_long = 114.51667
+        grid_spacing = 2./60.
+        lat_amount = 15
+        long_amount = 15
+
+       
+        boundary_polygon = [[250000,7660000],[280000,7660000],
+                             [280000,7630000],[250000,7630000]]
+        calculate_boundary_points('a_test_example',boundary_polygon,
+                                  ll_lat, ll_long, grid_spacing, 
+                                  lat_amount, long_amount,
+                                  verbose=self.verbose)
+        
+    def X_test_URS_points_neededII(self):
+        ll_lat = -21.5
+        ll_long = 114.5
+        grid_spacing = 1./60.
+        lat_amount = 30
+        long_amount = 30
+
+        # change this so lats and longs are inputed, then converted
+        
+        #boundary_polygon = [[7660000,250000],[7660000,280000],
+        #                     [7630000,280000],[7630000,250000]]
+        calculate_boundary_points(boundary_polygon, ll_lat, ll_long,
+                          grid_spacing, 
+                          lat_amount, long_amount,
+                          verbose=self.verbose)
+        
+    def test_URS_points_northern_hemisphere(self):
+               
+        LL_LAT = 8.0
+        LL_LONG = 97.0
+        GRID_SPACING = 2.0/60.0
+        LAT_AMOUNT = 2
+        LONG_AMOUNT = 2
+        ZONE = 47
+
+        # 
+        points = []
+        for i in range(2):
+            for j in range(2):
+                points.append((degminsec2decimal_degrees(8,1+i*2,0),
+                               degminsec2decimal_degrees(97,1+i*2,0)))
+        #print "points", points
+        geo_poly = Geospatial_data(data_points=points,
+                                     points_are_lats_longs=True)
+        poly_lat_long = geo_poly.get_data_points(as_lat_long=False,
+                                       isSouthHemisphere=False)
+        #print "seg_lat_long",  poly_lat_long
+        
+      #   geo=URS_points_needed_to_file('test_example_poly3', poly_lat_long,
+#                                   ZONE,
+#                                   LL_LAT, LL_LONG,
+#                                   GRID_SPACING,
+#                                   LAT_AMOUNT, LONG_AMOUNT,
+#                                   isSouthernHemisphere=False,
+#                                   export_csv=True,
+#                                   verbose=self.verbose)
+        
+        geo=calculate_boundary_points(poly_lat_long,
+                                  ZONE,
+                                  LL_LAT, LL_LONG,
+                                  GRID_SPACING,
+                                  LAT_AMOUNT, LONG_AMOUNT,
+                                  isSouthHemisphere=False,
+                                  verbose=self.verbose)
+        
+        results = frozenset(geo.get_data_points(as_lat_long=True,
+                                                isSouthHemisphere=False))
+        #print 'results',results
+
+        # These are a set of points that have to be in results
+        points = [] 
+        for i in range(2):
+            for j in range(2):
+                points.append((degminsec2decimal_degrees(8,i*2,0),
+                               degminsec2decimal_degrees(97,i*2,0)))
+        #print "answer points", points
+        answer = frozenset(points)
+        
+        for point in points:
+            found = False
+            for result in results:
+                if num.allclose(point, result):
+                    found = True
+                    break
+            if not found:
+                assert False
+        
+
+    def covered_in_other_tests_test_URS_points_needed_poly1(self):
+        # Values used for FESA 2007 results
+        # domain in southern hemisphere zone 51        
+        LL_LAT = -50.0
+        LL_LONG = 80.0
+        GRID_SPACING = 2.0/60.0
+        LAT_AMOUNT = 4800
+        LONG_AMOUNT = 3600
+        ZONE = 51
+        
+        poly1 = [[296361.89, 8091928.62],
+                 [429495.07,8028278.82],
+                 [447230.56,8000674.05],
+                 [429661.2,7982177.6],
+                 [236945.9,7897453.16],
+                 [183493.44,7942782.27],
+                 [226583.04,8008058.96]]
+
+        save_boundary_as_urs('test_example_poly2', poly1,
+                                  ZONE,
+                                  LL_LAT, LL_LONG,
+                                  GRID_SPACING,
+                                  LAT_AMOUNT, LONG_AMOUNT,
+                                  verbose=self.verbose)
+        
+
+
+    def covered_in_other_tests_test_URS_points_needed_poly2(self):
+        # Values used for 2004 validation work
+        # domain in northern hemisphere zone 47        
+        LL_LAT = 0.0
+        LL_LONG = 90.0
+        GRID_SPACING = 2.0/60.0
+        LAT_AMOUNT = (15-LL_LAT)/GRID_SPACING
+        LONG_AMOUNT = (100-LL_LONG)/GRID_SPACING 
+        ZONE = 47
+        
+        poly2 = [[419336.424,810100.845],
+                 [342405.0281,711455.8026],
+                 [274649.9152,723352.9603],
+                 [272089.092,972393.0131],
+                 [347633.3754,968551.7784],
+                 [427979.2022,885965.2313],
+                 [427659.0993,875721.9386],
+                 [429259.6138,861317.3083],
+                 [436301.8775,840830.723]]
+        
+        save_boundary_as_urs('test_example_poly2', poly2,
+                                  ZONE,
+                                  LL_LAT, LL_LONG,
+                                  GRID_SPACING,
+                                  LAT_AMOUNT, LONG_AMOUNT,
+                                  isSouthernHemisphere=False,
+                                  verbose=self.verbose) 
+        
+
+
+        
+################################################################################
+
+if __name__ == "__main__":
+    suite = unittest.makeSuite(Test_Urs,'test')
+    runner = unittest.TextTestRunner() #verbosity=2)
+    runner.run(suite)
Index: /trunk/anuga_core/source/anuga/file/urs.py
===================================================================
--- /trunk/anuga_core/source/anuga/file/urs.py	(revision 7765)
+++ /trunk/anuga_core/source/anuga/file/urs.py	(revision 7766)
@@ -1,4 +1,14 @@
-##
-# @brief 
+from struct import pack, unpack
+import array as p_array
+import numpy as num
+
+
+from anuga.coordinate_transforms.geo_reference import Geo_reference
+
+from anuga.geospatial_data.geospatial_data import ensure_absolute, \
+                                                    Geospatial_data
+
+from anuga.coordinate_transforms.redfearn import redfearn
+
 class Read_urs:
     """
@@ -88,2 +98,284 @@
         self.mux_file.close()
    
+
+
+
+### PRODUCING THE POINTS NEEDED FILE ###
+
+# Ones used for FESA 2007 results
+#LL_LAT = -50.0
+#LL_LONG = 80.0
+#GRID_SPACING = 1.0/60.0
+#LAT_AMOUNT = 4800
+#LONG_AMOUNT = 3600
+
+
+##
+# @brief 
+# @param file_name 
+# @param boundary_polygon 
+# @param zone 
+# @param ll_lat 
+# @param ll_long 
+# @param grid_spacing 
+# @param lat_amount 
+# @param long_amount 
+# @param isSouthernHemisphere 
+# @param export_csv 
+# @param use_cache 
+# @param verbose True if this function is to be verbose.
+# @return 
+
+def save_boundary_as_urs(file_name, boundary_polygon, zone,
+                              ll_lat, ll_long,
+                              grid_spacing,
+                              lat_amount, long_amount,
+                              isSouthernHemisphere=True,
+                              export_csv=False, use_cache=False,
+                              verbose=False):
+    """
+    Given the info to replicate the URS grid and a polygon output
+    a file that specifies the cloud of boundary points for URS.
+
+    This creates a .urs file.  This is in the format used by URS;
+    1st line is the number of points,
+    each line after represents a point,in lats and longs.
+
+    Note: The polygon cannot cross zones or hemispheres.
+
+    A work-a-round for different zones or hemispheres is to run this twice,
+    once for each zone, and then combine the output.
+
+    file_name - name of the urs file produced for David.
+    boundary_polygon - a list of points that describes a polygon.
+                      The last point is assumed ot join the first point.
+                      This is in UTM (lat long would be better though)
+
+     This is info about the URS model that needs to be inputted.
+
+    ll_lat - lower left latitude, in decimal degrees
+    ll-long - lower left longitude, in decimal degrees
+    grid_spacing - in deciamal degrees
+    lat_amount - number of latitudes
+    long_amount- number of longs
+
+    Don't add the file extension.  It will be added.
+    """
+
+    geo = calculate_boundary_points(boundary_polygon, zone, ll_lat, ll_long,
+                            grid_spacing,
+                            lat_amount, long_amount, isSouthernHemisphere,
+                            use_cache, verbose)
+
+    if not file_name[-4:] == ".urs":
+        file_name += ".urs"
+
+    geo.export_points_file(file_name, isSouthHemisphere=isSouthernHemisphere)
+
+    if export_csv:
+        if file_name[-4:] == ".urs":
+            file_name = file_name[:-4] + ".csv"
+        geo.export_points_file(file_name)
+
+    return geo
+
+
+##
+# @brief 
+# @param boundary_polygon
+# @param zone
+# @param ll_lat
+# @param ll_long
+# @param grid_spacing
+# @param lat_amount
+# @param long_amount
+# @param isSouthHemisphere
+# @param use_cache
+# @param verbose
+def calculate_boundary_points(boundary_polygon, zone, ll_lat,
+                      ll_long, grid_spacing,
+                      lat_amount, long_amount, isSouthHemisphere=True,
+                      use_cache=False, verbose=False):
+    args = (boundary_polygon,
+            zone, ll_lat,
+            ll_long, grid_spacing,
+            lat_amount, long_amount, isSouthHemisphere)
+    kwargs = {}
+
+    if use_cache is True:
+        try:
+            from anuga.caching import cache
+        except:
+            msg = 'Caching was requested, but caching module' \
+                  'could not be imported'
+            raise msg
+
+        geo = cache(_calculate_boundary_points,
+                    args, kwargs,
+                    verbose=verbose,
+                    compression=False)
+    else:
+        geo = apply(_calculate_boundary_points, args, kwargs)
+
+    return geo
+
+
+##
+# @brief 
+# @param boundary_polygon An iterable of points that describe a polygon.
+# @param zone
+# @param ll_lat Lower left latitude, in decimal degrees
+# @param ll_long Lower left longitude, in decimal degrees
+# @param grid_spacing Grid spacing in decimal degrees.
+# @param lat_amount
+# @param long_amount
+# @param isSouthHemisphere
+def _calculate_boundary_points(boundary_polygon,
+                       zone, ll_lat,
+                       ll_long, grid_spacing,
+                       lat_amount, long_amount,
+                       isSouthHemisphere):
+    """
+    boundary_polygon - a list of points that describes a polygon.
+                      The last point is assumed ot join the first point.
+                      This is in UTM (lat long would b better though)
+
+    ll_lat - lower left latitude, in decimal degrees
+    ll-long - lower left longitude, in decimal degrees
+    grid_spacing - in decimal degrees
+
+    """
+
+    msg = "grid_spacing can not be zero"
+    assert not grid_spacing == 0, msg
+
+    a = boundary_polygon
+
+    # List of segments.  Each segment is two points.
+    segs = [i and [a[i-1], a[i]] or [a[len(a)-1], a[0]] for i in range(len(a))]
+
+    # convert the segs to Lat's and longs.
+    # Don't assume the zone of the segments is the same as the lower left
+    # corner of the lat long data!!  They can easily be in different zones
+    lat_long_set = frozenset()
+    for seg in segs:
+        points_lat_long = points_needed(seg, ll_lat, ll_long, grid_spacing,
+                                        lat_amount, long_amount, zone,
+                                        isSouthHemisphere)
+        lat_long_set |= frozenset(points_lat_long)
+
+    if lat_long_set == frozenset([]):
+        msg = "URS region specified and polygon does not overlap."
+        raise ValueError, msg
+
+    # Warning there is no info in geospatial saying the hemisphere of
+    # these points.  There should be.
+    geo = Geospatial_data(data_points=list(lat_long_set),
+                          points_are_lats_longs=True)
+
+    return geo
+
+
+##
+# @brief Get the points that are needed to interpolate any point a a segment.
+# @param seg Two points in the UTM.
+# @param ll_lat Lower left latitude, in decimal degrees
+# @param ll_long Lower left longitude, in decimal degrees
+# @param grid_spacing 
+# @param lat_amount 
+# @param long_amount 
+# @param zone 
+# @param isSouthHemisphere 
+# @return A list of points.
+def points_needed(seg, ll_lat, ll_long, grid_spacing,
+                  lat_amount, long_amount, zone,
+                  isSouthHemisphere):
+    """
+    seg is two points, in UTM
+    return a list of the points, in lats and longs that are needed to
+    interpolate any point on the segment.
+    """
+
+    from math import sqrt
+
+    geo_reference = Geo_reference(zone=zone)
+    geo = Geospatial_data(seg, geo_reference=geo_reference)
+    seg_lat_long = geo.get_data_points(as_lat_long=True,
+                                       isSouthHemisphere=isSouthHemisphere)
+
+    # 1.415 = 2^0.5, rounded up....
+    sqrt_2_rounded_up = 1.415
+    buffer = sqrt_2_rounded_up * grid_spacing
+
+    max_lat = max(seg_lat_long[0][0], seg_lat_long[1][0]) + buffer
+    max_long = max(seg_lat_long[0][1], seg_lat_long[1][1]) + buffer
+    min_lat = min(seg_lat_long[0][0], seg_lat_long[1][0]) - buffer
+    min_long = min(seg_lat_long[0][1], seg_lat_long[1][1]) - buffer
+
+    first_row = (min_long - ll_long) / grid_spacing
+
+    # To round up
+    first_row_long = int(round(first_row + 0.5))
+
+    last_row = (max_long - ll_long) / grid_spacing # round down
+    last_row_long = int(round(last_row))
+
+    first_row = (min_lat - ll_lat) / grid_spacing
+    # To round up
+    first_row_lat = int(round(first_row + 0.5))
+
+    last_row = (max_lat - ll_lat) / grid_spacing # round down
+    last_row_lat = int(round(last_row))
+
+    max_distance = 157147.4112 * grid_spacing
+    points_lat_long = []
+
+    # Create a list of the lat long points to include.
+    for index_lat in range(first_row_lat, last_row_lat + 1):
+        for index_long in range(first_row_long, last_row_long + 1):
+            lat = ll_lat + index_lat*grid_spacing
+            long = ll_long + index_long*grid_spacing
+
+            #filter here to keep good points
+            if keep_point(lat, long, seg, max_distance):
+                points_lat_long.append((lat, long)) #must be hashable
+
+    # Now that we have these points, lets throw ones out that are too far away
+    return points_lat_long
+       
+
+
+##
+# @brief 
+# @param lat
+# @param long
+# @param seg Two points in UTM.
+# @param max_distance
+def keep_point(lat, long, seg, max_distance):
+    """
+    seg is two points, UTM
+    """
+
+    from math import sqrt
+
+    _ , x0, y0 = redfearn(lat, long)
+    x1 = seg[0][0]
+    y1 = seg[0][1]
+    x2 = seg[1][0]
+    y2 = seg[1][1]
+    x2_1 = x2-x1
+    y2_1 = y2-y1
+    try:
+        d = abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1))/sqrt( \
+            (x2_1)*(x2_1)+(y2_1)*(y2_1))
+    except ZeroDivisionError:
+        if sqrt((x2_1)*(x2_1)+(y2_1)*(y2_1)) == 0 \
+           and abs((x2_1)*(y1-y0)-(x1-x0)*(y2_1)) == 0:
+            return True
+        else:
+            return False
+
+    return d <= max_distance
+
+
+
