Index: /inundation/ga/storm_surge/pyvolution/data_manager.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 895)
+++ /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 896)
@@ -1010,5 +1010,6 @@
                minlat = None, maxlat =None,
                minlon = None, maxlon =None,
-               mint = None, maxt = None, mean_stage = 0):
+               mint = None, maxt = None, mean_stage = 0,
+               origin = (0,0,0)):
     """Convert 'Ferret' NetCDF format for wave propagation to
     sww format native to pyvolution.
@@ -1021,7 +1022,9 @@
     assumed to be given in the GDA94 datum.
 
-
     min's and max's: If omitted - full extend is used.
     To include a value min may equal it, while max must exceed it.
+
+    origin is a 3-tuple with geo referenced
+    UTM coordinates (zone, easting, northing)  
     """ 
 
@@ -1195,16 +1198,20 @@
     volumes = array(volumes)        
 
-    origin = (min(x), min(y))
-    outfile.minx = origin[0]
-    outfile.miny = origin[1]    
-
-    #print x - origin[0]
-    #print y - origin[1]
-
-    #print len(x)
-    #print number_of_points
-    
-    outfile.variables['x'][:] = x - origin[0]
-    outfile.variables['y'][:] = y - origin[1]
+    if origin == None:
+        zone = refzone        
+        xllcorner = min(x)
+        yllcorner = min(y)
+    else:
+        zone = origin[0]
+        xllcorner = origin[1]
+        yllcorner = origin[2]                
+
+    
+    outfile.xllcorner = xllcorner
+    outfile.yllcorner = yllcorner
+    outfile.zone = zone
+
+    outfile.variables['x'][:] = x - xllcorner
+    outfile.variables['y'][:] = y - yllcorner
     outfile.variables['z'][:] = 0.0
     outfile.variables['elevation'][:] = 0.0    
Index: /inundation/ga/storm_surge/pyvolution/test_shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 895)
+++ /inundation/ga/storm_surge/pyvolution/test_shallow_water.py	(revision 896)
@@ -2399,6 +2399,124 @@
         from util import mean
         domain1.reduction = mean 
-        domain1.smooth = True #FIXME: This may be a requirement for the 
-	                      #interpolation process 
+        domain1.smooth = False  #Exact result
+
+        domain1.default_order = 2
+	domain1.store = True	
+        domain1.set_datadir('.')
+        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
+	
+	#FIXME: This is extremely important!
+        #How can we test if they weren't stored?
+        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+
+		
+        #Bed-slope and friction at vertices (and interpolated elsewhere)
+        domain1.set_quantity('elevation', 0)
+        domain1.set_quantity('friction', 0)	
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain1)
+        Bd = Dirichlet_boundary([0.3,0,0])
+        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
+        #Initial condition
+        domain1.set_quantity('stage', 0)
+        domain1.check_integrity()
+
+        finaltime = 5
+        #Evolution  (full domain - large steps)
+        for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
+            pass
+            #domain1.write_time()
+
+        cv1 = domain1.quantities['stage'].centroid_values
+	
+	
+	#Create an triangle shaped domain (reusing coordinates from domain 1), 
+	#formed from the lower and right hand  boundaries and 
+	#the sw-ne diagonal 
+	#from domain 1. Call it domain2
+	
+	points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
+	           [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
+		   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
+	
+	vertices = [ [1,2,0], 
+	             [3,4,1], [2,1,4], [4,5,2], 
+	             [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]	
+
+        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
+	             (4,2):'right', (6,2):'right', (8,2):'right',
+		     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}    
+		     
+        domain2 = Domain(points, vertices, boundary)
+	
+	domain2.reduction = domain1.reduction
+        domain2.smooth = False
+        domain2.default_order = 2
+	
+        #Bed-slope and friction at vertices (and interpolated elsewhere)
+        domain2.set_quantity('elevation', 0)
+        domain2.set_quantity('friction', 0)
+        domain2.set_quantity('stage', 0)        
+
+        # Boundary conditions
+        Br = Reflective_boundary(domain2)        
+	Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
+				      domain2)
+        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
+        domain2.check_integrity()        
+	
+
+
+        #Evolution (small steps)
+        for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
+            pass
+
+	
+	#Use output from domain1 as spatio-temporal boundary for domain2
+        #and verify that results at right hand side are close.	
+
+        cv2 = domain2.quantities['stage'].centroid_values
+
+        #print take(cv1, (12,14,16))  #Right
+	#print take(cv2, (4,6,8))
+        #print take(cv1, (0,6,12))  #Bottom	
+	#print take(cv2, (0,1,4))
+        #print take(cv1, (0,8,16))  #Diag	
+	#print take(cv2, (0,3,8))
+
+	assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
+	assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
+        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
+
+        #Cleanup
+        os.remove(domain1.filename + '.' + domain1.format)        
+	
+	
+
+    def test_spatio_temporal_boundary_2(self):
+        """Test that boundary values can be read from file and interpolated
+        in both time and space.
+	This is a more basic test, verifying that boundary object 
+	produces the expected results
+
+
+        """
+        import time
+       
+	#Create sww file of simple propagation from left to right 
+	#through rectangular domain
+	
+        from mesh_factory import rectangular
+
+        #Create basic mesh
+        points, vertices, boundary = rectangular(3, 3)
+
+        #Create shallow water domain
+        domain1 = Domain(points, vertices, boundary)
+	
+        from util import mean
+        domain1.reduction = mean 
+        domain1.smooth = True #To mimic MOST output
 
         domain1.default_order = 2
@@ -2426,12 +2544,11 @@
         finaltime = 5
         #Evolution  (full domain - large steps)
-        for t in domain1.evolve(yieldstep = 0.671, finaltime = finaltime):
+        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
             pass
             #domain1.write_time()
 
-        cv1 = domain1.quantities['stage'].centroid_values
-	
-	
-	#Create an triangle shaped domain (reusing coordinates from domain 1), 
+	
+	#Create an triangle shaped domain (coinciding with some
+        #coordinates from domain 1), 
 	#formed from the lower and right hand  boundaries and 
 	#the sw-ne diagonal 
@@ -2453,123 +2570,5 @@
 	
 	domain2.reduction = domain1.reduction
-        domain2.smooth = True
-        domain2.default_order = 2
-	
-        #Bed-slope and friction at vertices (and interpolated elsewhere)
-        domain2.set_quantity('elevation', 0)
-        domain2.set_quantity('friction', 0)
-        domain2.set_quantity('stage', 0)        
-
-        # Boundary conditions
-        Br = Reflective_boundary(domain2)        
-	Bf = Spatio_temporal_boundary(domain1.filename + '.' + domain1.format,
-				      domain2)
-        domain2.set_boundary({'right':Br, 'bottom':Br, 'diagonal':Bf})
-        domain2.check_integrity()        
-	
-
-
-        #Evolution (small steps)
-        for t in domain2.evolve(yieldstep = 0.0711, finaltime = finaltime):
-            pass
-
-	
-	#Use output from domain1 as spatio-temporal boundary for domain2
-        #and verify that results at right hand side are close.	
-
-        cv2 = domain2.quantities['stage'].centroid_values
-
-        #print take(cv1, (12,14,16))  #Right
-	#print take(cv2, (4,6,8))
-        #print take(cv1, (0,6,12))  #Bottom	
-	#print take(cv2, (0,1,4))
-        #print take(cv1, (0,8,16))  #Diag	
-	#print take(cv2, (0,3,8))
-
-	assert allclose( take(cv1, (0,8,16)), take(cv2, (0,3,8))) #Diag
-	assert allclose( take(cv1, (0,6,12)), take(cv2, (0,1,4))) #Bottom
-        assert allclose( take(cv1, (12,14,16)), take(cv2, (4,6,8))) #RHS
-
-        #Cleanup
-        os.remove(domain1.filename + '.' + domain1.format)        
-	
-	
-
-    def test_spatio_temporal_boundary_2(self):
-        """Test that boundary values can be read from file and interpolated
-        in both time and space.
-	This is a more basic test, verifying that boundary object 
-	produces the expected results
-
-
-        """
-        import time
-       
-	#Create sww file of simple propagation from left to right 
-	#through rectangular domain
-	
-        from mesh_factory import rectangular
-
-        #Create basic mesh
-        points, vertices, boundary = rectangular(3, 3)
-
-        #Create shallow water domain
-        domain1 = Domain(points, vertices, boundary)
-	
-        from util import mean
-        domain1.reduction = mean 
-        domain1.smooth = True #FIXME: This may be a requirement for the 
-	                      #interpolation process 
-
-        domain1.default_order = 2
-	domain1.store = True	
-        domain1.set_datadir('.')
-        domain1.set_name('spatio_temporal_boundary_source' + str(time.time()))
-	
-	#FIXME: This is extremely important!
-        #How can we test if they weren't stored?
-        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
-
-		
-        #Bed-slope and friction at vertices (and interpolated elsewhere)
-        domain1.set_quantity('elevation', 0)
-        domain1.set_quantity('friction', 0)	
-
-        # Boundary conditions
-        Br = Reflective_boundary(domain1)
-        Bd = Dirichlet_boundary([0.3,0,0])	
-        domain1.set_boundary({'left': Bd, 'top': Bd, 'right': Br, 'bottom': Br})
-        #Initial condition
-        domain1.set_quantity('stage', 0)
-        domain1.check_integrity()
-
-        finaltime = 5
-        #Evolution  (full domain - large steps)
-        for t in domain1.evolve(yieldstep = 1, finaltime = finaltime):
-            pass
-            #domain1.write_time()
-
-	
-	#Create an triangle shaped domain (reusing coordinates from domain 1), 
-	#formed from the lower and right hand  boundaries and 
-	#the sw-ne diagonal 
-	#from domain 1. Call it domain2
-	
-	points = [ [0,0], [1.0/3,0], [1.0/3,1.0/3],
-	           [2.0/3,0], [2.0/3,1.0/3], [2.0/3,2.0/3],
-		   [1,0], [1,1.0/3], [1,2.0/3], [1,1]]
-	
-	vertices = [ [1,2,0], 
-	             [3,4,1], [2,1,4], [4,5,2], 
-	             [6,7,3], [4,3,7], [7,8,4], [5,4,8], [8,9,5]]	
-
-        boundary = { (0,1):'bottom', (1,1):'bottom', (4,1): 'bottom',
-	             (4,2):'right', (6,2):'right', (8,2):'right',
-		     (0,0):'diagonal', (3,0):'diagonal', (8,0):'diagonal'}    
-		     
-        domain2 = Domain(points, vertices, boundary)
-	
-	domain2.reduction = domain1.reduction
-        domain2.smooth = True
+        domain2.smooth = False
         domain2.default_order = 2
 	
@@ -2594,5 +2593,6 @@
 	points = concatenate( (reshape(x, shp), reshape(y, shp)), axis=1)	
 	#The diagonal points of domain 1 are 0, 5, 10, 15
-	#print points[0], points[5], points[15]
+        
+	#print points[0], points[5], points[10], points[15]
 	assert allclose( take(points, [0,5,10,15]), 
 	                 [[0,0], [1.0/3, 1.0/3], [2.0/3, 2.0/3], [1,1]])  
Index: /inundation/ga/storm_surge/pyvolution/test_util.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 895)
+++ /inundation/ga/storm_surge/pyvolution/test_util.py	(revision 896)
@@ -358,6 +358,6 @@
         from util import mean
         domain1.reduction = mean 
-        domain1.smooth = True #NOTE: This is a strict requirement for the 
-	                      #interpolation process 
+        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
+                              # only one value.
 
         domain1.default_order = 2
