Index: /inundation/ga/storm_surge/pyvolution/data_manager.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 1139)
+++ /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 1140)
@@ -218,6 +218,6 @@
 
         self.precision = Float32 #Use single precision
-        if hasattr(domain,'max_size'):
-            self.max_size = domain.max_size#file size max is 2Gig
+        if hasattr(domain, 'max_size'):
+            self.max_size = domain.max_size #file size max is 2Gig
         else:
             self.max_size = max_size
@@ -246,8 +246,9 @@
             #Reference point
             #Start time in seconds since the epoch (midnight 1/1/1970)
+            #FIXME: Use Georef
             fid.starttime = domain.starttime
             fid.xllcorner = domain.xllcorner
             fid.yllcorner = domain.yllcorner
-            fid.zone = str(domain.zone) #FIXME: ?
+            fid.zone = domain.zone
 
             # dimension definitions
@@ -1067,10 +1068,13 @@
     
     """
-    from Numeric import array, Float 
+    from Numeric import array, Float, concatenate, NewAxis, zeros 
 
     #FIXME: Should be variable
     datum = 'WGS84'
     false_easting = 500000
-    false_northing = 10000000    
+    false_northing = 10000000
+    NODATA_value = -9999
+    #
+
     
     if quantity is None:
@@ -1096,4 +1100,8 @@
     x = fid.variables['x'][:]
     y = fid.variables['y'][:]
+    volumes = fid.variables['volumes'][:]
+
+    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
+    assert len(vertex_points.shape) == 2
 
     ymin = min(y); ymax = max(y)
@@ -1113,11 +1121,12 @@
 
     #Get quantity and reduce if applicable
-    q = fid.variables[quantity]
-
+    q = fid.variables[quantity][:]
+
+    
     if len(q.shape) == 2:
-        q_reduced = array( number_of_points, Float ) 
-        
+        q_reduced = zeros( number_of_points, Float ) 
+
         for k in range(number_of_points):
-            q_reduced[k] = reduction( q[:,k] ) 
+            q_reduced[k] = reduction( q[:,k] )
 
         q = q_reduced
@@ -1133,6 +1142,6 @@
     prjid.write('Units         METERS\n')
     prjid.write('Spheroid      %s\n' %datum)
-    prjid.write('Xshift        %f\n' %false_easting)
-    prjid.write('Yshift        %f\n' %false_northing)
+    prjid.write('Xshift        %d\n' %false_easting)
+    prjid.write('Yshift        %d\n' %false_northing)
     prjid.write('Parameters\n')
     prjid.close()
@@ -1143,24 +1152,57 @@
 
     from Numeric import zeros, Float 
-    points = zeros ( (ncols*nrows, 2), Float )
-
-
-    for i in xrange(ncols):
-        xg = i*cellsize
-        for j in xrange(nrows):
-            yg = j*cellsize
-            k = i*nrows + j
+    grid_points = zeros ( (ncols*nrows, 2), Float )
+
+
+    for i in xrange(nrows):
+        yg = i*cellsize        
+        for j in xrange(ncols):
+            xg = j*cellsize
+            k = i*ncols + j
 
             #print k, xg, yg
-            points[k,0] = xg
-            points[k,1] = yg
+            grid_points[k,0] = xg
+            grid_points[k,1] = yg
 
     #Interpolate
     from least_squares import Interpolation
+    from util import inside_polygon
     
-
-
+    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
+                           precrop = False)
+
+    #Interpolate using quantity values
+    grid_values = interp.interpolate(q).flat
+
+    #Write
+    ascid = open(ascfile, 'w')
+
+    ascid.write('ncols         %d\n' %ncols)
+    ascid.write('nrows         %d\n' %nrows)
+    ascid.write('xllcorner     %d\n' %xllcorner)
+    ascid.write('yllcorner     %d\n' %yllcorner)
+    ascid.write('cellsize      %f\n' %cellsize)
+    ascid.write('NODATA_value  %d\n' %NODATA_value)                    
+
+
+    #Get bounding polygon from mesh
+    P = interp.mesh.get_boundary_polygon()
+
+    #print grid_points
+    #print grid_values
     
+    for i in range(nrows):
+        for j in range(ncols):
+            index = (nrows-i-1)*ncols+j
+            
+            if inside_polygon(grid_points[index], P):
+                ascid.write('%f ' %grid_values[index])                
+            else:
+                ascid.write('%d ' %NODATA_value)
+                
+        ascid.write('\n')
+        
     #Close
+    ascid.close()
     fid.close()
     
Index: /inundation/ga/storm_surge/pyvolution/domain.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/domain.py	(revision 1139)
+++ /inundation/ga/storm_surge/pyvolution/domain.py	(revision 1140)
@@ -16,5 +16,5 @@
     def __init__(self, coordinates, vertices, boundary = None,
                  conserved_quantities = None, other_quantities = None,
-                 tagged_elements = None, zone=None, xllcorner=0, yllcorner=0):
+                 tagged_elements = None, zone=0, xllcorner=0, yllcorner=0):
 
         Mesh.__init__(self, coordinates, vertices, boundary, tagged_elements)
Index: /inundation/ga/storm_surge/pyvolution/least_squares.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 1139)
+++ /inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 1140)
@@ -31,5 +31,5 @@
 class ShapeError(exceptions.Exception): pass
 
-from general_mesh import General_mesh
+#from general_mesh import General_mesh
 from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, ArrayType
 from mesh import Mesh
Index: /inundation/ga/storm_surge/pyvolution/test_data_manager.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_data_manager.py	(revision 1139)
+++ /inundation/ga/storm_surge/pyvolution/test_data_manager.py	(revision 1140)
@@ -583,5 +583,5 @@
 
 
-    def test_sww2asc(self):
+    def test_sww2asc_elevation(self):
         """Test that sww information can be converted correctly to asc/prj
         format readable by e.g. ArcView
@@ -593,9 +593,19 @@
 
         #Setup
-        self.domain.filename = 'datatest' + str(id(self))
+        self.domain.filename = 'datatest'
+        
+        prjfile = self.domain.filename + '.prj'
+        ascfile = self.domain.filename + '.asc'        
+        swwfile = self.domain.filename + '.sww'
+        
         self.domain.set_datadir('.')
         self.domain.format = 'sww'
         self.domain.smooth = True
         self.domain.set_quantity('elevation', lambda x,y: -x-y)
+
+        self.domain.xllcorner = 308500
+        self.domain.yllcorner = 6189000
+        self.domain.zone = 56
+        
         
         sww = get_dataobject(self.domain)
@@ -606,5 +616,5 @@
         sww.store_timestep('stage')
 
-
+        cellsize = 0.25
         #Check contents
         #Get NetCDF
@@ -612,46 +622,356 @@
         fid = NetCDFFile(sww.filename, 'r')
 
-
         # Get the variables
-        x = fid.variables['x']
-        y = fid.variables['y']
-        z = fid.variables['elevation']
-        time = fid.variables['time']
-        stage = fid.variables['stage']
-        #print x[:]
-        #print y[:]
-        #print z[:]
-        #print stage[:]        
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
 
         #Export to ascii/prj files
         sww2asc(self.domain.filename, 
                 quantity = 'elevation',                         
-                cellsize = 0.25,
-                origin = (56,0,0))
-
-
-        raise 'Under construction (Ole)'
-
-        #Check values
-	#Q = self.domain.quantities['stage']
-        #Q0 = Q.vertex_values[:,0]
-        #Q1 = Q.vertex_values[:,1]
-        #Q2 = Q.vertex_values[:,2]
-        # 
-        #A = stage[1,:]
-        #assert allclose(A[0], min(Q2[0], Q1[1]))
-        #assert allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
-        #assert allclose(A[2], Q0[3])
-        #assert allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
-        #
-        ##Center point
-        #assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\
-        #                          Q0[5], Q2[6], Q1[7]))
-
+                cellsize = cellsize)
+
+
+        #Check prj (meta data)
+        prjid = open(prjfile)
+
+        lines = prjid.readlines()
+        prjid.close()
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'projection'
+        assert L[1].strip().lower() == 'utm'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'zone'
+        assert L[1].strip().lower() == '56'
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'datum'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'zunits'
+        assert L[1].strip().lower() == 'no'                        
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'units'
+        assert L[1].strip().lower() == 'meters'                
+
+        L = lines[5].strip().split()
+        assert L[0].strip().lower() == 'spheroid'
+        assert L[1].strip().lower() == 'wgs84'
+
+        L = lines[6].strip().split()
+        assert L[0].strip().lower() == 'xshift'
+        assert L[1].strip().lower() == '500000'
+
+        L = lines[7].strip().split()
+        assert L[0].strip().lower() == 'yshift'
+        assert L[1].strip().lower() == '10000000'        
+
+        L = lines[8].strip().split()
+        assert L[0].strip().lower() == 'parameters'
+        
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()        
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '5'        
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'        
+
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()            
+            y = (4-j) * cellsize
+            for i in range(5):
+                assert allclose(float(L[i]), -i*cellsize - y)
+        
 
         fid.close()
 
         #Cleanup
-        os.remove(sww.filename)
+        os.remove(prjfile)
+        os.remove(ascfile)        
+        os.remove(swwfile)
+
+
+    def test_sww2asc_stage_reduction(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+
+        This tests the reduction of quantity stage using min
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+        self.domain.filename = 'datatest'
+        
+        prjfile = self.domain.filename + '.prj'
+        ascfile = self.domain.filename + '.asc'        
+        swwfile = self.domain.filename + '.sww'
+        
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+
+        self.domain.xllcorner = 308500
+        self.domain.yllcorner = 6189000
+        self.domain.zone = 56
+        
+        
+        sww = get_dataobject(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep('stage')
+
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep('stage')
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, 'r')
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+
+        #Export to ascii/prj files
+        sww2asc(self.domain.filename, 
+                quantity = 'stage',                         
+                cellsize = cellsize,
+                reduction = min)
+
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()        
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '5'        
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'        
+
+        
+        #Check grid values (where applicable)
+        for j in range(5):
+            if j%2 == 0:
+                L = lines[6+j].strip().split()            
+                jj = 4-j
+                for i in range(5):
+                    if i%2 == 0:
+                        index = jj/2 + i/2*3
+                        val0 = stage[0,index]
+                        val1 = stage[1,index]
+
+                        #print i, j, index, ':', L[i], val0, val1
+                        assert allclose(float(L[i]), min(val0, val1))
+
+
+        fid.close()
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)        
+        #os.remove(swwfile)
+
+
+
+
+    def test_sww2asc_missing_points(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+
+        This test includes the writing of missing values
+        """
+        
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup mesh not coinciding with rectangle.
+        #This will cause missing values to occur in gridded data
+
+
+        points = [                        [1.0, 1.0],
+                              [0.5, 0.5], [1.0, 0.5],
+                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]            
+
+        vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
+        
+        #Create shallow water domain
+        domain = Domain(points, vertices)
+        domain.default_order=2
+
+
+        #Set some field values
+        domain.set_quantity('elevation', lambda x,y: -x-y)        
+        domain.set_quantity('friction', 0.03)
+
+
+        ######################
+        # Boundary conditions
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'exterior': B} )
+
+
+        ######################
+        #Initial condition - with jumps
+
+        bed = domain.quantities['elevation'].vertex_values
+        stage = zeros(bed.shape, Float)
+
+        h = 0.3
+        for i in range(stage.shape[0]):
+            if i % 2 == 0:
+                stage[i,:] = bed[i,:] + h
+            else:
+                stage[i,:] = bed[i,:]
+
+        domain.set_quantity('stage', stage)
+        domain.distribute_to_vertices_and_edges()
+
+        domain.filename = 'datatest'
+        
+        prjfile = domain.filename + '.prj'
+        ascfile = domain.filename + '.asc'        
+        swwfile = domain.filename + '.sww'
+        
+        domain.set_datadir('.')
+        domain.format = 'sww'
+        domain.smooth = True
+
+
+        domain.xllcorner = 308500
+        domain.yllcorner = 6189000
+        domain.zone = 56
+        
+        
+        sww = get_dataobject(domain)
+        sww.store_connectivity()
+        sww.store_timestep('stage')
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(swwfile, 'r')
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+
+        #Export to ascii/prj files
+        sww2asc(domain.filename, 
+                quantity = 'elevation',                         
+                cellsize = cellsize)
+
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()        
+
+        L = lines[0].strip().split()
+        assert L[0].strip().lower() == 'ncols'
+        assert L[1].strip().lower() == '5'
+
+        L = lines[1].strip().split()
+        assert L[0].strip().lower() == 'nrows'
+        assert L[1].strip().lower() == '5'        
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert allclose(float(L[1].strip().lower()), 6189000)
+
+        L = lines[4].strip().split()
+        assert L[0].strip().lower() == 'cellsize'
+        assert allclose(float(L[1].strip().lower()), cellsize)
+
+        L = lines[5].strip().split()
+        assert L[0].strip() == 'NODATA_value'
+        assert L[1].strip().lower() == '-9999'        
+
+
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()            
+            y = (4-j) * cellsize
+            for i in range(5):
+                if i+j >= 4:
+                    assert allclose(float(L[i]), -i*cellsize - y)
+                else:
+                    #Missing values
+                    assert allclose(float(L[i]), -9999)
+
+
+
+        fid.close()
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)        
+        os.remove(swwfile)
 
 
@@ -1175,5 +1495,6 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Data_Manager,'test')
+    suite = unittest.makeSuite(Test_Data_Manager,'test')    
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2asc_stage')
     runner = unittest.TextTestRunner()
     runner.run(suite)
