Index: /inundation/coordinate_transforms/geo_reference.py
===================================================================
--- /inundation/coordinate_transforms/geo_reference.py	(revision 2593)
+++ /inundation/coordinate_transforms/geo_reference.py	(revision 2594)
@@ -7,5 +7,6 @@
 
 import types, sys
-from Numeric import array,Float,ArrayType
+from Numeric import array,Float,ArrayType,reshape
+from utilities.numerical_tools import ensure_numeric
 
 DEFAULT_ZONE = sys.maxint
@@ -131,33 +132,35 @@
         If the points do not have a geo ref, assume 'absolute' values
         """
+
+        # FIXME: Use ensure_numeric
         
         is_list = False
         if type(points) == types.ListType:
             is_list = True
-            if len(points)>0 and type(points[0]) not in [types.ListType,types.TupleType]:
-                #a single point is being passed.  make it a list of lists
-                points = [points]
-        elif type(points) == ArrayType:
-            if len(points.shape) == 1:
-                points = [points]
-        
-        # FIXME: currently untested
-        # Always returns a list of lists
-        if points_geo_ref is self:
-            return points        
-            
-        #convert into array
-        points = array(points).astype(Float) 
-        #add point geo ref to points
-        if not points_geo_ref is None:
-            points[:,0] += points_geo_ref.xllcorner 
-            points[:,1] += points_geo_ref.yllcorner 
-        
-        #subtract primary geo ref from points
-        points[:,0] -= self.xllcorner 
-        points[:,1] -= self.yllcorner
+
+        points = ensure_numeric(points)
+        if len(points.shape) == 1:
+            #One point has been passed
+            msg = 'Single point must have two elements'
+            assert len(points) == 2, msg
+            points = reshape(points, (1,2))
+
+            
+        if points_geo_ref is not self:
+            #add point geo ref to points
+            if not points_geo_ref is None:
+                points[:,0] += points_geo_ref.xllcorner 
+                points[:,1] += points_geo_ref.yllcorner 
+        
+            #subtract primary geo ref from points
+            points[:,0] -= self.xllcorner 
+            points[:,1] -= self.yllcorner
+
+            
         if is_list:
             points = points.tolist()
+            
         return points
+    
     
     def get_absolute(self, points):
@@ -188,4 +191,17 @@
         return points
 
+
+    def reconcile_zones(self, other):
+
+        if self.zone == other.zone or self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE:
+            pass
+        elif self.zone == DEFAULT_ZONE:
+            self.zone = other.zone
+        elif other.zone == DEFAULT_ZONE:
+            other.zone = self.zone            
+        else:    
+            msg = 'Both geospatial_data objects must be in the same \
+            ZONE to allow reconciliation. I got zone %d and %d' %(self.zone, other.zone)
+            raise msg
     
     #def easting_northing2geo_reffed_point(self, x, y):
Index: /inundation/coordinate_transforms/test_geo_reference.py
===================================================================
--- /inundation/coordinate_transforms/test_geo_reference.py	(revision 2593)
+++ /inundation/coordinate_transforms/test_geo_reference.py	(revision 2594)
@@ -159,5 +159,8 @@
         g = Geo_reference(56,x,y)
         lofl = array([[3.0,323.0]])
-        new_lofl = g.change_points_geo_ref(lofl)
+
+        
+        #print "5 lofl before",lofl         
+        new_lofl = g.change_points_geo_ref(lofl.copy())
         #print "5 lofl",lofl 
         #print "5 new_lofl",new_lofl
@@ -165,4 +168,6 @@
         self.failUnless(type(new_lofl) == ArrayType, ' failed')
         self.failUnless(type(new_lofl) == type(lofl), ' failed')
+
+
         for point,new_point in map(None,lofl,new_lofl):
             self.failUnless(point[0]-x==new_point[0], ' failed')
@@ -174,5 +179,5 @@
         g = Geo_reference(56,x,y)
         lofl = array([355.0,3.0])
-        new_lofl = g.change_points_geo_ref(lofl)
+        new_lofl = g.change_points_geo_ref(lofl.copy())        
         #print "lofl",lofl 
         #print "new_lofl",new_lofl
@@ -223,4 +228,33 @@
      
         self.failUnless(g == new_g, 'test___cmp__ failed')   
+
+
+    def test_reconcile(self):
+        g1 = Geo_reference(56,2,5)
+        g2 = Geo_reference(50,4,5)
+        g3 = Geo_reference(50,66,6)        
+        g_default = Geo_reference()                
+
+
+        g2.reconcile_zones(g3)
+        assert g2.get_zone() == g3.get_zone()
+
+        g_default.reconcile_zones(g3)
+        assert g_default.get_zone() == g3.get_zone()
+
+        g_default = Geo_reference()                
+        g3.reconcile_zones(g_default)
+        assert g_default.get_zone() == g3.get_zone()                
+
+        try:
+            g1.reconcile_zones(g2)
+        except:
+            pass
+        else:
+            msg = 'Should have raised an exception'
+            raise msg
+        
+        
+        
         
         
Index: /inundation/geospatial_data/geospatial_data.py
===================================================================
--- /inundation/geospatial_data/geospatial_data.py	(revision 2593)
+++ /inundation/geospatial_data/geospatial_data.py	(revision 2594)
@@ -176,10 +176,10 @@
         Default: absolute is True.
         """
+
         if absolute is True:
-            xll = self.geo_reference.xllcorner
-            yll = self.geo_reference.yllcorner
-            return self.data_points + [xll, yll]
+            return self.geo_reference.get_absolute(self.data_points)
         else: 
-            return self.data_points
+            return self.data_points        
+        
     
     def get_attributes(self, attribute_name = None):
@@ -202,5 +202,84 @@
         return self.attributes[attribute_name]
 
+
     def __add__(self, other):
+        """ 
+        returns the addition of 2 geospatical objects,
+        objects are concatencated to the end of each other
+            
+        NOTE: doesn't add if objects contain different 
+        attributes  
+        """
+
+
+        # find objects zone and checks if the same
+        geo_ref1 = self.get_geo_reference()
+        zone1 = geo_ref1.get_zone()
+        
+        geo_ref2 = other.get_geo_reference()
+        zone2 = geo_ref2.get_zone()
+
+        geo_ref1.reconcile_zones(geo_ref2)
+
+
+        # sets xll and yll as the smallest from self and other
+        # FIXME (Duncan and Ole): use lower left corner derived from absolute coordinates
+        if self.geo_reference.xllcorner <= other.geo_reference.xllcorner:
+            xll = self.geo_reference.xllcorner
+        else:
+            xll = other.geo_reference.xllcorner
+
+        if self.geo_reference.yllcorner <= other.geo_reference.yllcorner:
+            yll = self.geo_reference.yllcorner
+        else:
+            yll = other.geo_reference.yllcorner
+        new_geo_ref = Geo_reference(geo_ref1.get_zone(), xll, yll)
+
+        xll = yll = 0. 
+        relative_points1 = self.get_data_points(absolute=False)
+        relative_points2 = other.get_data_points(absolute=False)
+            
+        #print 'R1', relative_points1
+        #print 'R2', relative_points2        
+
+        
+        new_relative_points1 = new_geo_ref.change_points_geo_ref(relative_points1, geo_ref1)
+        new_relative_points2 = new_geo_ref.change_points_geo_ref(relative_points2, geo_ref2)
+
+        #print 'NR1', new_relative_points1
+        #print 'NR2', new_relative_points2        
+        
+        #Now both point sets are relative to new_geo_ref and zones have been reconciled
+
+
+        # Concatenate points
+        new_points = concatenate((new_relative_points1,
+                                  new_relative_points2),
+                                 axis = 0)
+
+      
+        # Concatenate attributes
+        new_attributes = {}
+        for x in self.attributes.keys():
+            if other.attributes.has_key(x):
+
+                attrib1 = self.attributes[x]
+                attrib2 = other.attributes[x]
+                new_attributes[x] = concatenate((attrib1, attrib2))
+
+            else:
+                msg = 'Both geospatial_data objects must have the same \n'
+                msg += 'attributes to allow addition.'
+                raise msg
+
+        
+        # Instantiate new data object and return    
+        return Geospatial_data(new_points,
+                               new_attributes,
+                               new_geo_ref)
+
+
+
+    def xxx__add__(self, other):
         """ 
         returns the addition of 2 geospatical objects,
Index: /inundation/geospatial_data/test_geospatial_data.py
===================================================================
--- /inundation/geospatial_data/test_geospatial_data.py	(revision 2593)
+++ /inundation/geospatial_data/test_geospatial_data.py	(revision 2594)
@@ -234,11 +234,63 @@
 
         P = G.get_data_points(absolute=True)
+
+        P_relative = G.get_data_points(absolute=False)
+        
+        assert allclose(P_relative, P - [0.1, 2.0])
+
         assert allclose(P, concatenate( (P1,P2) ))
-
-        P_relative = G.get_data_points(absolute=False)
-        assert allclose(P_relative, P - [0.1, 2.0])                      
         assert allclose(P, [[2.0, 4.1], [4.0, 7.3],
                             [5.1, 9.1], [6.1, 6.3]])
-
+        
+
+
+        
+
+    def test_add_with_geo_absolute (self):
+        """
+        Difference in Geo_reference resolved
+        """
+        points1 = array([[2.0, 4.1], [4.0, 7.3]])
+        points2 = array([[5.1, 9.1], [6.1, 6.3]])        
+        attributes1 = [2, 4]
+        attributes2 = [5, 76]
+        geo_ref1= Geo_reference(55, 1.0, 2.0)
+        geo_ref2 = Geo_reference(55, 2.0, 3.0)
+
+        
+                                
+        G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()],
+                             attributes1, geo_ref1)
+        
+        G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()],
+                             attributes2, geo_ref2)
+
+        #Check that absolute values are as expected
+        P1 = G1.get_data_points(absolute=True)
+        assert allclose(P1, points1)
+
+        P1 = G1.get_data_points(absolute=False)
+        assert allclose(P1, points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()])        
+
+        P2 = G2.get_data_points(absolute=True)
+        assert allclose(P2, points2)
+
+        P2 = G2.get_data_points(absolute=False)
+        assert allclose(P2, points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()])                
+        
+        G = G1 + G2
+        
+        assert allclose(G.get_geo_reference().get_xllcorner(), 1.0)
+        assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
+
+        P = G.get_data_points(absolute=True)
+
+        P_relative = G.get_data_points(absolute=False)
+        
+        assert allclose(P_relative, [[1.0, 2.1], [3.0, 5.3], [4.1, 7.1], [5.1, 4.3]])
+
+        assert allclose(P, concatenate( (points1,points2) ))
+                           
+        
 
 
