Index: /inundation/coordinate_transforms/__init__.py
===================================================================
--- /inundation/coordinate_transforms/__init__.py	(revision 2253)
+++ /inundation/coordinate_transforms/__init__.py	(revision 2253)
@@ -0,0 +1,3 @@
+from redfearn import *
+from point import *
+ 
Index: /inundation/coordinate_transforms/geo_reference.py
===================================================================
--- /inundation/coordinate_transforms/geo_reference.py	(revision 2253)
+++ /inundation/coordinate_transforms/geo_reference.py	(revision 2253)
@@ -0,0 +1,174 @@
+"""
+"""
+
+
+#FIXME: Ensure that all attributes of a georef are treated everywhere
+#and unit test
+
+import types
+from Numeric import array,Float,ArrayType
+
+DEFAULT_ZONE = 56
+TITLE = '#geo reference' + "\n" #this title is referred to in the .xya format
+
+class Geo_reference:
+    """
+    """
+
+    def __init__(self,
+                 zone = DEFAULT_ZONE,
+                 xllcorner = 0.0,
+                 yllcorner = 0.0,
+                 datum = 'wgs84',
+                 projection = 'UTM',
+                 units = 'm',
+                 false_easting = 500000,
+                 false_northing = 10000000, #Default for southern hemisphere 
+                 NetCDFObject=None,
+                 ASCIIFile=None,
+                 read_title=None):
+        """
+        input:
+        NetCDFObject - a handle to the netCDF file 
+        ASCIIFile - a handle to the text file
+        read_title - the title of the georeference text, if it was read in.
+        """
+
+        self.false_easting = false_easting
+        self.false_northing = false_northing        
+        self.datum = datum
+        self.projection = projection
+        self.zone = zone        
+        self.units = units
+        self.xllcorner = xllcorner
+        self.yllcorner = yllcorner        
+            
+        if NetCDFObject is not None:
+            self.read_NetCDF(NetCDFObject)
+  
+        if ASCIIFile is not None:
+            self.read_ASCII(ASCIIFile, read_title=read_title)
+
+    def get_xllcorner(self):
+        return self.xllcorner
+    
+    def get_yllcorner(self):
+        return self.yllcorner
+    
+    def get_zone(self):
+        return self.zone
+        
+    def write_NetCDF(self, outfile):
+        outfile.xllcorner = self.xllcorner 
+        outfile.yllcorner = self.yllcorner
+        outfile.zone = self.zone
+
+        outfile.false_easting = self.false_easting
+        outfile.false_northing = self.false_northing
+
+        outfile.datum = self.datum        
+        outfile.projection = self.projection
+        outfile.units = self.units
+
+    def read_NetCDF(self, infile):
+        self.xllcorner = infile.xllcorner[0]
+        self.yllcorner = infile.yllcorner[0] 
+        self.zone = infile.zone[0]
+
+        
+        assert (type(self.xllcorner) == types.FloatType or\
+                type(self.xllcorner) == types.IntType)
+        assert (type(self.yllcorner) == types.FloatType or\
+                type(self.yllcorner) == types.IntType)
+        assert (type(self.zone) == types.IntType)
+
+        #FIXME (Ole) Read in the rest...
+        
+        
+    def write_ASCII(self, fd):
+        fd.write(TITLE)
+        fd.write(str(self.zone) + "\n") 
+        fd.write(str(self.xllcorner) + "\n") 
+        fd.write(str(self.yllcorner) + "\n")
+
+    def read_ASCII(self, fd,read_title=None):
+        if read_title == None:
+            read_title = fd.readline() # remove the title line
+        #print "gr read_title[0:2].upper()",read_title[0:2].upper()
+        #print "gr TITLE[0:2].upper()",TITLE[0:2].upper() 
+        if read_title[0:2].upper() != TITLE[0:2].upper():
+            raise SyntaxError
+        self.zone = int(fd.readline())
+        self.xllcorner = float(fd.readline())
+        self.yllcorner = float(fd.readline())
+            
+        
+        assert (type(self.xllcorner) == types.FloatType)
+        assert (type(self.yllcorner) == types.FloatType)
+        assert (type(self.zone) == types.IntType)
+        
+        #false_easting = infile.false_easting[0]
+        #false_northing = infile.false_northing[0]
+        
+    def change_points_geo_ref(self, points,
+                              points_geo_ref=None):
+        """
+        Change the reference of a list or Numeric array of points.
+        If the points do not have a geo ref, assume 'absolute' values
+        """
+        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 pont is being passed.  make it a list of lists
+                points = [points]
+        elif type(points) == ArrayType:
+            if len(points.shape) == 1:
+                points = [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
+        if is_list:
+            points = points.tolist()
+        # return points
+        return points
+
+    
+    #def easting_northing2geo_reffed_point(self, x, y):
+    #    return [x-self.xllcorner, y - self.xllcorner]
+
+    #def easting_northing2geo_reffed_points(self, x, y):
+    #    return [x-self.xllcorner, y - self.xllcorner]
+
+    def get_origin(self):
+        return (self.zone, self.xllcorner, self.yllcorner)
+    
+    def __repr__(self):
+        return "(%i %f,%f)" % (self.zone, self.xllcorner, self.yllcorner)
+    
+    def __cmp__(self,other):
+        # FIXME (DSG) add a tolerence
+        if other is None:
+            return 1
+        cmp = 0
+        if not (self.xllcorner == self.xllcorner):
+            cmp = 1
+        if not (self.yllcorner == self.yllcorner):
+            cmp = 1
+        if not (self.zone == self.zone):
+            cmp = 1
+        return cmp
+        
+#-----------------------------------------------------------------------
+
+if __name__ == "__main__":
+    pass
Index: /inundation/coordinate_transforms/point.py
===================================================================
--- /inundation/coordinate_transforms/point.py	(revision 2253)
+++ /inundation/coordinate_transforms/point.py	(revision 2253)
@@ -0,0 +1,128 @@
+"""point.py - Represents a generic point on a sphere as a Python object.
+
+   See documentation of class Point for details.
+   Ole Nielsen, ANU 2002
+"""
+
+
+from math import cos, sin, pi
+
+def acos(c):
+  """acos -  Safe inverse cosine
+  
+     Input argument c is shrunk to admissible interval
+     to avoid case where a small rounding error causes
+     a math domain error.
+  """
+  from math import acos
+
+  if c > 1: c=1
+  if c < -1: c=-1
+
+  return acos(c)
+
+
+
+
+class Point:
+    """Definition of a generic point on the sphere.
+  
+    Defines a point in terms of latitude and longitude
+    and computes distances to other points on the sphere.
+
+    Initialise as
+      Point(lat, lon), where lat and lon are in decimal degrees (dd.dddd)
+      
+
+    Public Methods: 
+        distance_to(P)
+        bearing_to(P)
+        dist(P)
+
+    Author: Ole Nielsen, ANU 2002        
+    """
+     
+    # class constants
+    R = 6372000  # Approximate radius of Earth (m)
+    degrees2radians = pi/180.0
+    
+
+    def __init__(self, latitude, longitude):
+
+        assert(latitude >= -90 and latitude <= 90.0)
+        assert(longitude >= -180 and longitude <= 180.0)
+
+        self.latitude =  float(latitude)
+        self.longitude = float(longitude) 
+
+        lat = latitude * self.degrees2radians    # Converted to radians
+        lon = longitude * self.degrees2radians   # Converted to radians
+        self.coslat = cos(lat)
+        self.coslon = cos(lon)
+        self.sinlat = sin(lat)
+        self.sinlon = sin(lon)
+
+    def BearingTo(self,P):
+        """ Bearing (in degrees) to point P"""
+        AZ = self.AZ(P)
+        return int(round(AZ/self.degrees2radians))
+
+    
+    def DistanceTo(self,P):
+        """ Distance to point P"""           
+        GCA = self.GCA(P)
+        return self.R*GCA
+
+
+    def Dist(self,P):
+        """ Very cheap and rough approximation to distance"""
+        return max(abs(self.latitude-P.latitude),abs(self.longitude-P.longitude))
+
+
+    #--------------------------------------------------------------------------
+
+    def __repr__(self):
+        d = 2
+        lat = round(self.latitude,d)
+        lon = round(self.longitude,d)
+        return ' (' + str(lat)+ ', '+ str(lon)+')'
+
+
+    def GCA(self,P):
+        """ Compute the Creat Circle Angle (GCA) between current point and P.
+        """
+        
+        alpha = P.coslon*self.coslon + P.sinlon*self.sinlon
+        # The original formula is alpha = cos(self.lon - P.lon)
+        # but rewriting lets us make us of precomputed trigonometric values.
+        
+        x = alpha*self.coslat*P.coslat + self.sinlat*P.sinlat
+        return acos(x)
+
+    
+    def AZ(self,P):
+        """ Compute Azimuth bearing (AZ) from current point to P"""
+
+        # Compute cosine(AZ), where AZ is the azimuth angle
+        GCA = self.GCA(P)
+        c = P.sinlat - self.sinlat*cos(GCA)
+        c = c/self.coslat/sin(GCA)
+
+        AZ = acos(c)
+
+        # Reverse direction if bearing is westward, i.e. sin(self.lon - P.lon) > 0
+        # Without this correction the bearing due west, say, will be 90 degrees
+        # because the formulas work in the positive direction which is east.
+        #
+        # Precomputed trigonometric values are used to rewrite the formula:
+
+        if self.sinlon*P.coslon - self.coslon*P.sinlon > 0: AZ = 2*pi - AZ
+    
+        return AZ
+
+
+
+#-----------------------------------------------------------------------
+
+if __name__ == "__main__":
+    pass
Index: /inundation/coordinate_transforms/redfearn.py
===================================================================
--- /inundation/coordinate_transforms/redfearn.py	(revision 2253)
+++ /inundation/coordinate_transforms/redfearn.py	(revision 2253)
@@ -0,0 +1,160 @@
+"""Implementation of Redfearn's formula to compute UTM projections from latitude and longitude
+
+"""
+
+def degminsec2decimal_degrees(dd,mm,ss):
+    assert abs(mm) == mm
+    assert abs(ss) == ss    
+
+    if dd < 0:
+        sign = -1
+    else:
+        sign = 1
+    
+    return sign * (abs(dd) + mm/60. + ss/3600.)      
+    
+def decimal_degrees2degminsec(dec):
+
+    if dec < 0:
+        sign = -1
+    else:
+        sign = 1
+
+    dec = abs(dec)    
+    dd = int(dec)
+    f = dec-dd
+
+    mm = int(f*60)
+    ss = (f*60-mm)*60
+    
+    return sign*dd, mm, ss
+
+def redfearn(lat, lon, false_easting=None, false_northing=None):
+    """Compute UTM projection using Redfearn's formula
+
+    lat, lon is latitude and longitude in decimal degrees
+
+    If false easting and northing are specified they will override
+    the standard
+    """
+
+
+    from math import pi, sqrt, sin, cos, tan
+    
+
+
+    #GDA Specifications
+    a = 6378137.0                       #Semi major axis
+    inverse_flattening = 298.257222101  #1/f
+    K0 = 0.9996                         #Central scale factor    
+    zone_width = 6                      #Degrees
+
+    longitude_of_central_meridian_zone0 = -183    
+    longitude_of_western_edge_zone0 = -186
+
+    if false_easting is None:
+        false_easting = 500000
+
+    if false_northing is None:
+        if lat < 0:
+            false_northing = 10000000  #Southern hemisphere
+        else:
+            false_northing = 0         #Northern hemisphere)
+        
+    
+    #Derived constants
+    f = 1.0/inverse_flattening
+    b = a*(1-f)       #Semi minor axis
+
+    e2 = 2*f - f*f#    = f*(2-f) = (a^2-b^2/a^2   #Eccentricity
+    e = sqrt(e2)
+    e2_ = e2/(1-e2)   # = (a^2-b^2)/b^2 #Second eccentricity
+    e_ = sqrt(e2_)
+    e4 = e2*e2
+    e6 = e2*e4
+
+    #Foot point latitude
+    n = (a-b)/(a+b) #Same as e2 - why ?
+    n2 = n*n
+    n3 = n*n2
+    n4 = n2*n2
+
+    G = a*(1-n)*(1-n2)*(1+9*n2/4+225*n4/64)*pi/180
+
+
+    phi = lat*pi/180     #Convert latitude to radians
+
+    sinphi = sin(phi)   
+    sin2phi = sin(2*phi)
+    sin4phi = sin(4*phi)
+    sin6phi = sin(6*phi)
+
+    cosphi = cos(phi)
+    cosphi2 = cosphi*cosphi
+    cosphi3 = cosphi*cosphi2
+    cosphi4 = cosphi2*cosphi2
+    cosphi5 = cosphi*cosphi4    
+    cosphi6 = cosphi2*cosphi4
+    cosphi7 = cosphi*cosphi6
+    cosphi8 = cosphi4*cosphi4        
+
+    t = tan(phi)
+    t2 = t*t
+    t4 = t2*t2
+    t6 = t2*t4
+    
+    #Radius of Curvature
+    rho = a*(1-e2)/(1-e2*sinphi*sinphi)**1.5
+    nu = a/(1-e2*sinphi*sinphi)**0.5
+    psi = nu/rho
+    psi2 = psi*psi
+    psi3 = psi*psi2
+    psi4 = psi2*psi2
+
+
+
+    #Meridian distance
+
+    A0 = 1 - e2/4 - 3*e4/64 - 5*e6/256
+    A2 = 3.0/8*(e2+e4/4+15*e6/128)
+    A4 = 15.0/256*(e4+3*e6/4)
+    A6 = 35*e6/3072
+    
+    term1 = a*A0*phi
+    term2 = -a*A2*sin2phi
+    term3 = a*A4*sin4phi
+    term4 = -a*A6*sin6phi
+
+    m = term1 + term2 + term3 + term4 #OK
+
+    #Zone
+    zone = int((lon - longitude_of_western_edge_zone0)/zone_width)
+    central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
+
+    omega = (lon-central_meridian)*pi/180 #Relative longitude (radians)
+    omega2 = omega*omega
+    omega3 = omega*omega2
+    omega4 = omega2*omega2
+    omega5 = omega*omega4
+    omega6 = omega3*omega3
+    omega7 = omega*omega6
+    omega8 = omega4*omega4
+    
+    
+    #Northing
+    term1 = nu*sinphi*cosphi*omega2/2  
+    term2 = nu*sinphi*cosphi3*(4*psi2+psi-t2)*omega4/24
+    term3 = nu*sinphi*cosphi5*\
+            (8*psi4*(11-24*t2)-28*psi3*(1-6*t2)+\
+             psi2*(1-32*t2)-psi*2*t2+t4-t2)*omega6/720
+    term4 = nu*sinphi*cosphi7*(1385-3111*t2+543*t4-t6)*omega8/40320
+    northing = false_northing + K0*(m + term1 + term2 + term3 + term4)
+
+    #Easting
+    term1 = nu*omega*cosphi
+    term2 = nu*cosphi3*(psi-t2)*omega3/6
+    term3 = nu*cosphi5*(4*psi3*(1-6*t2)+psi2*(1+8*t2)-2*psi*t2+t4)*omega5/120
+    term4 = nu*cosphi7*(61-479*t2+179*t4-t6)*omega7/5040
+    easting = false_easting + K0*(term1 + term2 + term3 + term4)
+    
+    return zone, easting, northing
Index: /inundation/coordinate_transforms/test_geo_reference.py
===================================================================
--- /inundation/coordinate_transforms/test_geo_reference.py	(revision 2253)
+++ /inundation/coordinate_transforms/test_geo_reference.py	(revision 2253)
@@ -0,0 +1,218 @@
+#!/usr/bin/env python
+#
+
+import unittest
+import tempfile
+import os
+
+from geo_reference import *
+from Numeric import allclose,array
+
+class geo_referenceTestCase(unittest.TestCase):
+    def setUp(self):
+        pass
+        
+    def tearDown(self):
+        pass
+
+    
+    
+    def test_get_origin(self):
+        g = Geo_reference(56,1.9,1.9)
+        (z,x,y) = g.get_origin()
+
+        self.failUnless(z == g.get_zone(), ' failed')
+        self.failUnless(x == g.get_xllcorner(), ' failed')
+        self.failUnless(y == g.get_yllcorner(), ' failed') 
+        
+    def test_read_write_NetCDF(self):
+        from Scientific.IO.NetCDF import NetCDFFile
+        g = Geo_reference(56,1.9,1.9)
+        file_name = tempfile.mktemp(".geo_referenceTest")
+        
+        out_file = NetCDFFile(file_name, 'w')
+        g.write_NetCDF(out_file)
+        out_file.close()
+        
+        in_file = NetCDFFile(file_name, 'r')
+        new_g = Geo_reference(NetCDFObject=in_file)
+        in_file.close()
+        os.remove(file_name)
+
+        self.failUnless(g == new_g, 'test_read_write_NetCDF failed')  
+        
+    def test_read_write_ASCII(self):
+        from Scientific.IO.NetCDF import NetCDFFile
+        g = Geo_reference(56,1.9,1.9)
+        file_name = tempfile.mktemp(".geo_referenceTest")
+        fd = open(file_name,'w')
+        g.write_ASCII(fd)
+        fd.close()
+        
+        fd = open(file_name,'r')
+        new_g = Geo_reference(ASCIIFile=fd)
+        fd.close()
+        os.remove(file_name)
+
+        self.failUnless(g == new_g, 'test_read_write_ASCII failed')  
+    
+    def test_read_write_ASCII2(self):
+        from Scientific.IO.NetCDF import NetCDFFile
+        g = Geo_reference(56,1.9,1.9)
+        file_name = tempfile.mktemp(".geo_referenceTest")
+        fd = open(file_name,'w')
+        g.write_ASCII(fd)
+        fd.close()       
+        fd = open(file_name,'r')
+        line = fd.readline()
+        new_g = Geo_reference(ASCIIFile=fd, read_title=line)
+        fd.close()
+        os.remove(file_name)
+
+        self.failUnless(g == new_g, 'test_read_write_ASCII failed')
+        
+    def test_read_write_ASCII2(self):
+        from Scientific.IO.NetCDF import NetCDFFile
+        g = Geo_reference(56,1.9,1.9)
+        file_name = tempfile.mktemp(".geo_referenceTest")
+        fd = open(file_name,'w')
+        g.write_ASCII(fd)
+        fd.close()       
+        fd = open(file_name,'r')
+        line = fd.readline()
+        line = "fail !!"
+        try:
+            new_g = Geo_reference(ASCIIFile=fd, read_title=line)
+            fd.close()
+            os.remove(file_name)
+        except SyntaxError:
+            fd.close()
+            os.remove(file_name)
+        else:
+            self.failUnless(0 ==1,
+                        'bad text file did not raise error!')
+            
+    def test_change_points_geo_ref(self):
+        x = 3.0
+        y = 3.0
+        g = Geo_reference(56,x,y)
+        lofl = [[3.0,3.0], [6.0,6.0]]
+        new_lofl = g.change_points_geo_ref(lofl)
+        #print "lofl",lofl 
+        #print "new_lofl",new_lofl
+
+        self.failUnless(type(new_lofl) == types.ListType, ' 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')
+            self.failUnless(point[1]-y==new_point[1], ' failed')
+        
+    def test_change_points_geo_ref2(self):
+        x = 3.0
+        y = 3.0
+        g = Geo_reference(56,x,y)
+        lofl = [[3.0,3.0]]
+        new_lofl = g.change_points_geo_ref(lofl)
+        #print "lofl",lofl 
+        #print "new_lofl",new_lofl
+
+        self.failUnless(type(new_lofl) == types.ListType, ' 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')
+            self.failUnless(point[1]-y==new_point[1], ' failed')
+        
+    def test_change_points_geo_ref3(self):
+        x = 3.0
+        y = 3.0
+        g = Geo_reference(56,x,y)
+        lofl = [3.0,3.0]
+        new_lofl = g.change_points_geo_ref(lofl)
+        #print "lofl",lofl 
+        #print "new_lofl",new_lofl
+
+        self.failUnless(type(new_lofl) == types.ListType, ' 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')
+            self.failUnless(point[1]-y==new_point[1], ' failed')
+        
+    
+    def test_change_points_geo_ref4(self):
+        x = 3.0
+        y = 3.0
+        g = Geo_reference(56,x,y)
+        lofl = array([[3.0,3.0], [6.0,6.0]])
+        new_lofl = g.change_points_geo_ref(lofl)
+        #print "4 lofl",lofl 
+        #print "4 new_lofl",new_lofl
+
+        self.failUnless(type(new_lofl) == ArrayType, ' failed')
+        self.failUnless(type(new_lofl) == type(lofl), ' failed')
+        lofl[:,0] -= x
+        lofl[:,1] -= y
+        assert allclose(lofl,new_lofl)
+        
+    def test_change_points_geo_ref5(self):
+        x = 3.0
+        y = 3.0
+        g = Geo_reference(56,x,y)
+        lofl = array([[3.0,3.0]])
+        new_lofl = g.change_points_geo_ref(lofl)
+        #print "5 lofl",lofl 
+        #print "5 new_lofl",new_lofl
+
+        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')
+            self.failUnless(point[1]-y==new_point[1], ' failed')
+        
+    def test_change_points_geo_ref6(self):
+        x = 3.0
+        y = 3.0
+        g = Geo_reference(56,x,y)
+        lofl = array([3.0,3.0])
+        new_lofl = g.change_points_geo_ref(lofl)
+        #print "lofl",lofl 
+        #print "new_lofl",new_lofl
+
+        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')
+            self.failUnless(point[1]-y==new_point[1], ' failed')
+     
+    def test_change_points_geo_ref7(self):
+        x = 3.0
+        y = 3.0
+        point_x = 9.0
+        point_y = -60.0
+        g = Geo_reference(56,x,y)
+        points_geo_ref = Geo_reference(56,point_x,point_y)
+        lofl = [[3.0,3.0], [6.0,6.0]]
+        new_lofl = g.change_points_geo_ref(lofl,points_geo_ref=points_geo_ref)
+        #print "lofl",lofl 
+        #print "new_lofl",new_lofl
+
+        self.failUnless(type(new_lofl) == types.ListType, ' failed')
+        self.failUnless(type(new_lofl) == type(lofl), ' failed')
+        for point,new_point in map(None,lofl,new_lofl):
+            self.failUnless(point[0]+point_x-x==new_point[0], ' failed')
+            self.failUnless(point[1]+point_y-y==new_point[1], ' failed')
+               
+    def test___cmp__(self):
+        g = Geo_reference(56,1.9,1.9,)
+        new_g = Geo_reference(56,1.9,1.9)
+     
+        self.failUnless(g == new_g, 'test___cmp__ failed')   
+        
+        
+        
+#-------------------------------------------------------------
+if __name__ == "__main__":
+
+    suite = unittest.makeSuite(geo_referenceTestCase,'test')
+    runner = unittest.TextTestRunner() #verbosity=2)
+    runner.run(suite)
+    
Index: /inundation/coordinate_transforms/test_point.py
===================================================================
--- /inundation/coordinate_transforms/test_point.py	(revision 2253)
+++ /inundation/coordinate_transforms/test_point.py	(revision 2253)
@@ -0,0 +1,134 @@
+
+import unittest
+from point import Point
+from math import fabs
+
+#-------------------------------------------------------------
+
+class TestCase(unittest.TestCase):
+
+    def setUp(self):
+        
+        self.eps = 0.001    # Accept 0.1 % relative error
+        
+        self.RSISE = Point(-35.27456,149.12065)
+        self.Home = Point(-35.25629,149.12494)     # 28 Scrivener Street, ACT
+        self.Syd = Point(-33.93479,151.16794)      # Sydney Airport
+        self.Nadi = Point(-17.75330,177.45148)     # Nadi Airport 
+        self.Kobenhavn = Point(55.70248, 12.58364) # Kobenhavn, Denmark
+
+
+    def testBearingNorth(self):
+
+        eps = 1.0e-12
+        
+        p1 = Point(0.0,0.0)
+        p2 = Point(1.0,0.0)
+        p3 = Point(0.0,1.0)
+
+        # Assert that bearing is correct within double precision
+        self.failUnless((fabs(p1.BearingTo(p2)-0) < eps),\
+                        'Computed northward bearing: %d, Should have been: %d'\
+                         %(p1.BearingTo(p2), 0))
+
+
+    def testBearingSouth(self):
+
+        eps = 1.0e-12
+        
+        p1 = Point(0.0,0.0)
+        p2 = Point(1.0,0.0)
+        p3 = Point(0.0,1.0)
+
+        # Assert that bearing is correct within double precision
+        self.failUnless((fabs(p2.BearingTo(p1)-180) < eps),\
+                        'Computed southhward bearing: %d, Should have been: %d'\
+                         %(p2.BearingTo(p1), 180))
+
+
+    def testBearingWest(self):
+
+        eps = 1.0e-12
+        
+        p1 = Point(0.0,0.0)
+        p2 = Point(1.0,0.0)
+        p3 = Point(0.0,1.0)
+
+        # Assert that bearing is correct within double precision
+        self.failUnless((fabs(p3.BearingTo(p1)-270) < eps),\
+                        'Computed westward bearing: %d, Should have been: %d'\
+                         %(p3.BearingTo(p1), 270))
+        
+    def testBearingEast(self):
+
+        eps = 1.0e-12
+        
+        p1 = Point(0.0,0.0)
+        p2 = Point(1.0,0.0)
+        p3 = Point(0.0,1.0)
+
+        # Assert that bearing is correct within double precision
+        self.failUnless((fabs(p1.BearingTo(p3)-90) < eps),\
+                        'Computed eastward bearing: %d, Should have been: %d'\
+                         %(p1.BearingTo(p3), 90))
+        
+
+
+    def testRSISE2Home(self):
+        D = 2068   # True Distance to Home
+        B = 11     # True Bearing to Home
+        self.failUnless((fabs(self.RSISE.DistanceTo(self.Home) - D)/D < self.eps),\
+                        'Dist to Home failed')
+        self.failUnless((self.RSISE.BearingTo(self.Home) - B == 0),\
+                        'Computed bearing to Home: %d, Should have been: %d'\
+                         %(self.RSISE.BearingTo(self.Home), B))
+
+        
+
+    def testRSISE2Sydney(self):
+        D = 239.5 * 1000   # True Distance to Sydney Airport
+        B = 52             # True Bearing to Sydney Airport        
+        self.failUnless((fabs(self.RSISE.DistanceTo(self.Syd) - D)/D < self.eps),\
+                        'Dist to Sydney failed')
+        self.failUnless((self.RSISE.BearingTo(self.Syd) - B == 0),\
+                        'Computed bearing to Sydney: %d, Should have been: %d'\
+                         %(self.RSISE.BearingTo(self.Syd), B))
+
+
+    def testRSISE2Nadi(self):
+        D = 3406.1 * 1000   # True Distance to Nadi Airport
+        B = 63              # True Bearing to Nadi Airport        
+        self.failUnless((fabs(self.RSISE.DistanceTo(self.Nadi) - D)/D < self.eps),\
+                        'Dist to Nadi failed')
+        
+        self.failUnless((self.RSISE.BearingTo(self.Nadi) - B == 0),\
+                        'Computed bearing to Nadi: %d, Should have been: %d'\
+                         %(self.RSISE.BearingTo(self.Nadi), B))
+
+
+    def testRSISE2Kobenhavn(self):        
+        D = 16025 * 1000   # True Distance to Kobenhavn
+        B = 319            # True Bearing to Kobenhavn        
+        self.failUnless((fabs(self.RSISE.DistanceTo(self.Kobenhavn) - D)/D < self.eps),\
+                        'Computed Distance to Kobenhavn: %d, Should have been: %d' \
+                        %(self.RSISE.DistanceTo(self.Kobenhavn), D))
+        self.failUnless((self.RSISE.BearingTo(self.Kobenhavn) - B == 0),\
+                        'Computed Bearing to Kobenhavn: %d, Should have been: %d' \
+                        %(self.RSISE.BearingTo(self.Kobenhavn), B))
+
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+
+    mysuite = unittest.makeSuite(TestCase,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(mysuite)
+
+
+
+
+
+
+
+
Index: /inundation/coordinate_transforms/test_redfearn.py
===================================================================
--- /inundation/coordinate_transforms/test_redfearn.py	(revision 2253)
+++ /inundation/coordinate_transforms/test_redfearn.py	(revision 2253)
@@ -0,0 +1,148 @@
+
+#Test of redfearns formula. Tests can be verified at
+#
+#http://www.cellspark.com/UTM.html
+#http://www.ga.gov.au/nmd/geodesy/datums/redfearn_geo_to_grid.jsp
+
+
+import unittest
+from redfearn import *
+from Numeric import allclose
+
+
+#-------------------------------------------------------------
+
+class TestCase(unittest.TestCase):
+
+    def test_decimal_degrees_conversion(Self):
+        lat = degminsec2decimal_degrees(-37,39,10.15610)
+        lon = degminsec2decimal_degrees(143,55,35.38390) 
+        assert allclose(lat, -37.65282114)
+        assert allclose(lon, 143.9264955)
+
+        dd,mm,ss = decimal_degrees2degminsec(-37.65282114)
+        assert dd==-37
+        assert mm==39
+        assert allclose(ss, 10.15610)
+
+        dd,mm,ss = decimal_degrees2degminsec(143.9264955)
+        assert dd==143
+        assert mm==55
+        assert allclose(ss, 35.38390) 
+
+
+    def test_UTM_1(self):
+        #latitude:  -37 39' 10.15610" 
+        #Longitude: 143 55' 35.38390" 
+        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid) 
+        #Zone:   54    
+        #Easting:  758173.797  Northing: 5828674.340 
+        #Latitude:   -37  39 ' 10.15610 ''  Longitude: 143  55 ' 35.38390 '' 
+        #Grid Convergence:  1  47 ' 19.36 ''  Point Scale: 1.00042107 
+
+        lat = degminsec2decimal_degrees(-37,39,10.15610)
+        lon = degminsec2decimal_degrees(143,55,35.38390) 
+        assert allclose(lat, -37.65282114)
+        assert allclose(lon, 143.9264955)
+
+
+        zone, easting, northing = redfearn(lat,lon)
+
+        assert zone == 54
+        assert allclose(easting, 758173.797)
+        assert allclose(northing, 5828674.340)
+
+
+    def test_UTM_2(self):
+        #TEST 2
+
+        #Latitude:  -37 57 03.7203
+        #Longitude: 144 25 29.5244
+        #Zone:   55    
+        #Easting:  273741.297  Northing: 5796489.777 
+        #Latitude:   -37  57 ' 3.72030 ''  Longitude: 144  25 ' 29.52440 '' 
+        #Grid Convergence:  -1  35 ' 3.65 ''  Point Scale: 1.00023056 
+
+        lat = degminsec2decimal_degrees(-37,57,03.7203)
+        lon = degminsec2decimal_degrees(144,25,29.5244) 
+        #print lat, lon
+
+        zone, easting, northing = redfearn(lat,lon)
+
+        assert zone == 55
+        assert allclose(easting, 273741.297)
+        assert allclose(northing, 5796489.777)
+
+        
+    def test_UTM_3(self):
+        #Test 3
+        lat = degminsec2decimal_degrees(-60,0,0)
+        lon = degminsec2decimal_degrees(130,0,0) 
+
+        zone, easting, northing = redfearn(lat,lon)
+        #print zone, easting, northing
+
+        assert zone == 52
+        assert allclose(easting, 555776.267)
+        assert allclose(northing, 3348167.264)
+
+
+    def test_UTM_4(self):
+        #Test 4 (Kobenhavn, Northern hemisphere)
+        lat = 55.70248
+        dd,mm,ss = decimal_degrees2degminsec(lat)
+
+        lon = 12.58364
+        dd,mm,ss = decimal_degrees2degminsec(lon)
+
+        zone, easting, northing = redfearn(lat,lon)
+
+
+        assert zone == 33
+        assert allclose(easting, 348157.631)
+        assert allclose(northing, 6175612.993) 
+
+
+    def test_UTM_5(self):
+        #Test 5 (Wollongong)
+
+        lat = degminsec2decimal_degrees(-34,30,0.)
+        lon = degminsec2decimal_degrees(150,55,0.) 
+        
+        zone, easting, northing = redfearn(lat,lon)
+
+        #print zone, easting, northing
+
+        assert zone == 56
+        assert allclose(easting, 308728.009)
+        assert allclose(northing, 6180432.601)
+
+
+    #def test_UTM_6(self):
+    #    """Test 6 (Don's Wollongong file's ref point)
+    #    """
+    #
+    #    lat = -34.490286785873 
+    #    lon = 150.79712139578
+    #
+    #    dd,mm,ss = decimal_degrees2degminsec(lat)
+    #    print dd,mm,ss
+    #    dd,mm,ss = decimal_degrees2degminsec(lon)        
+    #    print dd,mm,ss
+    #     
+    #    zone, easting, northing = redfearn(lat,lon)
+    #
+    #    print zone, easting, northing
+    #
+    #    assert zone == 56
+    #    #assert allclose(easting, 297717.36468927) #out by 10m
+    #    #assert allclose(northing, 6181725.1724276)
+
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+
+    mysuite = unittest.makeSuite(TestCase,'test')
+    runner = unittest.TextTestRunner()
+    runner.run(mysuite)
Index: /inundation/load_mesh/__init__.py
===================================================================
--- /inundation/load_mesh/__init__.py	(revision 2253)
+++ /inundation/load_mesh/__init__.py	(revision 2253)
@@ -0,0 +1,1 @@
+#nothing
Index: /inundation/load_mesh/loadASCII.py
===================================================================
--- /inundation/load_mesh/loadASCII.py	(revision 2253)
+++ /inundation/load_mesh/loadASCII.py	(revision 2253)
@@ -0,0 +1,1323 @@
+"""
+The format for a .xya file is:
+1st line:     [attribute names]
+other lines:  x y [attributes]
+
+for example:
+elevation, friction
+0.6, 0.7, 4.9, 0.3
+1.9, 2.8, 5, 0.3
+2.7, 2.4, 5.2, 0.3
+
+The first two columns are always implicitly assumed to be x, y coordinates.
+Use the same delimiter for the attribute names and the data
+
+An xya file can optionally end with
+#geo reference
+56
+466600.0
+8644444.0
+
+When the 1st # is the zone,
+2nd # the xllcorner and 
+3rd # the yllcorner
+The format for a Points dictionary is:
+
+  ['pointlist'] a 2 column array describing points. 1st column x, 2nd column y.
+  ['attributelist'], a dictionary of 1D arrays, representing attribute values
+  at the point.  The dictionary key is the attribute header.
+    eg
+    dic['pointlist'] = [[1.0,2.0],[3.0,5.0]]
+    dic['attributelist']['elevation'] = [[7.0,5.0]
+
+    
+    The dict format for IO with mesh files is;
+    (the triangulation)
+    vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
+    vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
+    vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
+    segments: [[v1,v2],[v3,v4],...] (lists of integers) 
+    segment_tags : [tag,tag,...] list of strings
+    triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
+    triangle_tags: [s1,s2,...] A list of list of strings (probably not neccecary.  a list of string should be ok)
+    triangle_neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
+        
+    (the outline)   
+    points: [[x1,y1],[x2,y2],...] (lists of doubles)
+    point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
+    outline_segments: [[point1,point2],[p3,p4],...] (lists of integers) 
+    outline_segment_tags : [tag1,tag2,...] list of strings
+    holes : [[x1,y1],...](List of doubles, one inside each hole region)
+    regions : [ [x1,y1],...] (List of 4 doubles)
+    region_tags : [tag1,tag2,...] (list of strings)
+    region_max_areas: [ma1,ma2,...] (A list of doubles)
+    {Convension: A -ve max area means no max area}
+
+    Points files are .xya for ascii and .pts for NetCDF
+    Mesh files are .tsh for ascii and .msh for NetCDF
+
+    Note: point att's have no titles - that's currently ok, since least
+    squares adds the point info to vertices, and they have titles
+
+    The format for a .tsh file is (FIXME update this)
+    
+    First line:  <# of vertices> <# of attributes>
+    Following lines:  <vertex #> <x> <y> [attributes]
+    One line:  <# of triangles> 
+    Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
+    One line:  <# of segments> 
+    Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
+"""
+##FIXME (DSG-DSG) Is the dict format mentioned above a list of a numeric array?
+#  Needs to be defined
+##FIXME (DSG-DSG) if the ascii file isn't .xya give a better error message.
+
+
+from string import  find, rfind
+from Numeric import array, Float, Int16, Int32, Character,reshape, concatenate, take
+from os.path import splitext
+
+from coordinate_transforms.geo_reference import Geo_reference,TITLE
+
+import exceptions
+class TitleAmountError(exceptions.Exception): pass
+
+NOMAXAREA=-999
+
+# This is used by pmesh
+def import_mesh_file(ofile):
+    """
+    read a mesh file, either .tsh or .msh
+    
+    Note: will throw an IOError if it can't load the file.
+    Catch these!
+    """
+    try:
+        if ofile[-4:]== ".tsh":
+            dict = _read_tsh_file(ofile)
+        elif ofile[-4:]== ".msh":
+            dict = _read_msh_file(ofile)
+        else:      
+            msg = 'Extension %s is unknown' %ofile[-4:]
+            raise IOError, msg
+    except (SyntaxError,IndexError, ValueError): #FIXME No test for ValueError
+            msg = 'File could not be opened'
+            raise IOError, msg  
+    return dict
+
+
+def export_mesh_file(ofile,mesh_dict):
+    """
+    write a file, ofile, with the format
+    
+    First line:  <# of vertices> <# of attributes>
+    Following lines:  <vertex #> <x> <y> [attributes]
+    One line:  <# of triangles> 
+    Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
+    One line:  <# of segments> 
+    Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
+    """
+    #FIXME DSG-anyone: automate 
+    if not mesh_dict.has_key('points'):
+        mesh_dict['points'] = []
+    if not mesh_dict.has_key('point_attributes'):
+        mesh_dict['point_attributes'] = []
+    if not mesh_dict.has_key('outline_segments'):
+        mesh_dict['outline_segments'] = []
+    if not mesh_dict.has_key('outline_segment_tags'):
+        mesh_dict['outline_segment_tags'] = []
+    if not mesh_dict.has_key('holes'):
+        mesh_dict['holes'] = []
+    if not mesh_dict.has_key('regions'):
+        mesh_dict['regions'] = []
+        
+    if not mesh_dict.has_key('region_tags'):
+        mesh_dict['region_tags'] = []
+    if not mesh_dict.has_key('region_max_areas'):
+        mesh_dict['region_max_areas'] = []
+    if not mesh_dict.has_key('vertices'):
+        mesh_dict['vertices'] = []
+    if not mesh_dict.has_key('vertex_attributes'):
+        mesh_dict['vertex_attributes'] = []
+    if not mesh_dict.has_key('vertex_attribute_titles'):
+        mesh_dict['vertex_attribute_titles'] = []
+    if not mesh_dict.has_key('segments'):
+        mesh_dict['segments'] = []
+    if not mesh_dict.has_key('segment_tags'):
+        mesh_dict['segment_tags'] = []
+    if not mesh_dict.has_key('triangles'):
+        mesh_dict['triangles'] = []
+    if not mesh_dict.has_key('triangle_tags'):
+        mesh_dict['triangle_tags'] = []
+    if not mesh_dict.has_key('triangle_neighbors'):
+        mesh_dict['triangle_neighbors'] = []
+    #print "DSG************"
+    #print "mesh_dict",mesh_dict 
+    #print "DSG************"
+    
+    if (ofile[-4:] == ".tsh"):
+        _write_tsh_file(ofile,mesh_dict)
+    elif (ofile[-4:] == ".msh"):
+        _write_msh_file(ofile, mesh_dict)
+    else:
+        msg = 'Unknown file type %s ' %ofile
+        raise IOError, msg 
+
+def _read_tsh_file(ofile):
+    """
+    Read the text file format for meshes
+    """
+    fd = open(ofile,'r')
+    dict = _read_triangulation(fd)
+    dict_mesh = _read_outline(fd)
+    for element in dict_mesh.keys():
+        dict[element] = dict_mesh[element]
+    fd.close()
+    return dict
+
+    
+def _read_triangulation(fd):
+    """
+    Read the generated triangulation, NOT the outline. 
+    """
+    delimiter = " "
+    ######### loading the point info
+    line = fd.readline()
+    #print line
+    fragments = line.split()
+    #for fragment in fragments:
+    #print fragment
+    if fragments ==[]:
+        NumOfVertices = 0
+        NumOfVertAttributes = 0
+    else:
+        NumOfVertices = fragments[0]
+        NumOfVertAttributes = fragments[1]
+    points = []
+    pointattributes = []
+    for index in range(int(NumOfVertices)):        
+        #print index
+        fragments = fd.readline().split()
+        #print fragments
+        fragments.pop(0) #pop off the index
+        
+        # pop the x y off so we're left with a list of attributes       
+        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
+        points.append(vert)
+        apointattributes  = []
+        #print fragments
+        for fragment in fragments:
+            apointattributes.append(float(fragment))
+        pointattributes.append(apointattributes)
+        
+    ######### loading the point title info
+    line = fd.readline()
+    #print "point title comments",line
+    vertTitle = []
+    for index in range(int(NumOfVertAttributes)):        
+        #print index
+        fragments = fd.readline().strip() 
+        vertTitle.append(fragments)
+        #print vertTitle
+        
+    ######### loading the triangle info
+    line = fd.readline()
+    #print "triangle comments",line
+    fragments = line.split()
+    #for fragment in fragments:
+    #    print fragment
+    NumOfTriangles = fragments[0]
+    triangles = []
+    triangleattributes = []
+    triangleneighbors = []
+    for index in range(int(NumOfTriangles)):        
+        #print index
+        line = fd.readline()
+        line.strip() # so we can get the region string
+        fragments = line.split()
+        #print "triangle info", fragments
+        fragments.pop(0) #pop off the index
+        
+        tri = [int(fragments[0]),int(fragments[1]),int(fragments[2])]
+        triangles.append(tri)
+        neighbors = [int(fragments[3]),int(fragments[4]),int(fragments[5])]
+        triangleneighbors.append(neighbors)
+        for x in range(7): # remove index [<vertex #>] [<neigbouring tri #>]
+            line = line[find(line,delimiter):] # remove index
+            line = line.lstrip()
+        stringtag = line.strip()
+        triangleattributes.append(stringtag)
+        
+    ######### loading the segment info
+    line = fd.readline()
+    #print "seg comment line",line
+    fragments = line.split()
+    #for fragment in fragments:
+    #    print fragment
+    NumOfSegments = fragments[0]
+    segments = []
+    segmenttags = []
+    for index in range(int(NumOfSegments)):        
+        #print index
+        line = fd.readline()
+        line.strip() # to get the segment string
+        fragments = line.split()
+        #print fragments
+        fragments.pop(0) #pop off the index
+        seg = [int(fragments[0]),int(fragments[1])]
+        segments.append(seg)
+        line = line[find(line,delimiter):] # remove index
+        line = line.lstrip()
+        line = line[find(line,delimiter):] # remove x
+        line = line.lstrip()
+        line = line[find(line,delimiter):] # remove y
+        stringtag = line.strip()
+        segmenttags.append(stringtag)
+
+            
+    meshDict = {}
+    meshDict['vertices'] = points
+    meshDict['vertex_attributes'] = pointattributes
+    meshDict['triangles'] = triangles
+    meshDict['triangle_tags'] = triangleattributes
+    meshDict['triangle_neighbors'] = triangleneighbors 
+    meshDict['segments'] = segments 
+    meshDict['segment_tags'] = segmenttags
+    meshDict['vertex_attribute_titles'] = vertTitle
+    
+    return meshDict
+    
+def _read_outline(fd):
+    """
+    Note, if a file has no mesh info, it can still be read - the meshdic
+    returned will be 'empty'.
+    """
+    delimiter = " " # warning: split() calls are using default whitespace
+    
+    ######### loading the point info
+    line = fd.readline()
+    #print 'read_outline - line',line
+    fragments = line.split()
+    #for fragment in fragments:
+    if fragments ==[]:
+        NumOfVertices = 0
+        NumOfVertAttributes = 0
+    else:
+        NumOfVertices = fragments[0]
+        NumOfVertAttributes = fragments[1]
+    points = []
+    pointattributes = []
+    for index in range(int(NumOfVertices)):        
+        #print index
+        fragments = fd.readline().split() 
+        #print fragments
+        fragments.pop(0) #pop off the index
+        # pop the x y off so we're left with a list of attributes
+        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
+        points.append(vert)
+        apointattributes  = []
+        #print fragments
+        for fragment in fragments:
+            apointattributes.append(float(fragment))
+        pointattributes.append(apointattributes)
+        
+
+    ######### loading the segment info
+    line = fd.readline()
+    fragments = line.split()
+    #for fragment in fragments:
+    #    print fragment
+    if fragments ==[]:
+        NumOfSegments = 0
+    else:
+        NumOfSegments = fragments[0]
+    segments = []
+    segmenttags = []
+    for index in range(int(NumOfSegments)):  
+        #print index
+        line = fd.readline()
+        fragments = line.split()
+        #print fragments
+        fragments.pop(0) #pop off the index
+        
+        seg = [int(fragments[0]),int(fragments[1])]
+        segments.append(seg)
+        line = line[find(line,delimiter):] # remove index
+        line = line.lstrip()
+        line = line[find(line,delimiter):] # remove x
+        line = line.lstrip()
+        line = line[find(line,delimiter):] # remove y
+        stringtag = line.strip()
+        segmenttags.append(stringtag)  
+
+    ######### loading the hole info
+    line = fd.readline()
+    #print line
+    fragments = line.split()
+    #for fragment in fragments:
+    #    print fragment
+    if fragments ==[]:
+        numOfHoles = 0
+    else:
+        numOfHoles = fragments[0]
+    holes = []
+    for index in range(int(numOfHoles)):        
+        #print index
+        fragments = fd.readline().split()
+        #print fragments
+        fragments.pop(0) #pop off the index
+        hole = [float(fragments[0]),float(fragments[1])]
+        holes.append(hole)
+
+    
+    ######### loading the region info
+    line = fd.readline()
+    #print line
+    fragments = line.split()
+    #for fragment in fragments:
+    #    print fragment
+    if fragments ==[]:
+        numOfRegions = 0
+    else:
+        numOfRegions = fragments[0]
+    regions = []
+    regionattributes = []
+    for index in range(int(numOfRegions)):
+        line = fd.readline()
+        fragments = line.split()
+        #print fragments
+        fragments.pop(0) #pop off the index
+        region = [float(fragments[0]),float(fragments[1])]
+        regions.append(region)
+
+        line = line[find(line,delimiter):] # remove index
+        line = line.lstrip()
+        line = line[find(line,delimiter):] # remove x
+        line = line.lstrip()
+        line = line[find(line,delimiter):] # remove y
+        stringtag = line.strip()
+        regionattributes.append(stringtag)
+    regionmaxareas = []
+    line = fd.readline()
+    for index in range(int(numOfRegions)): # Read in the Max area info
+        line = fd.readline()
+        fragments = line.split()
+        #print fragments
+        # The try is here for format compatibility
+        try:
+            fragments.pop(0) #pop off the index
+            if len(fragments) == 0: #no max area
+                regionmaxareas.append(None)
+            else:
+                regionmaxareas.append(float(fragments[0]))
+        except IndexError, e:
+            regionmaxareas.append(None)
+
+    try:
+        geo_reference = Geo_reference(ASCIIFile=fd)
+    except:
+        #geo_ref not compulsory 
+        geo_reference = None
+        
+    meshDict = {}
+    meshDict['points'] = points
+    meshDict['point_attributes'] = pointattributes
+    meshDict['outline_segments'] = segments 
+    meshDict['outline_segment_tags'] = segmenttags
+    meshDict['holes'] = holes
+    meshDict['regions'] = regions
+    meshDict['region_tags'] = regionattributes
+    meshDict['region_max_areas'] = regionmaxareas
+    meshDict['geo_reference'] = geo_reference
+    
+    #print "meshDict",meshDict 
+    return meshDict
+
+def clean_line(line,delimiter):      
+    """Remove whitespace
+    """
+    #print ">%s" %line
+    line = line.strip()
+    #print "stripped>%s" %line
+    numbers = line.split(delimiter)
+    i = len(numbers) - 1
+    while i >= 0:
+        if numbers[i] == '':
+            numbers.pop(i)
+        else:
+            numbers[i] = numbers[i].strip()
+        
+        i += -1
+    #for num in numbers:
+    #    print "num>%s<" %num
+    return numbers
+
+def _write_ASCII_triangulation(fd,
+                               gen_dict):
+    vertices = gen_dict['vertices']
+    vertices_attributes = gen_dict['vertex_attributes']
+    try:
+        vertices_attribute_titles = gen_dict['vertex_attribute_titles']    
+    except KeyError, e:
+        #FIXME is this the best way? 
+        if vertices_attributes == [] or vertices_attributes[0] == []:
+             vertices_attribute_titles = []
+        else:
+            raise KeyError, e 
+        
+    triangles = gen_dict['triangles']
+    triangles_attributes = gen_dict['triangle_tags']
+    triangle_neighbors = gen_dict['triangle_neighbors']
+        
+    segments = gen_dict['segments']
+    segment_tags = gen_dict['segment_tags']
+     
+    numVert = str(len(vertices))
+    if (numVert == "0") or len(vertices_attributes) == 0:
+        numVertAttrib = "0"
+    else:
+        numVertAttrib = str(len(vertices_attributes[0]))
+    fd.write(numVert + " " + numVertAttrib + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Triangulation Vertices..." + "\n")
+
+    #<vertex #> <x> <y> [attributes]    
+    index = 0 
+    for vert in vertices:
+        attlist = ""
+        for att in vertices_attributes[index]:
+            attlist = attlist + str(att)+" "
+        attlist.strip()
+        fd.write(str(index) + " "
+                 + str(vert[0]) + " "
+                 + str(vert[1]) + " "
+                 + attlist + "\n")
+        index += 1
+
+    # write comments for title 
+    fd.write("# attribute column titles ...Triangulation Vertex Titles..." + "\n")
+    for title in vertices_attribute_titles:
+        fd.write(title + "\n")
+        
+    #<# of triangles>
+    n = len(triangles)
+    fd.write(str(n) + " # <# of triangles>, next lines <triangle #> [<vertex #>] [<neigbouring triangle #>] [attribute of region] ...Triangulation Triangles..." + "\n")
+        
+    # <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
+    for index in range(n):
+        neighbors = ""
+        tri = triangles[index]
+        for neighbor in triangle_neighbors[index]:
+            if neighbor:
+                neighbors += str(neighbor) + " "
+            else:
+                if neighbor == 0:
+                    neighbors +=  "0 "
+                else:
+                    neighbors +=  "-1 "
+        #Warning even though a list is past, only the first value
+        #is written.  There's an assumption that the list only
+        # contains one item. This assumption is made since the
+        #dict that's being passed around is also be used to communicate
+        # with triangle, and it seems to have the option of returning
+        # more than one value for triangle attributex
+        if triangles_attributes[index] == ['']:
+            att = ""
+        else:
+            att = str(triangles_attributes[index])
+        fd.write(str(index) + " "
+                 + str(tri[0]) + " " 
+                 + str(tri[1]) + " " 
+                 + str(tri[2]) + " " 
+                 + neighbors + " "
+                 + att + "\n")
+            
+    #One line:  <# of segments> 
+    fd.write(str(len(segments)) + 
+             " # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Triangulation Segments..." + "\n")
+        
+    #Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
+    for i in range(len(segments)):
+        seg = segments[i]
+        fd.write(str(i) + " "
+                 + str(seg[0]) + " " 
+                 + str(seg[1]) + " " 
+                 + str(segment_tags[i]) + "\n")
+
+
+def _write_tsh_file(ofile,mesh_dict):
+            fd = open(ofile,'w')
+            _write_ASCII_triangulation(fd,mesh_dict)
+            _write_ASCII_outline(fd,mesh_dict)    
+            fd.close()
+
+def _write_ASCII_outline(fd,
+                     dict):
+    points = dict['points']
+    point_attributes = dict['point_attributes']
+    segments = dict['outline_segments']
+    segment_tags = dict['outline_segment_tags']
+    holes = dict['holes']
+    regions = dict['regions']
+    region_tags = dict['region_tags']
+    region_max_areas = dict['region_max_areas']
+        
+    num_points = str(len(points))
+    if (num_points == "0"):
+        num_point_atts = "0"
+    else:
+        num_point_atts = str(len(point_attributes[0]))
+    fd.write(num_points + " " + num_point_atts + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Mesh Vertices..." + "\n")
+        
+    # <x> <y> [attributes]
+    for i,point in enumerate(points):
+        attlist = ""
+        for att in point_attributes[i]:
+            attlist = attlist + str(att)+" "
+        attlist.strip()
+        fd.write(str(i) + " "
+                 +str(point[0]) + " "
+                 + str(point[1]) + " "
+                 + attlist + "\n")
+            
+    #One line:  <# of segments> 
+    fd.write(str(len(segments)) + 
+             " # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Mesh Segments..." + "\n")
+        
+    #Following lines:  <vertex #>  <vertex #> [boundary tag]
+    for i,seg in enumerate(segments):
+        fd.write(str(i) + " "
+                 + str(seg[0]) + " " 
+                 + str(seg[1]) + " " 
+                 + str(segment_tags[i]) + "\n")
+
+    #One line:  <# of holes> 
+    fd.write(str(len(holes)) + 
+             " # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes..." + "\n")    
+    # <x> <y>
+    for i,h in enumerate(holes):
+        fd.write(str(i) + " "
+                 + str(h[0]) + " "
+                 + str(h[1]) + "\n")
+        
+    #One line:  <# of regions> 
+    fd.write(str(len(regions)) + 
+             " # <# of regions>, next lines <Region #> <x> <y> <tag>...Mesh Regions..." + "\n")    
+    # <index> <x> <y> <tag>
+    #print "regions",regions 
+    for i,r in enumerate(regions):
+        #print "r",r 
+        fd.write(str(i) + " " 
+                 + str(r[0]) + " "
+                 + str(r[1])+ " "
+                 + str(region_tags[i]) + "\n")
+        
+    # <index> [<MaxArea>|""]
+        
+    #One line:  <# of regions> 
+    fd.write(str(len(regions)) + 
+             " # <# of regions>, next lines <Region #> [Max Area] ...Mesh Regions..." + "\n") 
+    for i,r in enumerate(regions):
+        area = str(region_max_areas[i])
+            
+        fd.write(str(i) + " " + area + "\n")
+
+    # geo_reference info
+    if dict.has_key('geo_reference') and not dict['geo_reference'] is None:
+        dict['geo_reference'].write_ASCII(fd)
+
+def _write_msh_file(file_name, mesh):
+    """Write .msh NetCDF file   
+
+    WARNING: This function mangles the mesh data structure
+    """
+  
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    IntType = Int32
+    
+    # 
+    #the triangulation
+    mesh['vertices'] = array(mesh['vertices']).astype(Float)
+    mesh['vertex_attributes'] = array(mesh['vertex_attributes']).astype(Float)
+    mesh['vertex_attribute_titles'] = array(mesh['vertex_attribute_titles']).astype(Character) 
+    mesh['segments'] = array(mesh['segments']).astype(IntType)
+    mesh['segment_tags'] = array(mesh['segment_tags']).astype(Character)
+    mesh['triangles'] = array(mesh['triangles']).astype(IntType)
+    mesh['triangle_tags'] = array(mesh['triangle_tags']) #.astype(Character)
+    mesh['triangle_neighbors'] = array(mesh['triangle_neighbors']).astype(IntType) 
+    
+    #the outline 
+    mesh['points'] = array(mesh['points']).astype(Float)
+    mesh['point_attributes'] = array(mesh['point_attributes']).astype(Float)
+    mesh['outline_segments'] = array(mesh['outline_segments']).astype(IntType)
+    mesh['outline_segment_tags'] = array(mesh['outline_segment_tags']).astype(Character)
+    mesh['holes'] = array(mesh['holes']).astype(Float)
+    mesh['regions'] = array(mesh['regions']).astype(Float)
+    mesh['region_tags'] = array(mesh['region_tags']).astype(Character)
+    mesh['region_max_areas'] = array(mesh['region_max_areas']).astype(Float)
+    
+    #mesh = mesh_dict2array(mesh)
+    #print "new_mesh",new_mesh
+    #print "mesh",mesh 
+    
+    # NetCDF file definition
+    try:
+        outfile = NetCDFFile(file_name, 'w')
+    except IOError:
+        msg = 'File %s could not be created' %file_name
+        raise msg
+        
+    #Create new file
+    outfile.institution = 'Geoscience Australia'
+    outfile.description = 'NetCDF format for compact and portable storage ' +\
+                      'of spatial point data'
+    
+    # dimension definitions - fixed
+    outfile.createDimension('num_of_dimensions', 2) #This is 2d data
+    outfile.createDimension('num_of_segment_ends', 2) #Segs have two points
+    outfile.createDimension('num_of_triangle_vertices', 3)
+    outfile.createDimension('num_of_triangle_faces', 3)
+    outfile.createDimension('num_of_region_max_area', 1)
+
+    # Create dimensions, variables and set the variables
+    
+    # trianglulation
+    # vertices
+    if (mesh['vertices'].shape[0] > 0):
+        outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
+        outfile.createVariable('vertices', Float, ('num_of_vertices',
+                                                   'num_of_dimensions'))
+        outfile.variables['vertices'][:] = mesh['vertices']
+        if (mesh['vertex_attributes'].shape[0] > 0 and mesh['vertex_attributes'].shape[1] > 0):
+            outfile.createDimension('num_of_vertex_attributes',
+                                    mesh['vertex_attributes'].shape[1])
+            outfile.createDimension('num_of_vertex_attribute_title_chars',
+                                    mesh['vertex_attribute_titles'].shape[1])
+            outfile.createVariable('vertex_attributes',
+                                   Float,
+                                   ('num_of_vertices',
+                                    'num_of_vertex_attributes'))   
+            outfile.createVariable('vertex_attribute_titles',
+                                   Character,
+                                   ( 'num_of_vertex_attributes',
+                                     'num_of_vertex_attribute_title_chars' ))
+            outfile.variables['vertex_attributes'][:] = \
+                                                      mesh['vertex_attributes']
+            outfile.variables['vertex_attribute_titles'][:] = \
+                                     mesh['vertex_attribute_titles']
+    # segments
+    if (mesh['segments'].shape[0] > 0):       
+        outfile.createDimension('num_of_segments',
+                                mesh['segments'].shape[0])
+        outfile.createVariable('segments',
+                               IntType,
+                               ('num_of_segments', 'num_of_segment_ends'))
+        outfile.variables['segments'][:] = mesh['segments']
+        if (mesh['segment_tags'].shape[1] > 0):
+            outfile.createDimension('num_of_segment_tag_chars',
+                                    mesh['segment_tags'].shape[1])
+            outfile.createVariable('segment_tags',
+                                   Character,
+                                   ('num_of_segments',
+                                    'num_of_segment_tag_chars'))
+            outfile.variables['segment_tags'][:] = mesh['segment_tags']
+    # triangles    
+    if (mesh['triangles'].shape[0] > 0):
+        outfile.createDimension('num_of_triangles',
+                                mesh['triangles'].shape[0])
+        outfile.createVariable('triangles',
+                               IntType,
+                               ('num_of_triangles',
+                                'num_of_triangle_vertices'))
+        outfile.createVariable('triangle_neighbors',
+                               IntType,
+                               ('num_of_triangles',
+                                'num_of_triangle_faces'))
+        outfile.variables['triangles'][:] = mesh['triangles']
+        outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
+        if (mesh['triangle_tags'].shape[1] > 0):
+            outfile.createDimension('num_of_triangle_tag_chars',
+                                    mesh['triangle_tags'].shape[1]) 
+            outfile.createVariable('triangle_tags',
+                                   Character,
+                                   ('num_of_triangles',
+                                    'num_of_triangle_tag_chars'))
+            outfile.variables['triangle_tags'][:] = mesh['triangle_tags']
+    
+
+    # outline
+    # points
+    if (mesh['points'].shape[0] > 0):
+        outfile.createDimension('num_of_points', mesh['points'].shape[0])
+        outfile.createVariable('points', Float, ('num_of_points',
+                                                 'num_of_dimensions'))
+        outfile.variables['points'][:] = mesh['points'] 
+        if (mesh['point_attributes'].shape[0] > 0  and mesh['point_attributes'].shape[1] > 0):
+            outfile.createDimension('num_of_point_attributes',
+                                    mesh['point_attributes'].shape[1])
+            outfile.createVariable('point_attributes',
+                                   Float,
+                                   ('num_of_points',
+                                    'num_of_point_attributes'))
+            outfile.variables['point_attributes'][:] = mesh['point_attributes']
+    # outline_segments 
+    if (mesh['outline_segments'].shape[0] > 0):
+        outfile.createDimension('num_of_outline_segments',
+                                mesh['outline_segments'].shape[0])
+        outfile.createVariable('outline_segments',
+                               IntType,
+                               ('num_of_outline_segments',
+                                'num_of_segment_ends'))
+        outfile.variables['outline_segments'][:] = mesh['outline_segments']
+        if (mesh['outline_segment_tags'].shape[1] > 0):
+            outfile.createDimension('num_of_outline_segment_tag_chars',
+                                    mesh['outline_segment_tags'].shape[1])
+            outfile.createVariable('outline_segment_tags',
+                                   Character,
+                                   ('num_of_outline_segments',
+                                    'num_of_outline_segment_tag_chars'))
+            outfile.variables['outline_segment_tags'][:] = mesh['outline_segment_tags']
+    # holes
+    if (mesh['holes'].shape[0] > 0):
+        outfile.createDimension('num_of_holes', mesh['holes'].shape[0])
+        outfile.createVariable('holes', Float, ('num_of_holes',
+                                                'num_of_dimensions'))
+        outfile.variables['holes'][:] = mesh['holes']
+    # regions
+    if (mesh['regions'].shape[0] > 0):
+        outfile.createDimension('num_of_regions', mesh['regions'].shape[0])
+        outfile.createVariable('regions', Float, ('num_of_regions',
+                                                  'num_of_dimensions'))
+        outfile.createVariable('region_max_areas',
+                               Float,
+                               ('num_of_regions',))
+        outfile.variables['regions'][:] = mesh['regions']
+        outfile.variables['region_max_areas'][:] = mesh['region_max_areas']
+        if (mesh['region_tags'].shape[1] > 0):
+            outfile.createDimension('num_of_region_tag_chars',
+                                    mesh['region_tags'].shape[1])
+            outfile.createVariable('region_tags',
+                                   Character,
+                                   ('num_of_regions',
+                                    'num_of_region_tag_chars'))
+            outfile.variables['region_tags'][:] = mesh['region_tags']
+
+    # geo_reference info
+    if mesh.has_key('geo_reference') and not mesh['geo_reference'] == None:
+        mesh['geo_reference'].write_NetCDF(outfile)
+                                                 
+    outfile.close()
+
+
+    
+def _read_msh_file(file_name):
+    """
+    Read in an msh file.
+
+    """
+        
+    from Scientific.IO.NetCDF import NetCDFFile        	
+           
+    #Check contents
+    #Get NetCDF
+    
+    # see if the file is there.  Throw a QUIET IO error if it isn't
+    fd = open(file_name,'r')
+    fd.close()
+
+    #throws prints to screen if file not present
+    fid = NetCDFFile(file_name, 'r') 
+
+    mesh = {}
+    # Get the variables
+    # the triangulation
+    try:
+        mesh['vertices'] = fid.variables['vertices'][:]
+    except KeyError:
+        mesh['vertices'] = array([])
+    try:
+        mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:]
+    except KeyError:
+        mesh['vertex_attributes'] = []
+        for ob in mesh['vertices']:
+            mesh['vertex_attributes'].append([])
+    mesh['vertex_attribute_titles'] = []
+    try:
+        titles = fid.variables['vertex_attribute_titles'][:]
+        for i, title in enumerate(titles):
+            mesh['vertex_attribute_titles'].append(titles[i].tostring().strip())
+    except KeyError:
+        pass  
+    try:
+        mesh['segments'] = fid.variables['segments'][:]
+    except KeyError:
+        mesh['segments'] = array([])
+    mesh['segment_tags'] =[]
+    try:
+        tags = fid.variables['segment_tags'][:]
+        for i, tag in enumerate(tags):
+            mesh['segment_tags'].append(tags[i].tostring().strip())
+    except KeyError:
+        for ob in mesh['segments']:
+            mesh['segment_tags'].append('')
+    try:
+        mesh['triangles'] = fid.variables['triangles'][:]
+        mesh['triangle_neighbors'] = fid.variables['triangle_neighbors'][:]
+    except KeyError:
+        mesh['triangles'] = array([])
+        mesh['triangle_neighbors'] = array([])
+    mesh['triangle_tags'] =[]
+    try:
+        tags = fid.variables['triangle_tags'][:]
+        for i, tag in enumerate(tags):
+            mesh['triangle_tags'].append(tags[i].tostring().strip())
+    except KeyError:
+        for ob in mesh['triangles']:
+            mesh['triangle_tags'].append('')
+    
+    #the outline 
+    try:
+        mesh['points'] = fid.variables['points'][:]
+    except KeyError:
+        mesh['points'] = []
+    try:
+        mesh['point_attributes'] = fid.variables['point_attributes'][:]
+    except KeyError:
+        mesh['point_attributes'] = []
+        for point in mesh['points']:
+            mesh['point_attributes'].append([])
+    try:
+        mesh['outline_segments'] = fid.variables['outline_segments'][:]
+    except KeyError:
+        mesh['outline_segments'] = array([])
+    mesh['outline_segment_tags'] =[]
+    try:
+        tags = fid.variables['outline_segment_tags'][:]
+        for i, tag in enumerate(tags):
+            mesh['outline_segment_tags'].append(tags[i].tostring().strip())
+    except KeyError:
+        for ob in mesh['outline_segments']:
+            mesh['outline_segment_tags'].append('')
+    try:
+        mesh['holes'] = fid.variables['holes'][:]
+    except KeyError:
+        mesh['holes'] = array([])
+    try:
+        mesh['regions'] = fid.variables['regions'][:]
+    except KeyError:
+        mesh['regions'] = array([])
+    mesh['region_tags'] =[]
+    try:
+        tags = fid.variables['region_tags'][:]
+        for i, tag in enumerate(tags):
+            mesh['region_tags'].append(tags[i].tostring().strip())
+    except KeyError:
+        for ob in mesh['regions']:
+            mesh['region_tags'].append('')
+    try:
+        mesh['region_max_areas'] = fid.variables['region_max_areas'][:]
+    except KeyError:
+        mesh['region_max_areas'] = array([])
+    #mesh[''] = fid.variables[''][:]
+      
+    try:
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        mesh['geo_reference'] = geo_reference
+    except AttributeError, e:
+        #geo_ref not compulsory 
+        mesh['geo_reference'] = None
+    
+    fid.close()
+      
+    return mesh
+           
+
+## used by alpha shapes
+    
+def export_boundary_file( file_name, points, title, delimiter = ','):
+    """
+    export a file, ofile, with the format
+    
+    First line: Title variable
+    Following lines:  [point index][delimiter][point index] 
+
+    file_name - the name of the new file
+    points - List of point index pairs [[p1, p2],[p3, p4]..] 
+    title - info to write in the first line
+    """
+    
+    fd = open(file_name,'w')
+    
+    fd.write(title+"\n")
+    #[point index][delimiter][point index] 
+    for point in points:
+        fd.write( str(point[0]) + delimiter
+                  + str(point[1]) + "\n")
+    fd.close()
+
+
+###
+#  IMPORT/EXPORT POINTS FILES
+###
+
+def export_points_file(ofile,point_dict):
+    """
+    write a points file, ofile, as a binary (.msh) or text (.tsh) file
+    """
+    #this was done for all keys in the mesh file.
+    #if not mesh_dict.has_key('points'):
+    #    mesh_dict['points'] = []
+    if (ofile[-4:] == ".xya"):
+        _write_xya_file(ofile, point_dict)
+    elif (ofile[-4:] == ".pts"):
+        _write_pts_file(ofile, point_dict)
+    else:
+        msg = 'Unknown file type %s ' %ofile
+        raise IOError, msg 
+               
+def import_points_file(ofile, delimiter = None, verbose = False):
+    """
+    load an .xya or .pts file
+
+    Note: will throw an IOError if it can't load the file.
+    Catch these!
+    """
+    
+    if ofile[-4:]== ".xya":
+        try:
+            if delimiter == None:
+                try:
+                    fd = open(ofile)
+                    points_dict = _read_xya_file(fd, ',')
+                except SyntaxError:
+                    fd.close()
+                    fd = open(ofile)
+                    points_dict = _read_xya_file(fd, ' ')
+            else:
+                fd = open(ofile)
+                points_dict = _read_xya_file(fd, delimiter)
+            fd.close()
+        except (IndexError,ValueError,SyntaxError):
+            fd.close()    
+            msg = 'Could not open file %s ' %ofile
+            raise IOError, msg
+        except IOError:
+            # Catch this to add an error message
+            msg = 'Could not open file %s ' %ofile
+            raise IOError, msg
+        
+    elif ofile[-4:]== ".pts":
+        try:
+            points_dict = _read_pts_file(ofile, verbose)        
+        except IOError, e:    
+            msg = 'Could not open file %s ' %ofile
+            raise IOError, msg
+    else:      
+        msg = 'Extension %s is unknown' %ofile[-4:]
+        raise IOError, msg
+    return points_dict
+
+def extent_point_atts(point_atts):
+    """
+    Returns 4 points representing the extent
+    This losses attribute info.
+    """
+    point_atts['pointlist'] = extent(point_atts['pointlist'])
+    point_atts['attributelist'] = {}
+    return point_atts
+    
+def extent(points):
+    points = array(points).astype(Float)
+    max_x = min_x = points[0][0]
+    max_y = min_y = points[0][1]
+    for point in points[1:]:
+        x = point[0] 
+        if x > max_x: max_x = x
+        if x < min_x: min_x = x           
+        y = point[1] 
+        if y > max_y: max_y = y
+        if y < min_y: min_y = y
+    extent = array([[min_x,min_y],[max_x,min_y],[max_x,max_y],[min_x,max_y]])
+    #print "extent",extent 
+    return extent
+    
+def reduce_pts(infile, outfile, max_points, verbose = False):
+    """
+    reduces a points file by removong every second point until the # of points
+    is less than max_points.
+    """
+    point_atts = _read_pts_file(infile)
+    while point_atts['pointlist'].shape[0] > max_points:
+        if verbose: print "point_atts['pointlist'].shape[0]"
+        point_atts = half_pts(point_atts)
+    export_points_file(outfile, point_atts)
+ 
+def produce_half_point_files(infile, max_points, delimiter, verbose = False):
+    point_atts = _read_pts_file(infile)
+    root, ext = splitext(infile)
+    outfiles = []
+    if verbose: print "# of points", point_atts['pointlist'].shape[0]
+    while point_atts['pointlist'].shape[0] > max_points:
+        point_atts = half_pts(point_atts)
+        if verbose: print "# of points", point_atts['pointlist'].shape[0]
+        outfile = root + delimiter + str(point_atts['pointlist'].shape[0]) + ext
+        outfiles.append(outfile)
+        export_points_file(outfile,
+                  point_atts)
+    return outfiles
+    
+def point_atts2array(point_atts):
+    point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
+    
+    for key in point_atts['attributelist'].keys():
+        point_atts['attributelist'][key]= array(point_atts['attributelist'][key]).astype(Float)
+    return point_atts
+    
+def half_pts(point_atts):
+    point_atts2array(point_atts)
+    point_atts['pointlist'] = point_atts['pointlist'][::2]
+    
+    for key in point_atts['attributelist'].keys():
+        point_atts['attributelist'][key]=  point_atts['attributelist'][key] [::2]
+    return point_atts
+
+def concatinate_attributelist(dic):
+    """
+    giving a dic[attribute title] = attribute
+    return list of attribute titles, array of attributes 
+    """
+    point_attributes = array([]).astype(Float)
+    keys = dic.keys()
+    key = keys.pop(0)
+    point_attributes = reshape(dic[key],(dic[key].shape[0],1))
+    for key in keys:
+        #point_attributes = concatenate([point_attributes, dic[key]], axis=1)
+        reshaped = reshape(dic[key],(dic[key].shape[0],1))
+        point_attributes = concatenate([point_attributes, reshaped], axis=1)
+    return dic.keys(), point_attributes
+
+def _read_pts_file(file_name, verbose = False):
+    """Read .pts NetCDF file
+   
+    Return a dic of array of points, and dic of array of attribute
+
+    eg
+    dic['points'] = [[1.0,2.0],[3.0,5.0]]
+    dic['attributelist']['elevation'] = [[7.0,5.0]
+    """
+    #FIXME: (DSG)This format has issues.
+    # There can't be an attribute called points 
+    # consider format change
+
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    if verbose: print 'Reading ', file_name
+
+    
+    # see if the file is there.  Throw a QUIET IO error if it isn't
+    fd = open(file_name,'r')
+    fd.close()
+
+    #throws prints to screen if file not present
+    fid = NetCDFFile(file_name, 'r') 
+
+    point_atts = {}  
+    # Get the variables
+    point_atts['pointlist'] = array(fid.variables['points'])
+    keys = fid.variables.keys()
+    if verbose: print 'Got %d variables: %s' %(len(keys), keys)
+    try:
+        keys.remove('points')
+    except IOError, e:       
+        fid.close()    
+        msg = 'Expected keyword "points" but could not find it'
+        raise IOError, msg
+
+    attributes = {}
+    for key in keys:
+        if verbose: print "reading attribute '%s'" %key
+        
+        attributes[key] = array(fid.variables[key])
+        
+    point_atts['attributelist'] = attributes
+    
+    try:
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        point_atts['geo_reference'] = geo_reference
+    except AttributeError, e:
+        #geo_ref not compulsory 
+        point_atts['geo_reference'] = None
+    
+    
+    fid.close()
+    
+    #print "point_atts",point_atts 
+    return point_atts
+
+def _write_pts_file(file_name, point_atts):
+    """Write .pts NetCDF file   
+
+    WARNING: This function mangles the point_atts data structure
+    """
+    #FIXME: (DSG)This format has issues.
+    # There can't be an attribute called points 
+    # consider format change
+    
+    from Scientific.IO.NetCDF import NetCDFFile
+    point_atts2array(point_atts)
+    # NetCDF file definition
+    outfile = NetCDFFile(file_name, 'w')
+    
+    #Create new file
+    outfile.institution = 'Geoscience Australia'
+    outfile.description = 'NetCDF format for compact and portable storage ' +\
+                      'of spatial point data'
+    
+    # dimension definitions
+    shape = point_atts['pointlist'].shape[0]
+    outfile.createDimension('number_of_points', shape)  
+    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
+    
+    # variable definition
+    outfile.createVariable('points', Float, ('number_of_points',
+                                             'number_of_dimensions'))
+
+    #create variables  
+    outfile.variables['points'][:] = point_atts['pointlist'] #.astype(Float32)
+    for key in point_atts['attributelist'].keys():
+        outfile.createVariable(key, Float, ('number_of_points',))
+        outfile.variables[key][:] = point_atts['attributelist'][key] #.astype(Float32)
+        
+    if point_atts.has_key('geo_reference') and not point_atts['geo_reference'] == None:
+        point_atts['geo_reference'].write_NetCDF(outfile)
+        
+    outfile.close() 
+   
+def _read_xya_file(fd, delimiter):
+    #lines = fd.readlines()
+    points = []
+    pointattributes = []
+    #if len(lines) <= 1:
+    #    raise SyntaxError
+    title = fd.readline()
+    #title = lines.pop(0) # the first (title) line
+    att_names = clean_line(title,delimiter)
+
+    att_dict = {}
+    line = fd.readline()
+    numbers = clean_line(line,delimiter)
+    #for line in lines:
+    while len(numbers) > 1:
+        #print "line >%s" %line
+        #numbers = clean_line(line,delimiter)
+        #print "numbers >%s<" %numbers
+        #if len(numbers) < 2 and numbers != []:
+            
+            # A line without two numbers
+            # or a bad delimiter
+            #FIXME dsg-dsg change error that is raised.
+            #raise SyntaxError
+        if numbers != []:
+            try:
+                x = float(numbers[0])
+                y = float(numbers[1])
+                points.append([x,y])
+                numbers.pop(0)
+                numbers.pop(0)
+                #attributes = []
+                #print "att_names",att_names
+                #print "numbers",numbers 
+                if len(att_names) != len(numbers):
+                    fd.close()
+                    # It might not be a problem with the title
+                    #raise TitleAmountError
+                    raise IOError
+                for i,num in enumerate(numbers):
+                    num.strip()
+                    if num != '\n' and num != '':
+                        #attributes.append(float(num))
+                        att_dict.setdefault(att_names[i],[]).append(float(num))
+                    
+            except ValueError:
+                raise SyntaxError
+        line = fd.readline()
+        numbers = clean_line(line,delimiter)
+        
+    if line == '':
+        # end of file
+        geo_reference = None
+    else:
+        geo_reference = Geo_reference(ASCIIFile=fd,read_title=line)
+    
+    xya_dict = {}
+    xya_dict['pointlist'] = array(points).astype(Float)
+    
+    for key in att_dict.keys():
+        att_dict[key] = array(att_dict[key]).astype(Float)
+    xya_dict['attributelist'] = att_dict
+    xya_dict['geo_reference'] = geo_reference
+    #print "xya_dict",xya_dict 
+    return xya_dict
+
+#FIXME(dsg), turn this dict plus methods into a class?
+def take_points(dict,indices_to_keep):
+    dict = point_atts2array(dict)
+    #FIXME maybe the points data structure should become a class?
+    dict['pointlist'] = take(dict['pointlist'],indices_to_keep)
+
+    for key in dict['attributelist'].keys():
+        dict['attributelist'][key]= take(dict['attributelist'][key],
+                                         indices_to_keep)
+    return dict
+    
+def add_point_dictionaries (dict1, dict2):
+    """
+    """
+    dict1 = point_atts2array(dict1)
+    dict2 = point_atts2array(dict2)
+    
+    combined = {}  
+    combined['pointlist'] = concatenate((dict2['pointlist'],
+                                         dict1['pointlist']),axis=0)   
+    atts = {}   
+    for key in dict2['attributelist'].keys():
+        atts[key]= concatenate((dict2['attributelist'][key],
+                                dict1['attributelist'][key]), axis=0)
+    combined['attributelist']=atts
+    combined['geo_reference'] = dict1['geo_reference']
+    return combined
+                 
+def _write_xya_file( file_name, xya_dict, delimiter = ','):
+    """
+    export a file, ofile, with the xya format
+    
+    """
+    points = xya_dict['pointlist'] 
+    pointattributes = xya_dict['attributelist']
+    
+    fd = open(file_name,'w')
+ 
+    titlelist = ""
+    for title in pointattributes.keys():
+        titlelist = titlelist + title + delimiter
+    titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
+    fd.write(titlelist+"\n")
+    #<vertex #> <x> <y> [attributes]
+    for i,vert in enumerate( points):
+        
+        attlist = ","
+        for att in pointattributes.keys():
+            attlist = attlist + str(pointattributes[att][i])+ delimiter
+        attlist = attlist[0:-len(delimiter)] # remove the last delimiter
+        attlist.strip()
+        fd.write( str(vert[0]) + delimiter
+                  + str(vert[1])
+                  + attlist + "\n")
+    
+    # geo_reference info
+    if xya_dict.has_key('geo_reference') and \
+           not xya_dict['geo_reference'] is None:
+        xya_dict['geo_reference'].write_ASCII(fd)
+    fd.close()
+     
+if __name__ == "__main__":
+    pass
Index: /inundation/load_mesh/test_loadASCII.py
===================================================================
--- /inundation/load_mesh/test_loadASCII.py	(revision 2253)
+++ /inundation/load_mesh/test_loadASCII.py	(revision 2253)
@@ -0,0 +1,1054 @@
+#!/usr/bin/env python
+#
+
+import tempfile
+import unittest
+
+import os
+import tempfile
+
+from os.path import splitext
+
+from Numeric import array, allclose 
+    
+from loadASCII import *
+from coordinate_transforms.geo_reference import Geo_reference
+import loadASCII
+
+class loadASCIITestCase(unittest.TestCase):
+    def setUp(self):
+        self.dict ={}
+        self.dict['outline_segments'] = [(0, 1), (1, 2), (0, 2), (0, 3)]
+        self.dict['outline_segment_tags'] = ['50', '40', '30', '20']
+        self.dict['holes'] = [(0.2, 0.6)]
+        self.dict['point_attributes'] = [[5, 2], [4, 2], [3, 2], [2,2]]
+        self.dict['regions'] = [(0.3, 0.3),(0.3, 0.4)]
+        self.dict['region_tags'] = ['1.3', 'yeah']
+        self.dict['region_max_areas'] = [36.0,-7.1]
+        self.dict['points'] = [(0.0, 0.0), (0.0, 4.0), (4.0, 0.0), (1.0, 1.0)]
+        self.dict['vertices'] = [(0.0, 0.0), (0.0, 4.0),
+                                 (4.0, 0.0), (1.0, 1.0), (2.0, 2.0)]
+        self.dict['triangles'] = [(3, 2, 4), (1, 0, 3),
+                                             (3, 4,1), (2, 3, 0)]
+        self.dict['segments'] = [(0, 1), (1, 4), (2, 0),
+                                            (0, 3), (4, 2)]
+        self.dict['triangle_tags'] = ['1.3', '1.3',
+                                      '1.3', '1.3']
+        self.dict['vertex_attributes'] = [[1.2,2.], [1.2,2.],
+                                          [1.2,2.], [1.2,2.], [1.2,3.]]
+        self.dict['triangle_neighbors'] = [[-1, 2, 3], [3, 2, -1],
+                                           [-1, 1, 0], [1, -1, 0]]
+        self.dict['segment_tags'] = ['50', '40', '30', '20', '40']
+        self.dict['vertex_attribute_titles'] = ['bed elevation', 'height']
+        self.dict['geo_reference'] = Geo_reference(56,1.9,1.9)
+        
+        self.sparse_dict ={}
+        self.sparse_dict['outline_segments'] = []
+        self.sparse_dict['outline_segment_tags'] = []
+        self.sparse_dict['holes'] = []
+        self.sparse_dict['points'] = [(0.0, 0.0),(9,8)]
+        self.sparse_dict['point_attributes'] = [[],[]] # points don't have to have
+                                                    # attributes
+        self.sparse_dict['regions'] = []
+        self.sparse_dict['region_tags'] = []
+        self.sparse_dict['region_max_areas'] = []
+
+        self.sparse_dict['vertices'] = []
+        self.sparse_dict['triangles'] = []
+        self.sparse_dict['segments'] = []
+        self.sparse_dict['triangle_tags'] = []
+        self.sparse_dict['vertex_attributes'] = []
+        self.sparse_dict['triangle_neighbors'] = []
+        self.sparse_dict['segment_tags'] = []
+        self.sparse_dict['vertex_attribute_titles'] = []
+        
+        self.blank_dict ={}
+        self.blank_dict['outline_segments'] = []
+        self.blank_dict['outline_segment_tags'] = []
+        self.blank_dict['holes'] = []
+        self.blank_dict['points'] = []
+        self.blank_dict['point_attributes'] = [] 
+        self.blank_dict['regions'] = []
+        self.blank_dict['region_tags'] = []
+        self.blank_dict['region_max_areas'] = []
+        self.blank_dict['vertices'] = []
+        self.blank_dict['triangles'] = []
+        self.blank_dict['segments'] = []
+        self.blank_dict['triangle_tags'] = []
+        self.blank_dict['vertex_attributes'] = []
+        self.blank_dict['triangle_neighbors'] = []
+        self.blank_dict['segment_tags'] = []
+        self.blank_dict['vertex_attribute_titles'] = []
+        
+        self.tri_dict ={}
+        self.tri_dict['outline_segments'] = [[0,1]]
+        self.tri_dict['outline_segment_tags'] = ['']
+        self.tri_dict['holes'] = []
+        self.tri_dict['points'] = [(9,8),(7,8)]
+        self.tri_dict['point_attributes'] = [[],[]] 
+        self.tri_dict['regions'] = []
+        self.tri_dict['region_tags'] = []
+        self.tri_dict['region_max_areas'] = []
+        self.tri_dict['vertices'] = [[9,8],[7,8], [4,5]]
+        self.tri_dict['triangles'] = [[0,1,2]]
+        self.tri_dict['segments'] = [[0,1]]
+        self.tri_dict['triangle_tags'] = ['']
+        self.tri_dict['vertex_attributes'] = [[],[],[]]
+        self.tri_dict['triangle_neighbors'] = [[0,0,0]]
+        self.tri_dict['segment_tags'] = ['']
+        self.tri_dict['vertex_attribute_titles'] = []
+        
+        self.seg_dict ={}
+        self.seg_dict['outline_segments'] = [[0,1]]
+        self.seg_dict['outline_segment_tags'] = ['']
+        self.seg_dict['holes'] = []
+        self.seg_dict['points'] = [(9,8),(7,8)]
+        self.seg_dict['point_attributes'] = [[],[]]  
+        self.seg_dict['regions'] = [(5,4)]
+        self.seg_dict['region_tags'] = [''] 
+        self.seg_dict['region_max_areas'] = [-999]
+        self.seg_dict['vertices'] = [(9,8),(7,8)]
+        self.seg_dict['triangles'] = []
+        self.seg_dict['segments'] = [[0,1]]
+        self.seg_dict['triangle_tags'] = []
+        self.seg_dict['vertex_attributes'] = [[],[]]
+        self.seg_dict['triangle_neighbors'] = []
+        self.seg_dict['segment_tags'] = ['']
+        self.seg_dict['vertex_attribute_titles'] = []
+        
+        self.reg_dict ={}
+        self.reg_dict['outline_segments'] = [[0,1]]
+        self.reg_dict['outline_segment_tags'] = ['']
+        self.reg_dict['holes'] = []
+        self.reg_dict['points'] = [(9,8),(7,8)]
+        self.reg_dict['point_attributes'] = [[],[]] 
+        self.reg_dict['regions'] = [(5,4)]
+        self.reg_dict['region_tags'] = ['']
+        self.reg_dict['region_max_areas'] = []
+        self.reg_dict['vertices'] = [(9,8),(7,8)]
+        self.reg_dict['triangles'] = []
+        self.reg_dict['segments'] = [[0,1]]
+        self.reg_dict['triangle_tags'] = []
+        self.reg_dict['vertex_attributes'] = [[],[]]
+        self.reg_dict['triangle_neighbors'] = []
+        self.reg_dict['segment_tags'] = ['']
+        self.reg_dict['vertex_attribute_titles'] = []
+        
+        self.triangle_tags_dict ={}
+        self.triangle_tags_dict['outline_segments'] = [(0, 1), (1, 2), (0, 2), (0, 3)]
+        self.triangle_tags_dict['outline_segment_tags'] = ['50', '40', '30', '20']
+        self.triangle_tags_dict['holes'] = [(0.2, 0.6)]
+        self.triangle_tags_dict['point_attributes'] = [[5, 2], [4, 2], [3, 2], [2,2]]
+        self.triangle_tags_dict['regions'] = [(0.3, 0.3),(0.3, 0.4)]
+        self.triangle_tags_dict['region_tags'] = ['1.3', 'yeah']
+        self.triangle_tags_dict['region_max_areas'] = [36.0,-7.1]
+        self.triangle_tags_dict['points'] = [(0.0, 0.0), (0.0, 4.0), (4.0, 0.0), (1.0, 1.0)]
+        self.triangle_tags_dict['vertices'] = [(0.0, 0.0), (0.0, 4.0),
+                                 (4.0, 0.0), (1.0, 1.0), (2.0, 2.0)]
+        self.triangle_tags_dict['triangles'] = [(3, 2, 4), (1, 0, 3),
+                                             (3, 4,1), (2, 3, 0)]
+        self.triangle_tags_dict['segments'] = [(0, 1), (1, 4), (2, 0),
+                                            (0, 3), (4, 2)]
+        self.triangle_tags_dict['triangle_tags'] = ['yeah', '1.3',
+                                      '1.3', '']
+        self.triangle_tags_dict['vertex_attributes'] = [[1.2,2.], [1.2,2.],
+                                          [1.2,2.], [1.2,2.], [1.2,3.]]
+        self.triangle_tags_dict['triangle_neighbors'] = [[-1, 2, 3], [3, 2, -1],
+                                           [-1, 1, 0], [1, -1, 0]]
+        self.triangle_tags_dict['segment_tags'] = ['50', '40', '30', '20', '40']
+        self.triangle_tags_dict['vertex_attribute_titles'] = ['bed elevation', 'height']
+        self.triangle_tags_dict['geo_reference'] = Geo_reference(56,1.9,1.9)
+        
+    def tearDown(self):
+        pass
+
+  ############### .TSH ##########      
+    def test_export_mesh_file(self):
+        import os
+        import tempfile
+        
+        meshDict = self.dict
+        fileName = tempfile.mktemp(".tsh")
+        export_mesh_file(fileName, meshDict)
+        loadedDict = import_mesh_file(fileName)
+        
+        #print "*(*( meshDict"
+        #print meshDict
+        #print "*(*(  loadedDcit"
+        #print loadedDict
+        #print "*(*("
+        
+        self.failUnless(array(meshDict['vertices'])  ==
+                        array(loadedDict['vertices']),
+                         'test_export_mesh_file failed. Test 1')
+        self.failUnless(array(meshDict['triangles'])  ==
+                        array(loadedDict['triangles']),
+                         'test_export_mesh_file failed. Test 2')
+        self.failUnless(array(meshDict['segments'])  ==
+                        array(loadedDict['segments']),
+                         'test_export_mesh_file failed. Test 3')
+        self.failUnless(array(meshDict['triangle_tags'])  ==
+                        array(loadedDict['triangle_tags']),
+                         'test_export_mesh_file failed. Test 4')
+        
+        self.failUnless(meshDict['vertex_attributes']  ==
+                        loadedDict['vertex_attributes'],
+                         'test_export_mesh_file failed. Test 5')
+        self.failUnless(array(meshDict['triangle_neighbors'])  ==
+                        array(loadedDict['triangle_neighbors']),
+                         'test_export_mesh_file failed. Test 6')
+        self.failUnless(array(meshDict['segment_tags'])  ==
+                        array(loadedDict['segment_tags']),
+                         'test_export_mesh_file failed. Test 7')
+        self.failUnless(array(meshDict['vertex_attribute_titles'])  ==
+                        array(loadedDict['vertex_attribute_titles']),
+                         'test_export_mesh_file failed. Test 8')
+        self.failUnless(array(meshDict['geo_reference'])  ==
+                        array(loadedDict['geo_reference']),
+                         'test_export_mesh_file failed. Test 9')
+            
+        os.remove(fileName)
+ 
+    def test_read_write_tsh_file(self):
+        dict = self.dict.copy()
+        fileName = tempfile.mktemp(".tsh")
+        export_mesh_file(fileName,dict)
+        loaded_dict = import_mesh_file(fileName)
+        os.remove(fileName)
+        dict = self.dict
+        #print "*********************"
+        #print dict
+        #print "**loaded_dict*******************"
+        #print loaded_dict
+        #print "*********************"       
+        self.check_mesh_dicts(loaded_dict,dict, 'test_read_write_msh_file')
+        
+    def test_read_write_tsh_fileII(self):
+        dict = self.sparse_dict.copy()
+        fileName = tempfile.mktemp(".tsh")
+        export_mesh_file(fileName,dict)
+        loaded_dict = import_mesh_file(fileName)
+        dict = self.sparse_dict   
+        self.check_mesh_dicts(loaded_dict,dict, 'test_read_write_msh_file')
+        os.remove(fileName)
+        
+    def test_read_write_tsh_fileIII(self):
+        dict = self.blank_dict.copy()
+        fileName = tempfile.mktemp(".tsh")
+        export_mesh_file(fileName,dict)
+        loaded_dict = import_mesh_file(fileName)
+        os.remove(fileName)
+        dict = self.blank_dict
+        #print "*********************"
+        #print dict
+        #print "**loaded_dict*******************"
+        #print loaded_dict
+        #print "*********************"       
+        self.check_mesh_dicts(loaded_dict,dict, 'test_read_write_msh_fileIII')
+
+    def test_read_write_tsh_file4(self):
+        dict = self.seg_dict.copy()
+        fileName = tempfile.mktemp(".tsh")
+        export_mesh_file(fileName,dict)
+        loaded_dict = import_mesh_file(fileName)
+        os.remove(fileName)
+        dict = self.seg_dict   
+        self.check_mesh_dicts(loaded_dict,dict, 'test_read_write_msh_file4')
+
+    def test_read_write_tsh_file5(self):
+        dict = self.triangle_tags_dict.copy()
+        fileName = tempfile.mktemp(".tsh")
+        export_mesh_file(fileName,dict)
+        loaded_dict = import_mesh_file(fileName)
+        dict = self.triangle_tags_dict
+        #print "*********************"
+        #print dict
+        #print "**loaded_dict*******************"
+        #print loaded_dict
+        #print "*********************"       
+        self.check_mesh_dicts(loaded_dict,dict, 'test_read_write_msh_file5')  
+        os.remove(fileName)
+
+    def test_read_write_tsh_file6(self):
+        dict = self.tri_dict.copy()
+        fileName = tempfile.mktemp(".tsh")
+        export_mesh_file(fileName,dict)
+        loaded_dict = import_mesh_file(fileName)
+        dict = self.tri_dict    
+        self.check_mesh_dicts(loaded_dict,dict, 'test_read_write_msh_file6')  
+        os.remove(fileName)
+        
+########################## BAD .TSH ##########################
+
+    def test_load_bad_no_file_tsh(self):
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".tsh")
+        #print fileName
+        try:
+            dict = import_mesh_file(fileName)
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'imaginary file did not raise error!')
+         
+    def test_read_write_tsh_file_bad(self):
+        dict = self.tri_dict.copy()
+        fileName = tempfile.mktemp(".xxx")
+        try:
+            export_mesh_file(fileName,dict)
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad tsh file did not raise error!')       
+        
+    def test_import_tsh_bad(self):
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".tsh")
+        file = open(fileName,"w")
+        #   this is  a bad tsh file
+        file.write("elevn\n\
+1.0 what \n\
+0.0 the \n\
+1.0 !!! \n")
+        file.close()
+        #print fileName
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad tsh file did not raise error!')       
+        os.remove(fileName)
+
+    def test_import_tsh3(self):
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".tsh")
+        file = open(fileName,"w")
+        file.write("1.0 \n\
+showme1.0 0.0 10.0 \n\
+0.0 1.0\n\
+13.0 \n")
+        file.close()
+        #print fileName
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad tsh file did not raise error!')
+        
+        os.remove(fileName)         
+
+            
+  ############### .MSH ##########
+        
+    def test_read_write_msh_file(self):
+        dict = self.dict.copy()
+        fileName = tempfile.mktemp(".msh")
+        export_mesh_file(fileName,dict)
+        loaded_dict = loadASCII._read_msh_file(fileName)
+        os.remove(fileName)
+        dict = self.dict
+        #print "*********************"
+        #print dict
+        #print "**loaded_dict*******************"
+        #print loaded_dict
+        #print "*********************"
+        self.check_mesh_dicts(loaded_dict,dict,'test_read_write_msh_file') 
+
+    def test_read_write_msh_fileII(self):
+        dict = self.sparse_dict.copy()
+        fileName = tempfile.mktemp(".msh")
+        export_mesh_file(fileName,dict)
+        loaded_dict = loadASCII._read_msh_file(fileName)
+        os.remove(fileName)
+        dict = self.sparse_dict
+        #print "*********************"
+        #print dict
+        #print "**loaded_dict*******************"
+        #print loaded_dict
+        #print "*********************"       
+        self.check_mesh_dicts(loaded_dict,dict, 'test_read_write_msh_fileII')
+              
+    def test_read_write_msh_fileIII(self):
+        dict = self.blank_dict.copy()
+        fileName = tempfile.mktemp(".msh")
+        export_mesh_file(fileName,dict)
+        loaded_dict = loadASCII._read_msh_file(fileName)
+        os.remove(fileName)
+        dict = self.blank_dict
+        #print "*********************"
+        #print dict
+        #print "**loaded_dict*******************"
+        #print loaded_dict
+        #print "*********************"       
+        self.check_mesh_dicts(loaded_dict,dict, 'test_read_write_msh_fileIII')
+       
+    def test_read_write_msh_file4(self):
+        dict = self.seg_dict.copy()
+        fileName = tempfile.mktemp(".msh")
+        export_mesh_file(fileName,dict)
+        loaded_dict = loadASCII._read_msh_file(fileName)
+        os.remove(fileName)
+        dict = self.seg_dict
+        #print "*********************"
+        #print dict
+        #print "**loaded_dict*******************"
+        #print loaded_dict
+        #print "*********************" 
+        self.check_mesh_dicts(loaded_dict,dict, 'test_read_write_msh_fileIII')
+        
+    def test_read_write_msh_file5(self):
+        dict = self.triangle_tags_dict.copy()
+        fileName = tempfile.mktemp(".msh")
+        export_mesh_file(fileName,dict)
+        loaded_dict = loadASCII._read_msh_file(fileName)
+        os.remove(fileName)
+        dict = self.triangle_tags_dict
+        #print "msh_file5*********************"
+        #print dict
+        #print "**loaded_dict*******************"
+        #print loaded_dict
+        #print "*********************" 
+        self.check_mesh_dicts(loaded_dict,dict, 'test_read_write_msh_fileIII')
+                                       
+        
+    def test_read_write_msh_file6(self):
+        dict = self.tri_dict.copy()
+        fileName = tempfile.mktemp(".msh")
+        export_mesh_file(fileName,dict)
+        loaded_dict = loadASCII._read_msh_file(fileName)
+        os.remove(fileName)
+        dict = self.tri_dict
+        #print "*********************"
+        #print dict
+        #print "**loaded_dict*******************"
+        #print loaded_dict
+        #print "*********************" 
+        self.check_mesh_dicts(loaded_dict,dict, 'test_read_write_msh_fileIII')
+                         
+    def check_mesh_dicts(self, loaded_dict, dict, fail_string ):
+        assert allclose(array(loaded_dict['points']),
+                        array(dict['points']))
+                        
+	assert allclose(array(loaded_dict['point_attributes']),
+                        array(dict['point_attributes']))
+        assert allclose(array(loaded_dict['outline_segments']),
+                        array(dict['outline_segments']))
+        
+        self.failUnless(loaded_dict['outline_segment_tags']  ==
+                        dict['outline_segment_tags'],
+                        fail_string + ' failed!! Test 4')
+        
+            
+        assert allclose(array(loaded_dict['regions']),
+                        array(dict['regions']))
+        self.failUnless(loaded_dict['region_tags']  ==
+                        dict['region_tags'],
+                        fail_string + ' failed!! Test 5')
+        
+        assert allclose(array(loaded_dict['region_max_areas']),
+                        array(dict['region_max_areas']))
+ 
+        assert allclose(array(loaded_dict['holes']),
+                        array(dict['holes']))
+ 
+        assert allclose(array(dict['vertices']),
+                        array(loaded_dict['vertices']))
+ 
+        assert allclose(array(dict['triangles']),
+                        array(loaded_dict['triangles']))
+ 
+        assert allclose(array(dict['segments']),
+                        array(loaded_dict['segments']))
+        for ob, ldob in map(None,dict['triangle_tags'],
+                              loaded_dict['triangle_tags']):
+            self.failUnless(ob  == ldob,
+                        fail_string + ' failed!! Test triangle_tags')
+        self.failUnless(loaded_dict['vertex_attributes']  ==
+                        dict['vertex_attributes'],
+                        fail_string + ' failed!! Test vertex_attributes')
+ 
+        assert allclose(array(dict['triangle_neighbors']),
+                        array(loaded_dict['triangle_neighbors']))
+
+        for seg, ldseg in map(None,dict['segment_tags'],
+                              loaded_dict['segment_tags']):
+            self.failUnless(seg  == ldseg,
+                        fail_string + ' failed!! Test 8')
+        try:
+            assert allclose(array(dict['vertex_attribute_titles']),
+                            array(loaded_dict['vertex_attribute_titles']))
+        except TypeError:
+            self.failUnless(array(loaded_dict['vertex_attribute_titles'])  ==
+                        array(dict['vertex_attribute_titles']),
+                        fail_string + ' failed!! Test 8')
+        try:    
+            self.failUnless(loaded_dict['geo_reference']  ==
+                        dict['geo_reference'] ,
+                        fail_string + ' failed!! Test geo_reference')
+        except KeyError:          
+            self.failUnless(not dict.has_key('geo_reference' and
+                         loaded_dict['geo_reference'] == None)   ,
+                        fail_string + ' failed!! Test geo_reference')
+  
+########################## BAD .MSH ##########################          
+
+    def test_load_bad_no_file_msh(self):
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".msh")
+        #print fileName
+        try:
+            dict = import_mesh_file(fileName)
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'imaginary file did not raise error!')
+            
+    def test_import_msh_bad(self):
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".msh")
+        file = open(fileName,"w")
+        #   this is  a bad tsh file
+        file.write("elevn\n\
+1.0 what \n\
+0.0 the \n\
+1.0 !!! \n")
+        file.close()
+        #print fileName
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad msh file did not raise error!')        
+        os.remove(fileName)         
+   
+  ###################### .XYA ##############################
+        
+    def test_export_xya_file(self):
+        dict = {}
+        att_dict = {}
+	dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	att_dict['elevation'] = array([10.0, 0.0, 10.4])
+	att_dict['brightness'] = array([10.0, 0.0, 10.4])
+        dict['attributelist'] = att_dict
+        dict['geo_reference'] = Geo_reference(56,1.9,1.9)
+        
+        
+        fileName = tempfile.mktemp(".xya")
+        export_points_file(fileName, dict)
+        dict2 = import_points_file(fileName)
+        #print "fileName",fileName 
+        os.remove(fileName)
+        #print "dict2",dict2 
+        
+	assert allclose(dict2['pointlist'],[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	assert allclose(dict2['attributelist']['elevation'], [10.0, 0.0, 10.4])
+        answer = [10.0, 0.0, 10.4]
+	assert allclose(dict2['attributelist']['brightness'], answer)
+        #print "dict2['geo_reference']",dict2['geo_reference'] 
+        self.failUnless(dict['geo_reference'] == dict2['geo_reference'],
+                         'test_writepts failed. Test geo_reference')
+
+    def test_export_xya_file2(self):
+        dict = {}
+        att_dict = {}
+	dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	att_dict['elevation'] = array([10.0, 0.0, 10.4])
+	att_dict['brightness'] = array([10.0, 0.0, 10.4])
+        dict['attributelist'] = att_dict
+        
+        
+        fileName = tempfile.mktemp(".xya")
+        export_points_file(fileName, dict)
+        dict2 = import_points_file(fileName)
+        #print "fileName",fileName 
+        os.remove(fileName)
+        #print "dict2",dict2 
+        
+	assert allclose(dict2['pointlist'],[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	assert allclose(dict2['attributelist']['elevation'], [10.0, 0.0, 10.4])
+        answer = [10.0, 0.0, 10.4]
+	assert allclose(dict2['attributelist']['brightness'], answer)
+
+        
+    def test_loadxya(self):
+        """
+        comma delimited
+        """
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("elevation  , speed \n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.4, 40.0\n")
+        file.close()
+        #print fileName
+        dict = import_points_file(fileName,delimiter = ',')
+        os.remove(fileName)
+	assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
+	assert allclose(dict['attributelist']['speed'], [0.0, 10.0, 40.0])
+
+    #FIXME - redundant test?
+    def test_loadxy(self):
+        """
+        To test the mesh side of loading xya files.
+        Not the loading of xya files
+        """
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("elevation speed \n\
+1.0 0.0 10.0 0.0\n\
+0.0 1.0 0.0 10.0\n\
+1.0 0.0 10.4 40.0\n")
+        file.close()
+        #print fileName
+        dict = import_points_file(fileName)
+        os.remove(fileName)
+	assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
+	assert allclose(dict['attributelist']['speed'], [0.0, 10.0, 40.0])
+
+               
+    def test_loadxya2(self):
+        """
+        space delimited
+        """
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  elevation   speed \n\
+1.0 0.0 10.0 0.0\n\
+0.0 1.0 0.0 10.0\n\
+1.0 0.0 10.4 40.0\n")
+        file.close()
+        #print fileName
+        dict = import_points_file(fileName,delimiter = ' ')
+        os.remove(fileName)
+	assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
+	assert allclose(dict['attributelist']['speed'], [0.0, 10.0, 40.0])
+     
+    def test_loadxya3(self):
+        """
+        space delimited
+        """
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  elevation   speed \n\
+1.0 0.0 10.0 0.0\n\
+0.0 1.0 0.0 10.0\n\
+1.0 0.0 10.4 40.0\n\
+#geocrap\n\
+56\n\
+56.6\n\
+3\n")
+        file.close()
+        #print fileName
+        dict = import_points_file(fileName,delimiter = ' ')
+        os.remove(fileName)
+	assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
+	assert allclose(dict['attributelist']['speed'], [0.0, 10.0, 40.0])
+     
+        geo_reference = Geo_reference(56, 56.6, 3.0)
+        
+        self.failUnless(geo_reference == dict['geo_reference'],
+                         'test_writepts failed. Test geo_reference')
+        
+########################## BAD .XYA ##########################
+  
+    def test_loadxy_bad_no_file_xya(self):
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".xya")
+        #print fileName
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'imaginary file did not raise error!')
+  
+    def test_read_write_points_file_bad(self):
+        dict = self.tri_dict.copy()
+        fileName = tempfile.mktemp(".xxx")
+        try:
+            export_points_file(fileName,dict)
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad points file extension did not raise error!')
+     
+    def test_read_write_points_file_bad2(self):
+        dict = {}
+        att_dict = {}
+	dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	att_dict['elevation'] = array([10.0, 0.0, 10.4])
+	att_dict['brightness'] = array([10.0, 0.0, 10.4])
+        dict['attributelist'] = att_dict
+        dict['geo_reference'] = Geo_reference(56,1.9,1.9)
+        try:
+            export_points_file("_???/yeah.xya",dict)
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad points file extension did not raise error!')
+                   
+    def test_loadxy_bad(self):
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  elevation   \n\
+1.0 0.0 10.0 0.0\n\
+0.0 1.0 0.0 10.0\n\
+1.0 0.0 10.4 40.0\n")
+        file.close()
+        #print fileName
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad xya file did not raise error!')
+        os.remove(fileName)
+       
+    def test_loadxy_bad2(self):
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("elevation\n\
+1.0 0.0 10.0 \n\
+0.0 1.0\n\
+1.0 \n")
+        file.close()
+        #print fileName
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad xya file did not raise error!')    
+        os.remove(fileName)
+   
+    def test_loadxy_bad3(self):
+        """
+        specifying wrong delimiter
+        """
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  elevation  , speed \n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.4, 40.0\n")
+        file.close()
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad xya file did not raise error!')    
+        os.remove(fileName)
+     
+    def test_loadxy_bad4(self):
+        """
+        specifying wrong delimiter
+        """
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  elevation   speed \n\
+1.0 0.0 10.0 0.0\n\
+0.0 1.0 0.0 10.0\n\
+1.0 0.0 10.4 40.0\n\
+yeah")
+        file.close()
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad xya file did not raise error!')    
+        os.remove(fileName)
+  
+    def test_loadxy_bad4(self):
+        """
+        specifying wrong delimiter
+        """
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  elevation   speed \n\
+1.0 0.0 10.0 0.0\n\
+0.0 1.0 0.0 10.0\n\
+1.0 0.0 10.4 40.0\n\
+#geocrap")
+        file.close()
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad xya file did not raise error!')    
+        os.remove(fileName)
+
+    def test_loadxy_bad5(self):
+        """
+        specifying wrong delimiter
+        """
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("  elevation   speed \n\
+1.0 0.0 10.0 0.0\n\
+0.0 1.0 0.0 10.0\n\
+1.0 0.0 10.4 40.0\n\
+#geocrap\n\
+crap")
+        file.close()
+        try:
+            dict = import_points_file(fileName,delimiter = ' ')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'bad xya file did not raise error!')    
+        os.remove(fileName)        
+  ############### .PTS ##########
+            
+    def test_loadpts(self):
+        
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        fileName = tempfile.mktemp(".pts")
+        # NetCDF file definition
+        outfile = NetCDFFile(fileName, 'w')
+        
+        # dimension definitions
+        outfile.createDimension('number_of_points', 3)    
+        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
+    
+        # variable definitions
+        outfile.createVariable('points', Float, ('number_of_points',
+                                                 'number_of_dimensions'))
+        outfile.createVariable('elevation', Float, ('number_of_points',))
+    
+        # Get handles to the variables
+        points = outfile.variables['points']
+        elevation = outfile.variables['elevation']
+ 
+        points[0, :] = [1.0,0.0]
+        elevation[0] = 10.0 
+        points[1, :] = [0.0,1.0]
+        elevation[1] = 0.0  
+        points[2, :] = [1.0,0.0]
+        elevation[2] = 10.4    
+
+        outfile.close()
+        
+        dict = import_points_file(fileName)
+        os.remove(fileName)
+        answer =  [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
+        #print "dict['pointlist']",dict['pointlist']
+        #print "answer",answer
+	assert allclose(dict['pointlist'], [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	assert allclose(dict['attributelist']['elevation'], [10.0, 0.0, 10.4])
+ 
+    def test_writepts(self):
+        dict = {}
+        att_dict = {}
+	dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	att_dict['elevation'] = array([10.0, 0.0, 10.4])
+	att_dict['brightness'] = array([10.0, 0.0, 10.4])
+        dict['attributelist'] = att_dict
+        dict['geo_reference'] = Geo_reference(56,1.9,1.9)
+        
+        
+        fileName = tempfile.mktemp(".pts")
+        export_points_file(fileName, dict)
+        dict2 = import_points_file(fileName)
+        #print "fileName",fileName 
+        os.remove(fileName)
+        #print "dict2",dict2 
+        
+	assert allclose(dict2['pointlist'],[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	assert allclose(dict2['attributelist']['elevation'], [10.0, 0.0, 10.4])
+        answer = [10.0, 0.0, 10.4]
+	assert allclose(dict2['attributelist']['brightness'], answer)
+
+        #print "dict['geo_reference'] ",dict['geo_reference'] 
+        #print "dict2['geo_reference']",dict2['geo_reference']
+        
+        self.failUnless(dict['geo_reference'] == dict2['geo_reference'],
+                         'test_writepts failed. Test geo_reference')
+        
+ ########################## BAD .PTS ##########################          
+
+    def test_load_bad_no_file_pts(self):
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".pts")
+        #print fileName
+        try:
+            dict = import_points_file(fileName)
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,
+                        'imaginary file did not raise error!')
+            
+  ############### .PTS OTHER ##########
+
+    def test_concatinate_attributelist(self):
+        dic = {}
+        dic['one'] = array([1,2])
+        dic['2'] = array([2,7])
+        dic['three'] = array([3,79])
+        dic['4'] = array([4,47])
+        dic['five'] = array([5,17])
+        titles, block = concatinate_attributelist(dic)
+        #print "titles", titles
+        #print "array", block
+        self.failUnless(titles  == ['4', '2', 'five', 'three', 'one'],
+                         'test_concatinate_attributelist failed.')
+	assert allclose(block, [[4,2,5,3,1],[47,7,17,79,2]])
+
+    def test_half_pts(self):
+        dict = {}
+        att_dict = {}
+	dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	att_dict['elevation'] = array([10.0, 0.0, 10.4])
+	att_dict['brightness'] = array([10.0, 0.0, 10.4])
+        dict['attributelist'] = att_dict
+        out_dict = half_pts(dict)
+	assert allclose(out_dict['pointlist'],[[1.0, 0.0],[1.0, 0.0]])
+	assert allclose(out_dict['attributelist']['elevation'], [10.0, 10.4])
+	assert allclose(out_dict['attributelist']['brightness'], [10.0, 10.4])
+  
+    def test_extent_point_atts(self):
+        
+        dict = {}
+        att_dict = {}
+	dict['pointlist'] = array([[1.0, 10.0],[0.0, 1.0],[10.0, -10.0]])
+	att_dict['elevation'] = array([30.0, 0.0, 10.4])
+	att_dict['brightness'] = array([10.0, 0.0, 10.4])
+        dict['attributelist'] = att_dict
+        out_dict = extent_point_atts(dict)
+        
+        #print "out_dict['pointlist']",out_dict #['pointlist'] 
+	assert allclose(out_dict['pointlist'],[[0.0, -10.0],[10.0, -10.0],
+                                               [10.0,10.0],[0.0, 10.0]])
+
+        self.failUnless(dict['attributelist']  == {},
+                         'test_extent_point_atts failed. Test 1')
+        
+    def test_reduce_pts(self):
+        dict = {}
+        att_dict = {}
+	dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	att_dict['elevation'] = array([10.0, 0.0, 10.4])
+	att_dict['brightness'] = array([10.0, 0.0, 10.4])
+        dict['attributelist'] = att_dict
+        
+        inFileName = tempfile.mktemp(".pts")
+        export_points_file(inFileName, dict)
+        
+        outFileName = tempfile.mktemp(".pts")
+        
+        dict2 = reduce_pts(inFileName,outFileName, 1 )
+        os.remove(inFileName)
+
+        dict2 = import_points_file(outFileName)
+        os.remove(outFileName)
+        #print "dict2",dict2 
+        
+	assert allclose(dict2['pointlist'],[[1.0, 0.0]])
+	assert allclose(dict2['attributelist']['elevation'], [10.0])
+	assert allclose(dict2['attributelist']['brightness'], [10.0])
+       
+    def test_produce_half_point_files(self):
+        dict = {}
+        att_dict = {}
+	dict['pointlist'] = array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+	att_dict['elevation'] = array([10.0, 0.0, 10.4])
+	att_dict['brightness'] = array([10.0, 0.0, 10.4])
+        dict['attributelist'] = att_dict
+        
+        inFileName = tempfile.mktemp(".pts")
+        export_points_file(inFileName, dict)
+        
+        outFileName = tempfile.mktemp(".pts")
+        delimiter = '_p'
+        outfiles = produce_half_point_files(inFileName, 1,
+                                            delimiter, verbose = False )
+        os.remove(inFileName)
+        root, ext = splitext(inFileName)
+        outFileName = root + delimiter + ext
+        #print "outFileName",outfiles  
+        dict2 = import_points_file(outfiles[1])
+        for file in outfiles:
+            #print "del file",file 
+            os.remove(file)
+        
+	assert allclose(dict2['pointlist'],[[1.0, 0.0]])
+	assert allclose(dict2['attributelist']['elevation'], [10.0])
+	assert allclose(dict2['attributelist']['brightness'], [10.0])
+                      
+#-------------------------------------------------------------
+if __name__ == "__main__":
+
+    suite = unittest.makeSuite(loadASCIITestCase,'test')
+    #suite = unittest.makeSuite(loadASCIITestCase,'test_writepts')
+    runner = unittest.TextTestRunner(verbosity=0)
+    runner.run(suite)
+    
Index: undation/pmesh/test_triangmodule.py
===================================================================
--- /inundation/pmesh/test_triangmodule.py	(revision 2252)
+++ 	(revision )
@@ -1,305 +1,0 @@
-#!/usr/bin/env python
-#
-
-import sys
-
-import unittest
-import triang
-
-
-class triangTestCase(unittest.TestCase):
-    def setUp(self):
-        pass
-        
-    def tearDown(self):
-        pass
- 
-    def testrectangle(self):
-
-        points = []
-        seglist = []
-        segattlist = []
-        holelist = []
-        regionlist = []
-
-	points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0)]
-	# This will cause a seg fault
-        #points = [[0.0,0.0],[0.0,10.0],[3.0,0.0],[3.0,10.0]]
-
-        # automatically generate pointattlist
-	pointattlist = []
-        for point in points:
-            pointattlist.append([])
-        #print "pointattlist",pointattlist
-	#pointattlist = [[],[],[],[]]
-        trilist = []
-        mode = "Qzcn"
-        data = triang.genMesh(points,seglist,holelist,regionlist,
-                              pointattlist,segattlist,trilist, mode)
-        self.failUnless(data['generatedtrianglelist'] ==[(1, 0, 2), (2, 3, 1)],
-                        'generatedtrianglelist is wrong!')
-        #print "data['generatedsegmentlist']",data['generatedsegmentlist']
-        self.failUnless(data['generatedsegmentlist'] ==[(2, 0), (3, 2), (1, 3),
-                                                        (0, 1)]
-              ,          'generatedsegmentlist is wrong!')
-        self.failUnless(data['generatedpointlist'] ==[(0.0, 0.0), (0.0, 10.0),
-                                                      (3.0, 0.0), (3.0, 10.0)],
-                        ' is wrong!')
-        
-        #print "data['generatedtriangleneighborlist']",
-        #data['generatedtriangleneighborlist'] 
-        self.failUnless(data['generatedtriangleneighborlist'] ==[(-1, 1, -1),
-                                                                 (-1, 0, -1)],
-                        'generatedtriangleneighborlist is wrong!')
-
-        
-    def testrectangleII(self):
-
-        points = []
-        seglist = []
-        holelist = []
-        regionlist = []
-
-	points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0)]
-	pointattlist = [[],[],[],[]]
-	regionlist.append( (1.2,1.2,5.0) )
-	seglist = [(0,1),(1,3),(3,2),(2,0)]
-        segattlist = [0,0,0,0]
-        trilist = []
-        mode = "Qzp"
-        data = triang.genMesh(points,seglist,holelist,regionlist,
-                              pointattlist,segattlist,trilist, mode)
-        
-        self.failUnless(data['generatedtrianglelist'] ==[(1, 0, 2), (2, 3, 1)],
-                        'trianglelist is wrong!')
-        self.failUnless(data['generatedsegmentlist'] ==[(0, 1), (1, 3),
-                                                        (3, 2), (2, 0)],
-                        'segmentlist is wrong!')
-        self.failUnless(data['generatedpointlist'] ==[(0.0, 0.0), (0.0, 10.0),
-                                                      (3.0, 0.0), (3.0, 10.0)],
-                        ' is wrong!')
-   
-    def testsegmarker(self):
-
-        points = []
-        seglist = []
-        holelist = []
-        regionlist = []
-
-	points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0)]
-	pointattlist = [[],[],[],[]]
-	regionlist.append( (1.2,1.2,5.0) )
-	seglist = [(0,1),(1,3),(3,2),(2,0)]
-        segattlist = [1.0,2.0,3.0,4.0]
-        trilist = []
-        mode = "Qzp"
-        data = triang.genMesh(points,seglist,holelist,regionlist,
-                              pointattlist,segattlist,trilist, mode)
-        
-        self.failUnless(data['generatedtrianglelist'] ==[(1, 0, 2), (2, 3, 1)],
-                        'trianglelist is wrong!')
-        self.failUnless(data['generatedsegmentlist'] ==[(0, 1), (1, 3),
-                                                        (3, 2), (2, 0)],
-                        'segmentlist is wrong!')
-        self.failUnless(data['generatedpointlist'] ==[(0.0, 0.0), (0.0, 10.0),
-                                                      (3.0, 0.0), (3.0, 10.0)],
-                        ' is wrong!')
-        self.failUnless(data['generatedsegmentmarkerlist'] ==[1,2,3,4],
-                        ' is wrong!')
-             
-    def testbad_region(self):
-
-        points = []
-        seglist = []
-        holelist = []
-        regionlist = [(1.2,1.2)]
-
-	points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0)]
-	pointattlist = [[],[],[],[]]
-	seglist = [(0,1),(1,3),(3,2),(2,0)]
-        segattlist = [0,0,0,0]
-        trilist = []
-        mode = "Qzpn"
-        try:
-            data = triang.genMesh(points,seglist,holelist,
-                                  regionlist,pointattlist,segattlist,
-                                  trilist, mode)
-            
-        except TypeError:
-            pass
-        else:
-            self.failUnless(0 ==1,
-                        'bad region list did not raise error!')
-      
-    def testregion_with_maxarea(self):
-
-        points = []
-        seglist = []
-        holelist = []
-        regionlist = [(3,1,1.0)]
-
-	points = [(0.0,0.0),(6.0,0.0),(6.0,6.0),(0.0,6.0)]
-	pointattlist = [[],[],[],[]]
-	seglist = [(0,1),(1,2),(3,2),(3,0),(0,2)]
-        segattlist = [0,0,0,0,0]
-        trilist = []
-        mode = "Qzpna36a"
-        data = triang.genMesh(points,seglist,holelist,regionlist,
-                              pointattlist,segattlist, trilist, mode)
-        
-        self.failUnless(len(data['generatedtrianglelist']) == 2,
-                        'testregion_with_maxarea 1: # of tris is wrong!')
-        ## Another test case    
-        regionlist = [(3,1,1.0),(1,3,1.0,8.0)]
-        mode = "Qzp21na36a"
-        data = triang.genMesh(points,seglist,holelist,regionlist,
-                              pointattlist,segattlist, trilist, mode)
-        #print "len(data['generatedtrianglelist']",len(data['generatedtrianglelist'])
-        # for Duncan On unix this returns a 7 triangle result.
-        # for Duncan on Windows returns a 6 triangle result.
-        # for Ole on nautilus this returns 6
-        # for Duncan on nautilus this returns 7
-        # ??, it seems to be the results from triangle that is
-        # causing the different results, and we are treating
-        # triangle as a back box.
-       
-        self.failUnless(len(data['generatedtrianglelist']) >= 6,
-                        'testregion_with_maxarea 2: # of tris is wrong!')
-        ## Another test case    
-        regionlist = [(3,1,1.0,8.0),(1,3,1.0,8.0)]
-        mode = "Qzpna36a"
-        data = triang.genMesh(points,seglist,holelist,regionlist,
-                              pointattlist,segattlist, trilist, mode)
-        #print "len(data['generatedtrianglelist']",len(data['generatedtrianglelist'])
-        # On unix this returns a 10 triangle result.
-        # Windows returns a 8 triangle result.
-        self.failUnless(len(data['generatedtrianglelist']) >= 8,
-                        'testregion_with_maxarea 3: # of tris is wrong!')
-                  
-        ## Another test case    
-        regionlist = [(3,1,1.0),(1,3,1.0,8.0)]
-        mode = "Qzpna8a"
-        data = triang.genMesh(points,seglist,holelist,regionlist,
-                              pointattlist,segattlist, trilist, mode)
-        #print "len(data['generatedtrianglelist']",len(data['generatedtrianglelist'])
-        # On unix this returns a 10 triangle result.
-        # Windows returns a 8 triangle result.
-        
-        self.failUnless(len(data['generatedtrianglelist']) >= 8,
-                        'testregion_with_maxarea 4: # of tris is wrong!')
-                  
-    def testbad_point(self):
-
-        holelist = []
-        regionlist = []
-
-	points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0)]
-	pointattlist = [[],[],[],[]]
-	seglist = [(0,1),(1,3),(3,2),(2,0)]
-        segattlist = [0,0,0,0]
-        trilist = []
-        mode = "Qzpn"
-        try:
-            data = triang.genMesh(points,seglist,holelist,regionlist,
-                                  pointattlist,segattlist, trilist, mode)
-            
-        except TypeError:
-            pass
-        else:
-            self.failUnless(0 ==1,
-                        'bad point list did not raise error!')
-      
-    def testbad_hole(self):
-
-        holelist = [(9.0)]
-        regionlist = []
-
-	points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0)]
-	pointattlist = [[],[],[],[]]
-	seglist = [(0,1),(1,3),(3,2),(2,0)]
-        segattlist = [0,0,0,0]
-        trilist = []
-        mode = "Qzpn"
-        try:
-            data = triang.genMesh(points,seglist,holelist,regionlist,
-                                  pointattlist,segattlist, trilist, mode)
-            
-        except TypeError:
-            pass
-        else:
-            self.failUnless(0 ==1,
-                        'bad hole list did not raise error!')
-      
-    def testbad_segment(self):
-
-        holelist = []
-        regionlist = []
-
-	points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0)]
-	pointattlist = [[],[],[],[]]
-	seglist = [(0,1),(1,3),(3,2),(2)]
-        segattlist = [0,0,0,0]
-        trilist = []
-        mode = "Qzpn"
-        try:
-            data = triang.genMesh(points,seglist,holelist,regionlist,
-                                  pointattlist,segattlist, trilist, mode)
-            
-        except TypeError:
-            pass
-        else:
-            self.failUnless(0 ==1,
-                        'bad segment list did not raise error!')
-     
-    def testbad_segattlist(self):
-
-        holelist = []
-        regionlist = []
-
-	points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0)]
-	pointattlist = [[],[],[],[]]
-	seglist = [(0,1),(1,3),(3,2),(2,0)]
-        segattlist = [0,0]
-        trilist = []
-        mode = "Qzpn"
-        try:
-            data = triang.genMesh(points,seglist,holelist,regionlist,
-                                  pointattlist,segattlist, trilist, mode)
-            
-            self.failUnless(data['trianglelist'] ==[(1, 0, 2), (2, 3, 1)],
-                        'trianglelist is wrong!')
-        except TypeError:
-            pass
-        else:
-            self.failUnless(0 ==1,
-                        'bad segment attribute list did not raise error!')
-            
-    def testrectangle_regions(self):
-
-        points = []
-        seglist = []
-        holelist = []
-        regionlist = []
-
-	points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0)]
-	pointattlist = [[],[],[],[]]
-        #it seems that
-        #triangle only associates one region with a triangle
-	regionlist.append( (1.3,1.3,88.33) )
-	regionlist.append( (1.2,1.2,77,55) )
-	seglist = [(0,1),(1,3),(3,2),(2,0)]
-        segattlist = [0,0,0,0]
-        trilist = []
-        mode = "QAzpq"
-        data = triang.genMesh(points,seglist,holelist,regionlist,
-                              pointattlist,segattlist, trilist, mode)
-           
-        self.failUnless(data['generatedtriangleattributelist'] ==[[77.0], [77.0], [77.0], [77.0]],
-                        'triangleattributelist is wrong!')
-        
-        
-if __name__ == "__main__":
-
-    suite = unittest.makeSuite(triangTestCase,'test')
-    runner = unittest.TextTestRunner()
-    runner.run(suite)
Index: /inundation/pyvolution/test_all.py
===================================================================
--- /inundation/pyvolution/test_all.py	(revision 2252)
+++ /inundation/pyvolution/test_all.py	(revision 2253)
@@ -15,7 +15,5 @@
 #List files that should be excluded from the testing process.
 #E.g. if they are known to fail and under development
-#exclude = ['test_least_squares.py', 'test_cg_solve.py',
-          # 'test_interpolate_sww.py']
-#exclude = ['test_pmesh2domain.py'] 
+#exclude = ['test_geo_reference.py','test_redfearn.py', 'test_point.py'] 
 
 
@@ -71,5 +69,7 @@
 
     filenameToModuleName = lambda f: os.path.splitext(f)[0]
+    #print "files",files 
     moduleNames = map(filenameToModuleName, files)
+    #print "moduleNames",moduleNames
     modules = map(__import__, moduleNames)
     load = unittest.defaultTestLoader.loadTestsFromModule
