Index: inundation/ga/storm_surge/alpha_shape/alpha_shape.py
===================================================================
--- inundation/ga/storm_surge/alpha_shape/alpha_shape.py	(revision 579)
+++ inundation/ga/storm_surge/alpha_shape/alpha_shape.py	(revision 579)
@@ -0,0 +1,103 @@
+"""Alpha shape   
+"""
+
+import exceptions
+from Numeric import array, Float
+from load_mesh.loadASCII import load_xya_file, export_boundary_file
+
+class PointError(exceptions.Exception): pass
+
+OUTPUT_FILE_TITLE = "# The alpha shape boundary defined by point index pairs of edges"
+
+def alpha_shape_via_files(point_file, boundary_file, alpha= None):
+    
+    from load_mesh.loadASCII import load_xya_file
+    
+    point_dict = load_xya_file(point_file)
+    points = point_dict['pointlist']
+    #title_string = point_dict['title']
+    
+    alpha = Alpha_Shape(points)
+    alpha.write_boundary(boundary_file)
+
+class Alpha_Shape:
+
+    def __init__(self, points, alpha = None):
+
+        
+        """ Build interpolation matrix mapping from
+        function values at vertices to function values at data points
+
+        Inputs:
+        
+          points: List of coordinate pairs [[x1, y1],[x2, y2]..] 
+
+          alpha: alpha shape parameter
+          
+        """
+        self._set_points(points)
+        self._alpha_shape_algorithm()
+            
+
+    def _set_points(self, points):
+        """
+        """
+        
+        if len (points) <= 2:
+            raise PointError, "Too few points to find an alpha shape"
+        
+        #Convert input to Numeric arrays
+        self.points = array(points).astype(Float)
+
+        if len (points) <= 2:
+            raise PointError, "Too few points to find an alpha shape"
+        
+    
+    def write_boundary(self,file_name):
+        """
+        Write the boundary to a file
+        """
+        #print " this info will be in the file",boundary
+        export_boundary_file(file_name, self.get_boundary(),
+                             OUTPUT_FILE_TITLE, delimiter = ',')
+    
+    def get_boundary(self):
+        """
+        """
+        return self.boundary
+    
+    def _alpha_shape_algorithm(self):
+        
+        # A python style guide suggests using _[function name] to
+        # specify interal functions
+        #- this is the first time I've used it though - DSG
+        
+        # this produces a baaad boundary
+        boundary = []
+        for point_index in range(len(self.points)-1):
+            boundary.append([point_index,point_index +1])
+        boundary.append([len(self.points)-1,0])
+        self.boundary = boundary
+        
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    """
+    Load in a mesh and data points with attributes.
+    Fit the attributes to the mesh.
+    Save a new mesh file.
+    """
+    import os, sys
+    usage = "usage: %s point_file.xya boundary_file.?? [alpha]"\
+            %os.path.basename(sys.argv[0])
+
+    if len(sys.argv) < 3:
+        print usage
+    else:
+        point_file = sys.argv[1]
+        boundary_file = sys.argv[2]
+        if len(sys.argv) > 4:
+            alpha = sys.argv[3]
+        else:
+            alpha = None
+        alpha_shape_via_files(point_file, boundary_file, alpha)
+        
Index: inundation/ga/storm_surge/alpha_shape/test_alpha_shape.py
===================================================================
--- inundation/ga/storm_surge/alpha_shape/test_alpha_shape.py	(revision 579)
+++ inundation/ga/storm_surge/alpha_shape/test_alpha_shape.py	(revision 579)
@@ -0,0 +1,94 @@
+#!/usr/bin/env python
+
+#TEST
+
+import unittest
+from Numeric import allclose
+
+from alpha_shape import *
+
+class TestCase(unittest.TestCase):
+
+    def setUp(self):
+        pass
+
+         
+    def tearDown(self):
+        pass
+
+
+    def test_alpha(self):
+        a = [0.0, 0.0]
+        b = [1.0, 0.0]
+        c = [2.0, 0.0]
+        d = [2.0, 2.0]
+        e = [1.0, 2.0]
+        f = [0.0, 2.0]
+
+        alpha = Alpha_Shape([a,b,c,d,e,f])
+        result = alpha.get_boundary()
+        answer = [[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 0]]
+        #print "result",result
+        assert allclose(answer, result) 
+
+        
+    def test_not_enough_points(self):
+        a = [0.0, 0.0]
+        b = [1.0, 0.0]
+
+        try:
+            alpha = Alpha_Shape([a,b])
+        except PointError:
+            pass
+        else:
+            self.failUnless(0==1,
+                        'point list with 2 points did not raise an error!')
+
+
+    def test_alpha_stand_alone(self):
+        import os
+        import tempfile
+       
+        fileName = tempfile.mktemp(".xya")
+        file = open(fileName,"w")
+        file.write("title row \n\
+0.0, 0.0\n\
+1.0, 0.0\n\
+2.0, 0.0\n\
+2.0, 2.0\n\
+1.0, 2.0\n\
+0.0, 2.0\n")
+        file.close()
+
+        output_file_name = tempfile.mktemp(".bnd")
+        command = 'python .\\alpha_shape.py ' + fileName + ' ' + output_file_name
+        #print "command", command
+        os.system(command)
+        os.remove(fileName)
+        
+        file = open(output_file_name,"r")
+        lFile = file.read().split('\n')
+        file.close()
+        os.remove(output_file_name)
+        self.failUnless(lFile[1] == "0,1" and
+                        lFile[2] == "1,2" and 
+                        lFile[3] == "2,3" and 
+                        lFile[4] == "3,4" and 
+                        lFile[5] == "4,5" and 
+                        lFile[6] == "5,0" and 
+                        lFile[7] == ""
+                        ,
+                        'boundary file is wrong')
+
+
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    suite = unittest.makeSuite(TestCase,'test')
+    runner = unittest.TextTestRunner(verbosity=1)
+    runner.run(suite)
+
+    
+    
+
+
+
