Index: /inundation/ga/storm_surge/pyvolution/data_manager.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 1177)
+++ /inundation/ga/storm_surge/pyvolution/data_manager.py	(revision 1178)
@@ -50,4 +50,5 @@
 from Numeric import concatenate
 
+from coordinate_transforms.geo_reference import Geo_reference, DEFAULT_ZONE
 
 def make_filename(s):
@@ -248,7 +249,4 @@
             #FIXME: Use Georef
             fid.starttime = domain.starttime
-            fid.xllcorner = domain.xllcorner
-            fid.yllcorner = domain.yllcorner
-            fid.zone = domain.zone
 
             # dimension definitions
@@ -267,5 +265,6 @@
             fid.createVariable('y', self.precision, ('number_of_points',))
             fid.createVariable('elevation', self.precision, ('number_of_points',))
-
+            if domain.geo_reference is not None:
+                domain.geo_reference.write_NetCDF(fid)
 	    #FIXME: Backwards compatibility
 	    fid.createVariable('z', self.precision, ('number_of_points',))
@@ -1077,6 +1076,4 @@
     false_northing = 10000000
     NODATA_value = -9999
-    #
-
     
     if quantity is None:
@@ -1110,7 +1107,15 @@
     number_of_points = fid.dimensions['number_of_points']
     if origin is None:
-        xllcorner = fid.xllcorner[0]
-        yllcorner = fid.yllcorner[0]
-        zone = fid.zone[0]
+        
+        # get geo_reference
+        #sww files don't have to have a geo_ref
+        try:
+            geo_reference = Geo_reference(NetCDFObject=fid)
+        except AttributeError, e:
+            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
+        
+        xllcorner = geo_reference.get_xllcorner()
+        yllcorner = geo_reference.get_yllcorner()
+        zone = geo_reference.get_zone()
     else:
         zone = origin[0]
@@ -1872,8 +1877,5 @@
     ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
 
-    xllcorner = fid.xllcorner[0]
-    yllcorner = fid.yllcorner[0]
     starttime = fid.starttime[0]
-    zone = fid.zone
     volumes = fid.variables['volumes'][:] #Connectivity
     coordinates=transpose(asarray([x.tolist(),y.tolist()]))
@@ -1883,4 +1885,11 @@
     other_quantities = []
 
+    # get geo_reference
+    #sww files don't have to have a geo_ref
+    try:
+        geo_reference = Geo_reference(NetCDFObject=fid)
+    except AttributeError, e:
+        geo_reference = None
+        
     if verbose: print '    interpolating quantities'
     for quantity in fid.variables.keys():
@@ -1916,7 +1925,6 @@
     if not boundary is None:
         domain.boundary = boundary
-    domain.zone=zone
-    domain.xllcorner=float(xllcorner)
-    domain.yllcorner=float(yllcorner)
+    
+    domain.geo_reference = geo_reference
 
     domain.starttime=float(starttime)+float(t)
Index: /inundation/ga/storm_surge/pyvolution/domain.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/domain.py	(revision 1177)
+++ /inundation/ga/storm_surge/pyvolution/domain.py	(revision 1178)
@@ -16,7 +16,8 @@
     def __init__(self, coordinates, vertices, boundary = None,
                  conserved_quantities = None, other_quantities = None,
-                 tagged_elements = None, zone=0, xllcorner=0, yllcorner=0):
-
-        Mesh.__init__(self, coordinates, vertices, boundary, tagged_elements)
+                 tagged_elements = None, geo_reference = None):
+
+        Mesh.__init__(self, coordinates, vertices, boundary,
+                      tagged_elements, geo_reference)
 
         from Numeric import zeros, Float
@@ -73,7 +74,7 @@
         #Origin in UTM coordinates
         #FIXME: This should be set if read by a msh file
-        self.zone = zone
-        self.xllcorner = xllcorner
-        self.yllcorner = yllcorner
+        #self.zone = zone
+        #self.xllcorner = xllcorner
+        #self.yllcorner = yllcorner
 
 
Index: /inundation/ga/storm_surge/pyvolution/general_mesh.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/general_mesh.py	(revision 1177)
+++ /inundation/ga/storm_surge/pyvolution/general_mesh.py	(revision 1178)
@@ -1,3 +1,5 @@
-    
+
+from coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE
+
 class General_mesh:
     """Collection of triangular elements (purely geometric)
@@ -48,5 +50,6 @@
     """
 
-    def __init__(self, coordinates, triangles, origin = None):
+    def __init__(self, coordinates, triangles,
+                 geo_reference=None):
         """
         Build triangles from x,y coordinates (sequence of 2-tuples or
@@ -63,5 +66,8 @@
         self.coordinates = array(coordinates)
 
-        self.origin = origin
+        if geo_reference is None:
+            self.geo_reference = Geo_reference(DEFAULT_ZONE,0.0,0.0) 
+        else:
+            self.geo_reference = geo_reference
 
         #Input checks
Index: /inundation/ga/storm_surge/pyvolution/least_squares.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 1177)
+++ /inundation/ga/storm_surge/pyvolution/least_squares.py	(revision 1178)
@@ -28,5 +28,9 @@
 from sparse import Sparse, Sparse_CSR
 from cg_solve import conjugate_gradient, VectorShapeError
+
+from coordinate_transforms.geo_reference import Geo_reference
+
 import time
+
 
 try:
@@ -303,6 +307,14 @@
         #FIXME: Trying the normal mesh while testing precrop,
         #       The functionality of boundary_polygon is needed for that
+
+        #FIXME - geo ref does not have to go into mesh.
+        # Change the point co-ords to conform to the
+        # mesh co-ords early in the code
+        if mesh_origin == None:
+            geo = None
+        else:
+            geo = Geo_reference(mesh_origin[0],mesh_origin[1],mesh_origin[2])
         self.mesh = Mesh(vertex_coordinates, triangles,
-                         origin = mesh_origin)
+                         geo_reference = geo)
         #FIXME, remove if mesh checks it.
         self.mesh.check_integrity()
@@ -392,5 +404,9 @@
 
         #Shift data points to same origin as mesh (if specified)
-        mesh_origin = self.mesh.origin
+
+        #FIXME this will shift if there was no geo_ref.
+        #But all this should be removed amyhow.
+        #change coords before this point
+        mesh_origin = self.mesh.geo_reference.get_origin()
         if point_coordinates is not None:
             if data_origin is not None:
Index: /inundation/ga/storm_surge/pyvolution/mesh.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/mesh.py	(revision 1177)
+++ /inundation/ga/storm_surge/pyvolution/mesh.py	(revision 1178)
@@ -58,5 +58,5 @@
 
     def __init__(self, coordinates, triangles, boundary = None,
-                 tagged_elements = None, origin = None):
+                 tagged_elements = None, geo_reference = None):
         """
         Build triangles from x,y coordinates (sequence of 2-tuples or
@@ -67,5 +67,5 @@
         from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
 
-        General_mesh.__init__(self, coordinates, triangles, origin)
+        General_mesh.__init__(self, coordinates, triangles, geo_reference)
 
         N = self.number_of_elements
Index: /inundation/ga/storm_surge/pyvolution/pmesh2domain.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/pmesh2domain.py	(revision 1177)
+++ /inundation/ga/storm_surge/pyvolution/pmesh2domain.py	(revision 1178)
@@ -14,14 +14,16 @@
     from domain import Domain
 
-    vertex_coordinates, volumes, tag_dict, vertex_quantity_dict ,tagged_elements_dict = \
-                        pmesh_to_domain(fileName,
-                                                   setting_function=setting_function)
+    vertex_coordinates, volumes, tag_dict, vertex_quantity_dict \
+                        ,tagged_elements_dict, geo_reference = \
+           pmesh_to_domain(fileName,
+                           setting_function=setting_function)
         
-    assert issubclass(DomainClass, Domain), "DomainClass is not a subclass of Domain."
+    assert issubclass(DomainClass, Domain),"DomainClass is not a subclass of Domain."
 
 
     
     domain = DomainClass(vertex_coordinates, volumes, tag_dict,
-                         tagged_elements = tagged_elements_dict )
+                         tagged_elements = tagged_elements_dict,
+                         geo_reference = geo_reference )
 
 
@@ -35,5 +37,6 @@
     # This doesn't work on the domain instance.
     # This is still needed so -ve elevations don't cuase 'lakes'
-    # The fixme we discussed was to only create a quantity when its values are set.
+    # The fixme we discussed was to only create a quantity when its values
+    #are set.
     # I think that's the way to go still
     
@@ -70,20 +73,13 @@
     volumes = mesh_dict['triangles']
 
-    #if setting_function:
-    #    if not type(setting_function) is ListType:
-    #        setting_function = [setting_function]
-    #    for funct in setting_function:
-    #        mesh_dict = funct(mesh_dict, vertices = mesh_vertices,
-    #                        volumes = volumes)
-
-
     vertex_quantity_dict = {}
     point_atts = transpose(mesh_dict['vertex_attributes'])
     point_titles  = mesh_dict['vertex_attribute_titles']
+    geo_reference  = mesh_dict['geo_reference']
     for quantity, value_vector in map (None, point_titles, point_atts):
         vertex_quantity_dict[quantity] = value_vector
     tag_dict = pmesh_dict_to_tag_dict(mesh_dict)
     tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict)
-    return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict,tagged_elements_dict
+    return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict,tagged_elements_dict, geo_reference
 
 
Index: /inundation/ga/storm_surge/pyvolution/shallow_water.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 1177)
+++ /inundation/ga/storm_surge/pyvolution/shallow_water.py	(revision 1178)
@@ -63,12 +63,11 @@
 
     def __init__(self, coordinates, vertices, boundary = None,
-                 tagged_elements = None):
+                 tagged_elements = None, geo_reference = None):
 
         conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
         other_quantities = ['elevation', 'friction']
-
         Generic_domain.__init__(self, coordinates, vertices, boundary,
                                 conserved_quantities, other_quantities,
-                                tagged_elements)
+                                tagged_elements, geo_reference)
 
         from config import minimum_allowed_height, g
Index: /inundation/ga/storm_surge/pyvolution/test_data_manager.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_data_manager.py	(revision 1177)
+++ /inundation/ga/storm_surge/pyvolution/test_data_manager.py	(revision 1178)
@@ -10,4 +10,6 @@
 from shallow_water import *
 from config import epsilon
+
+from coordinate_transforms.geo_reference import Geo_reference
 
 class Test_Data_Manager(unittest.TestCase):
@@ -604,8 +606,5 @@
         self.domain.set_quantity('elevation', lambda x,y: -x-y)
 
-        self.domain.xllcorner = 308500
-        self.domain.yllcorner = 6189000
-        self.domain.zone = 56
-        
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
         
         sww = get_dataobject(self.domain)
@@ -746,7 +745,5 @@
         self.domain.set_quantity('elevation', lambda x,y: -x-y)
 
-        self.domain.xllcorner = 308500
-        self.domain.yllcorner = 6189000
-        self.domain.zone = 56
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
         
         
@@ -897,9 +894,5 @@
         domain.smooth = True
 
-
-        domain.xllcorner = 308500
-        domain.yllcorner = 6189000
-        domain.zone = 56
-        
+        domain.geo_reference = Geo_reference(56,308500,6189000)
         
         sww = get_dataobject(domain)
@@ -919,4 +912,9 @@
         time = fid.variables['time'][:]
 
+        try:
+            geo_reference = Geo_reference(NetCDFObject=fid)
+        except AttributeError, e:
+            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
+        
         #Export to ascii/prj files
         sww2asc(domain.filename, 
@@ -1345,6 +1343,8 @@
         yiel=0.01
 	points, vertices, boundary = rectangular(10,10)
+        
 	#Create shallow water domain
 	domain = Domain(points, vertices, boundary)
+        domain.geo_reference = Geo_reference(56,11,11)
 	domain.smooth = False
 	domain.visualise = False
@@ -1395,6 +1395,7 @@
         ###################
 
-	bits = ['xllcorner','yllcorner','vertex_coordinates']
-
+	bits = [ 'geo_reference.get_xllcorner()',
+                'geo_reference.get_yllcorner()',
+                'vertex_coordinates']
 	for quantity in ['elevation']+domain.quantities_to_be_stored:
 	    bits.append('quantities["%s"].get_integral()'%quantity)
@@ -1402,6 +1403,6 @@
 
 	for bit in bits:
-	#    print 'testing that domain.'+bit+' has been restored'
-        #    print bit
+	    #print 'testing that domain.'+bit+' has been restored'
+            #print bit
         #print 'done'
 	    assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
@@ -1452,5 +1453,7 @@
         ##################
 
-	bits = ['xllcorner','yllcorner','vertex_coordinates']
+	bits = [ 'geo_reference.get_xllcorner()',
+                'geo_reference.get_yllcorner()',
+                'vertex_coordinates']
 
 	for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
@@ -1528,6 +1531,7 @@
 	    filler = 0
 	    domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
-
-	bits = ['xllcorner','yllcorner','vertex_coordinates']
+        bits = [ 'geo_reference.get_xllcorner()',
+                'geo_reference.get_yllcorner()',
+                'vertex_coordinates']
 
 	for quantity in ['elevation']+domain.quantities_to_be_stored:
@@ -1569,5 +1573,5 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(Test_Data_Manager,'test')    
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2asc_stage')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2domain')
     runner = unittest.TextTestRunner()
     runner.run(suite)
Index: /inundation/ga/storm_surge/pyvolution/test_pmesh2domain.py
===================================================================
--- /inundation/ga/storm_surge/pyvolution/test_pmesh2domain.py	(revision 1177)
+++ /inundation/ga/storm_surge/pyvolution/test_pmesh2domain.py	(revision 1178)
@@ -22,4 +22,6 @@
 from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
      Transmissive_boundary
+
+from coordinate_transforms.geo_reference import Geo_reference
 
 
@@ -64,5 +66,10 @@
 2 2 0 0 \n \
 0 # <x> <y> ...Mesh Holes... \n \
-0 # <x> <y> <attribute>...Mesh Regions... \n")
+0 # <x> <y> <attribute>...Mesh Regions... \n \
+0 # <x> <y> <attribute>...Mesh Regions, area... \n \
+# Geo reference \n \
+56 \n \
+140 \n \
+120 \n")
          file.close()
 
@@ -112,4 +119,8 @@
          self.failUnless( len(domain.boundary)  == 4,
                           "test_pmesh2Domain Too many boundaries")
+         #FIXME change to use get_xllcorner
+         self.failUnless(domain.geo_reference.xllcorner  == 140.0,
+                          "bad geo_referece")
+         
     def old_test_tags_to_boundaries (self):
          meshDict = {}
