Changeset 7276 for anuga_core/source/anuga/geospatial_data

Timestamp:

Jun 30, 2009, 2:07:41 PM (16 years ago)

Author:

ole

Message:

Merged numpy branch back into the trunk.

In ~/sandpit/anuga/anuga_core/source
svn merge -r 6246:HEAD ../../branches/numpy .

In ~/sandpit/anuga/anuga_validation
svn merge -r 6417:HEAD ../branches/numpy_anuga_validation .

In ~/sandpit/anuga/misc
svn merge -r 6809:HEAD ../branches/numpy_misc .

For all merges, I used numpy version where conflicts existed

The suites test_all.py (in source/anuga) and validate_all.py passed using Python2.5 with numpy on my Ubuntu Linux box.

Location:

anuga_core/source/anuga/geospatial_data

Files:

• geospatial_data.py (modified) (86 diffs)
• test_geospatial_data.py (modified) (70 diffs)

anuga_core/source/anuga/geospatial_data/geospatial_data.py

@@ -1 @@
-"""Class Geospatial_data - Manipulation of locations on the planet and
-associated attributes.
-
+"""Class Geospatial_data
+
+Manipulation of locations on the planet and associated attributes.
 """
 
@@ -9 +9 @@
 from warnings import warn
 from string import lower
-from RandomArray import randint, seed, get_seed
 from copy import deepcopy
+import copy
+
 from Scientific.IO.NetCDF import NetCDFFile
-
-import Numeric as num
+import numpy as num
+from numpy.random import randint, seed
 
 from anuga.coordinate_transforms.lat_long_UTM_conversion import UTMtoLL
@@ -20 +21 @@
      TitleError, DEFAULT_ZONE, ensure_geo_reference, write_NetCDF_georeference
 from anuga.coordinate_transforms.redfearn import convert_from_latlon_to_utm
+from anuga.utilities.system_tools import clean_line
 from anuga.utilities.anuga_exceptions import ANUGAError
 from anuga.config import points_file_block_line_size as MAX_READ_LINES
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
+from anuga.config import netcdf_float
+
 
 DEFAULT_ATTRIBUTE = 'elevation'
@@ -56 +60 @@
                  load_file_now=True,
                  verbose=False):
-        """
-        Create instance from data points and associated attributes
+        """Create instance from data points and associated attributes
 
         data_points: x,y coordinates in meters. Type must be either a
-        sequence of 2-tuples or an Mx2 Numeric array of floats.  A file name
+        sequence of 2-tuples or an Mx2 numeric array of floats.  A file name
         with extension .txt, .cvs or .pts can also be passed in here.
 
@@ -131 +134 @@
 
         verbose:
-
         """
 
         if isinstance(data_points, basestring):
-            # assume data point is really a file name
+            # assume data_points is really a file name
             file_name = data_points
 
         self.set_verbose(verbose)
-        self.geo_reference = None #create the attribute
+        self.geo_reference = None
         self.file_name = file_name
 
@@ -148 +150 @@
 
         if file_name is None:
-            if latitudes is not None \
-               or longitudes is not None \
-               or points_are_lats_longs:
+            if (latitudes is not None or longitudes is not None
+                    or points_are_lats_longs):
                 data_points, geo_reference = \
                     _set_using_lat_long(latitudes=latitudes,
@@ -156 +157 @@
                                         geo_reference=geo_reference,
                                         data_points=data_points,
-                                        points_are_lats_longs=\
+                                        points_are_lats_longs=
                                             points_are_lats_longs)
             self.check_data_points(data_points)
@@ -188 +189 @@
                         # This message was misleading.
                         # FIXME (Ole): Are we blocking here or not?
-                        #print 'ASCII formats are not that great for '
-                        #print 'blockwise reading. Consider storing this'
-                        #print 'data as a pts NetCDF format'
+                        # print 'ASCII formats are not that great for '
+                        # print 'blockwise reading. Consider storing this'
+                        # print 'data as a pts NetCDF format'
 
     ##
@@ -203 +204 @@
 
     ##
-    # @brief Save points data in instance.
-    # @param data_points Points data to store in instance and check.
+    # @brief Check data points.
+    # @param data_points Points data to check and store in instance.
     # @note Throws ValueError exception if no data.
     def check_data_points(self, data_points):
-        """Checks data points
-        """
+        """Checks data points"""
 
         if data_points is None:
@@ -214 +214 @@
             msg = 'There is no data or file provided!'
             raise ValueError, msg
-
         else:
             self.data_points = ensure_numeric(data_points)
@@ -226 +225 @@
     # @note Throws exception if unable to convert dict keys to numeric.
     def set_attributes(self, attributes):
-        """Check and assign attributes dictionary
-        """
+        """Check and assign attributes dictionary"""
 
         if attributes is None:
@@ -234 +232 @@
 
         if not isinstance(attributes, DictType):
-            #Convert single attribute into dictionary
+            # Convert single attribute into dictionary
             attributes = {DEFAULT_ATTRIBUTE: attributes}
 
-        #Check input attributes
+        # Check input attributes
         for key in attributes.keys():
             try:
                 attributes[key] = ensure_numeric(attributes[key])
             except:
-                msg = 'Attribute %s could not be converted' %key
-                msg += 'to a numeric vector'
-                raise msg
+                msg = ("Attribute '%s' (%s) could not be converted to a"
+                       "numeric vector" % (str(key), str(attributes[key])))
+                raise Exception, msg
 
         self.attributes = attributes
 
-    #def set_geo_reference(self, geo_reference):
-    #    # FIXME (Ole): Backwards compatibility - deprecate
-    #    self.setgeo_reference(geo_reference)
-
     ##
     # @brief Set the georeference of geospatial data.
-    # @param geo_reference The georeference data    to set.
+    # @param geo_reference The georeference data to set.
     # @note Will raise exception if param not instance of Geo_reference.
     def set_geo_reference(self, geo_reference):
@@ -263 +257 @@
         """
 
-        from anuga.coordinate_transforms.geo_reference import Geo_reference
-
         if geo_reference is None:
             # Use default - points are in absolute coordinates
@@ -275 +267 @@
             # FIXME (Ole): This exception will be raised even
             # if geo_reference is None. Is that the intent Duncan?
-            msg = 'Argument geo_reference must be a valid Geo_reference '
-            msg += 'object or None.'
-            raise msg
+            msg = ('Argument geo_reference must be a valid Geo_reference '
+                   'object or None.')
+            raise Expection, msg
 
         # If a geo_reference already exists, change the point data according to
@@ -384 +376 @@
     # @param isSouthHemisphere If True, return lat/lon points in S.Hemi.
     # @return A set of data points, in appropriate form.
-    def get_data_points(self, absolute=True, geo_reference=None,
-                        as_lat_long=False, isSouthHemisphere=True):
+    def get_data_points(self,
+                        absolute=True,
+                        geo_reference=None,
+                        as_lat_long=False,
+                        isSouthHemisphere=True):
         """Get coordinates for all data points as an Nx2 array
 
@@ -466 +461 @@
     # @return The new geospatial object.
     def __add__(self, other):
-        """Returns the addition of 2 geospatical objects,
+        """Returns the addition of 2 geospatial objects,
         objects are concatencated to the end of each other
 
@@ -494 +489 @@
             if self.attributes is None:
                 if other.attributes is not None:
-                    msg = 'Geospatial data must have the same '
-                    msg += 'attributes to allow addition.'
+                    msg = ('Geospatial data must have the same '
+                           'attributes to allow addition.')
                     raise Exception, msg
 
@@ -505 +500 @@
                         attrib1 = self.attributes[x]
                         attrib2 = other.attributes[x]
-                        new_attributes[x] = num.concatenate((attrib1, attrib2))
+                        new_attributes[x] = num.concatenate((attrib1, attrib2),
+                                                            axis=0) #??default#
                     else:
-                        msg = 'Geospatial data must have the same \n'
-                        msg += 'attributes to allow addition.'
+                        msg = ('Geospatial data must have the same '
+                               'attributes to allow addition.')
                         raise Exception, msg
         else:
-            #other is None:
+            # other is None:
             new_points = self.get_data_points(absolute=True)
             new_attributes = self.attributes
@@ -525 +521 @@
     # @return The new geospatial object.
     def __radd__(self, other):
-        """Handle cases like None + Geospatial_data(...)
-        """
+        """Handle cases like None + Geospatial_data(...)"""
 
         return self + other
 
-    ###
-    #  IMPORT/EXPORT POINTS FILES
-    ###
+################################################################################
+#  IMPORT/EXPORT POINTS FILES
+################################################################################
 
     ##
@@ -542 +537 @@
     def import_points_file(self, file_name, delimiter=None, verbose=False):
         """ load an .txt, .csv or .pts file
+
         Note: will throw an IOError/SyntaxError if it can't load the file.
         Catch these!
@@ -553 +549 @@
 
         attributes = {}
-        if file_name[-4:]== ".pts":
+        if file_name[-4:] == ".pts":
             try:
                 data_points, attributes, geo_reference = \
@@ -560 +556 @@
                 msg = 'Could not open file %s ' % file_name
                 raise IOError, msg
-        elif file_name[-4:]== ".txt" or file_name[-4:]== ".csv":
+        elif file_name[-4:] == ".txt" or file_name[-4:]== ".csv":
             try:
                 data_points, attributes, geo_reference = \
@@ -566 +562 @@
             except IOError, e:
                 # This should only be if a file is not found
-                msg = 'Could not open file %s. ' % file_name
-                msg += 'Check the file location.'
+                msg = ('Could not open file %s. Check the file location.'
+                       % file_name)
                 raise IOError, msg
             except SyntaxError, e:
                 # This should only be if there is a format error
-                msg = 'Could not open file %s. \n' %file_name
-                msg += Error_message['IOError']
+                msg = ('Problem with format of file %s.\n%s'
+                       % (file_name, Error_message['IOError']))
                 raise SyntaxError, msg
         else:
-            msg = 'Extension %s is unknown' %file_name[-4:]
+            msg = 'Extension %s is unknown' % file_name[-4:]
             raise IOError, msg
 
@@ -590 +586 @@
     def export_points_file(self, file_name, absolute=True,
                            as_lat_long=False, isSouthHemisphere=True):
-
-        """
-        write a points file, file_name, as a text (.csv) or binary (.pts) file
+        """write a points file as a text (.csv) or binary (.pts) file
+
         file_name is the file name, including the extension
         The point_dict is defined at the top of this file.
@@ -621 +616 @@
                                 self.get_all_attributes(),
                                 self.get_geo_reference())
-
         elif file_name[-4:] == ".txt" or file_name[-4:] == ".csv":
             msg = "ERROR: trying to write a .txt file with relative data."
@@ -628 +622 @@
                             self.get_data_points(absolute=True,
                                                  as_lat_long=as_lat_long,
-                                  isSouthHemisphere=isSouthHemisphere),
+                                           isSouthHemisphere=isSouthHemisphere),
                             self.get_all_attributes(),
                             as_lat_long=as_lat_long)
-
         elif file_name[-4:] == ".urs" :
             msg = "ERROR: Can not write a .urs file as a relative file."
@@ -637 +630 @@
             _write_urs_file(file_name,
                             self.get_data_points(as_lat_long=True,
-                                  isSouthHemisphere=isSouthHemisphere))
-
+                                           isSouthHemisphere=isSouthHemisphere))
         else:
             msg = 'Unknown file type %s ' %file_name
@@ -659 +651 @@
         """
 
-        #FIXME: add the geo_reference to this
+        # FIXME: add the geo_reference to this
         points = self.get_data_points()
-        sampled_points = num.take(points, indices)
+        sampled_points = num.take(points, indices, axis=0)
 
         attributes = self.get_all_attributes()
@@ -668 +660 @@
         if attributes is not None:
             for key, att in attributes.items():
-                sampled_attributes[key] = num.take(att, indices)
+                sampled_attributes[key] = num.take(att, indices, axis=0)
 
         return Geospatial_data(sampled_points, sampled_attributes)
@@ -679 +671 @@
     # @note Points in each result object are selected randomly.
     def split(self, factor=0.5, seed_num=None, verbose=False):
-        """Returns two
-        geospatial_data object, first is the size of the 'factor'
+        """Returns two geospatial_data object, first is the size of the 'factor'
         smaller the original and the second is the remainder. The two
-        new object are disjoin set of each other.
+        new objects are disjoint sets of each other.
 
         Points of the two new object have selected RANDOMLY.
@@ -698 +689 @@
         """
 
-        i=0
+        i = 0
         self_size = len(self)
         random_list = []
         remainder_list = []
-        new_size = round(factor*self_size)
+        new_size = round(factor * self_size)
 
         # Find unique random numbers
         if verbose: print "make unique random number list and get indices"
 
-        total=num.array(range(self_size), num.Int)     #array default#
+        total = num.array(range(self_size), num.int)    #array default#
         total_list = total.tolist()
 
@@ -721 +712 @@
         # Set seed if provided, mainly important for unit test!
         # plus recalcule seed when no seed provided.
-        if seed_num != None:
-            seed(seed_num, seed_num)
+        if seed_num is not None:
+            seed(seed_num)
         else:
             seed()
 
-        if verbose: print "seed:", get_seed()
+        #if verbose: print "seed:", get_seed()
 
         random_num = randint(0, self_size-1, (int(new_size),))
         random_num = random_num.tolist()
 
-        #need to sort and reverse so the pop() works correctly
+        # need to sort and reverse so the pop() works correctly
         random_num.sort()
         random_num.reverse()
@@ -737 +728 @@
         if verbose: print "make random number list and get indices"
 
-        j=0
-        k=1
+        j = 0
+        k = 1
         remainder_list = total_list[:]
 
         # pops array index (random_num) from remainder_list
-        # (which starts as the
-        # total_list and appends to random_list
+        # (which starts as the total_list and appends to random_list)
         random_num_len = len(random_num)
         for i in random_num:
             random_list.append(remainder_list.pop(i))
             j += 1
-            #prints progress
+            # prints progress
             if verbose and round(random_num_len/10*k) == j:
                 print '(%s/%s)' % (j, random_num_len)
@@ -754 +744 @@
 
         # FIXME: move to tests, it might take a long time
-        # then create an array of random lenght between 500 and 1000,
+        # then create an array of random length between 500 and 1000,
         # and use a random factor between 0 and 1
         # setup for assertion
@@ -760 +750 @@
         test_total.extend(remainder_list)
         test_total.sort()
-        msg = 'The two random lists made from the original list when added ' \
-              'together DO NOT equal the original list'
+        msg = ('The two random lists made from the original list when added '
+               'together DO NOT equal the original list')
         assert total_list == test_total, msg
 
@@ -778 +768 @@
         file pointer position
         """
-        from Scientific.IO.NetCDF import NetCDFFile
-        #FIXME - what to do if the file isn't there
+
+        # FIXME - what to do if the file isn't there
 
         # FIXME (Ole): Shouldn't this go into the constructor?
@@ -807 +797 @@
             self.number_of_blocks = self.number_of_points/self.max_read_lines
             # This computes the number of full blocks. The last block may be
-            # smaller and won't be ircluded in this estimate.
+            # smaller and won't be included in this estimate.
 
             if self.verbose is True:
-                print 'Reading %d points (in ~%d blocks) from file %s. ' \
-                      % (self.number_of_points,
-                         self.number_of_blocks,
-                         self.file_name),
-                print 'Each block consists of %d data points' \
-                      % self.max_read_lines
-
+                print ('Reading %d points (in ~%d blocks) from file %s. '
+                       % (self.number_of_points, self.number_of_blocks,
+                          self.file_name)),
+                print ('Each block consists of %d data points'
+                       % self.max_read_lines)
         else:
             # Assume the file is a csv file
@@ -847 +835 @@
 
             if self.verbose is True:
-                if self.show_verbose >= self.start_row \
-                   and self.show_verbose < fin_row:
-                    print 'Reading block %d (points %d to %d) out of %d'\
-                          %(self.block_number,
-                            self.start_row,
-                            fin_row,
-                            self.number_of_blocks)
+                if (self.show_verbose >= self.start_row
+                    and self.show_verbose < fin_row):
+                    print ('Reading block %d (points %d to %d) out of %d'
+                           % (self.block_number, self.start_row,
+                              fin_row, self.number_of_blocks))
 
                     self.show_verbose += max(self.max_read_lines,
                                              self.verbose_block_size)
-
 
             # Read next block
@@ -869 +854 @@
 
             self.block_number += 1
-
         else:
             # Assume the file is a csv file
@@ -876 +860 @@
                  att_dict,
                  geo_ref,
-                 self.file_pointer) = \
-                        _read_csv_file_blocking(self.file_pointer,
-                                                self.header[:],
-                                                max_read_lines=\
-                                                    self.max_read_lines,
-                                                verbose=self.verbose)
+                 self.file_pointer) = _read_csv_file_blocking(self.file_pointer,
+                                                              self.header[:],
+                                                              max_read_lines=
+                                                           self.max_read_lines,
+                                                              verbose=
+                                                                  self.verbose)
 
                 # Check that the zones haven't changed.
@@ -888 +872 @@
                     self.blocking_georef = geo_ref
                 elif self.blocking_georef is not None:
-                    msg = 'Geo reference given, then not given.'
-                    msg += ' This should not happen.'
+                    msg = ('Geo reference given, then not given.'
+                           ' This should not happen.')
                     raise ValueError, msg
                 geo = Geospatial_data(pointlist, att_dict, geo_ref)
@@ -907 +891 @@
                 del self.file_pointer
                 # This should only be if there is a format error
-                msg = 'Could not open file %s. \n' % self.file_name
-                msg += Error_message['IOError']
+                msg = ('Could not open file %s.\n%s'
+                       % (self.file_name, Error_message['IOError']))
                 raise SyntaxError, msg
         return geo
-
 
 ##################### Error messages ###########
 Error_message = {}
 Em = Error_message
-Em['IOError'] = "NOTE: The format for a comma separated .txt/.csv file is:\n"
-Em['IOError'] += "        1st line:     [column names]\n"
-Em['IOError'] += "        other lines:  [x value], [y value], [attributes]\n"
-Em['IOError'] += "\n"
-Em['IOError'] += "           for example:\n"
-Em['IOError'] += "           x, y, elevation, friction\n"
-Em['IOError'] += "           0.6, 0.7, 4.9, 0.3\n"
-Em['IOError'] += "           1.9, 2.8, 5, 0.3\n"
-Em['IOError'] += "           2.7, 2.4, 5.2, 0.3\n"
-Em['IOError'] += "\n"
-Em['IOError'] += "The first two columns are assumed to be x, y coordinates.\n"
-Em['IOError'] += "The attribute values must be numeric.\n"
+Em['IOError'] = ('NOTE: The format for a comma separated .txt/.csv file is:\n'
+                 '        1st line:     [column names]\n'
+                 '        other lines:  [x value], [y value], [attributes]\n'
+                 '\n'
+                 '           for example:\n'
+                 '           x, y, elevation, friction\n'
+                 '           0.6, 0.7, 4.9, 0.3\n'
+                 '           1.9, 2.8, 5, 0.3\n'
+                 '           2.7, 2.4, 5.2, 0.3\n'
+                 '\n'
+                 'The first two columns are assumed to be x, y coordinates.\n'
+                 'The attribute values must be numeric.\n')
 
 ##
@@ -936 +919 @@
 # @param data_points ??
 # @param points_are_lats_longs ??
+# @note IS THIS USED???
 def _set_using_lat_long(latitudes,
                         longitudes,
@@ -941 +925 @@
                         data_points,
                         points_are_lats_longs):
-    """
-    if the points has lat long info, assume it is in (lat, long) order.
-    """
+    """If the points has lat long info, assume it is in (lat, long) order."""
 
     if geo_reference is not None:
-        msg = "A georeference is specified yet latitude and longitude " \
-              "are also specified!"
+        msg = ('A georeference is specified yet latitude and longitude '
+               'are also specified!')
         raise ValueError, msg
 
     if data_points is not None and not points_are_lats_longs:
-        msg = "Data points are specified yet latitude and longitude are " \
-              "also specified."
+        msg = ('Data points are specified yet latitude and longitude are '
+               'also specified.')
         raise ValueError, msg
 
@@ -988 +970 @@
     """Read .pts NetCDF file
 
-    Return a dic of array of points, and dic of array of attribute
+    Return a (dict_points, dict_attribute, geo_ref)
     eg
-    dic['points'] = [[1.0,2.0],[3.0,5.0]]
-    dic['attributelist']['elevation'] = [[7.0,5.0]]
+    dict['points'] = [[1.0,2.0],[3.0,5.0]]
+    dict['attributelist']['elevation'] = [[7.0,5.0]]
     """
-
-    from Scientific.IO.NetCDF import NetCDFFile
 
     if verbose: print 'Reading ', file_name
@@ -1055 +1035 @@
                                                  header,
                                                  max_read_lines=1e30)
-                                    #If the file is bigger that this, block..
+                                    # If the file is bigger that this, block..
                                     # FIXME (Ole) What's up here?
     except ANUGAError:
@@ -1095 +1075 @@
 # @param verbose True if this function is to be verbose.
 # @note Will throw IndexError, SyntaxError exceptions.
-def _read_csv_file_blocking(file_pointer, header,
+def _read_csv_file_blocking(file_pointer,
+                            header,
                             delimiter=CSV_DELIMITER,
                             max_read_lines=MAX_READ_LINES,
@@ -1136 +1117 @@
             if len(header) != len(numbers):
                 file_pointer.close()
-                msg = "File load error. " \
-                      "There might be a problem with the file header."
+                msg = ('File load error. '
+                       'There might be a problem with the file header.')
                 raise SyntaxError, msg
             for i,n in enumerate(numbers):
                 n.strip()
                 if n != '\n' and n != '':
-                    #attributes.append(float(n))
                     att_dict.setdefault(header[i],[]).append(float(n))
         except ValueError:
@@ -1150 +1130 @@
         raise StopIteration
 
-    pointlist = num.array(points, num.Float)
+    pointlist = num.array(points, num.float)
     for key in att_dict.keys():
-        att_dict[key] = num.array(att_dict[key], num.Float)
+        att_dict[key] = num.array(att_dict[key], num.float)
 
     # Do stuff here so the info is in lat's and longs
@@ -1183 +1163 @@
 # @note Will throw IOError and AttributeError exceptions.
 def _read_pts_file_header(fid, verbose=False):
-    """Read the geo_reference and number_of_points from a .pts file
-    """
+    '''Read the geo_reference and number_of_points from a .pts file'''
 
     keys = fid.variables.keys()
@@ -1205 +1184 @@
 
 ##
-# @brief Read the body of a .csf file, with blocking.
+# @brief Read the body of a .pts file, with blocking.
 # @param fid Handle to already open file.
 # @param start_row Start row index of points to return.
@@ -1212 +1191 @@
 # @return Tuple of (pointlist, attributes).
 def _read_pts_file_blocking(fid, start_row, fin_row, keys):
-    """Read the body of a .csv file.
-    """
+    '''Read the body of a .pts file.'''
 
     pointlist = num.array(fid.variables['points'][start_row:fin_row])
@@ -1239 +1217 @@
 
     WARNING: This function mangles the point_atts data structure
-    #F??ME: (DSG)This format has issues.
+    # F??ME: (DSG)This format has issues.
     # There can't be an attribute called points
     # consider format change
@@ -1251 +1229 @@
     """
 
-    from Scientific.IO.NetCDF import NetCDFFile
-
     # NetCDF file definition
     outfile = NetCDFFile(file_name, netcdf_mode_w)
@@ -1258 +1234 @@
     # Create new file
     outfile.institution = 'Geoscience Australia'
-    outfile.description = 'NetCDF format for compact and portable storage ' \
-                          'of spatial point data'
+    outfile.description = ('NetCDF format for compact and portable storage '
+                           'of spatial point data')
 
     # Dimension definitions
     shape = write_data_points.shape[0]
     outfile.createDimension('number_of_points', shape)
-    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
+    outfile.createDimension('number_of_dimensions', 2) # This is 2d data
 
     # Variable definition
-    outfile.createVariable('points', num.Float, ('number_of_points',
-                                                 'number_of_dimensions'))
-
-    #create variables
-    outfile.variables['points'][:] = write_data_points #.astype(Float32)
+    outfile.createVariable('points', netcdf_float,
+                           ('number_of_points', 'number_of_dimensions'))
+
+    # create variables
+    outfile.variables['points'][:] = write_data_points
 
     if write_attributes is not None:
         for key in write_attributes.keys():
-            outfile.createVariable(key, num.Float, ('number_of_points',))
-            outfile.variables[key][:] = write_attributes[key] #.astype(Float32)
+            outfile.createVariable(key, netcdf_float, ('number_of_points',))
+            outfile.variables[key][:] = write_attributes[key]
 
     if write_geo_reference is not None:
@@ -1296 +1272 @@
                     as_lat_long=False,
                     delimiter=','):
-    """Write a .csv file.
-    """
+    """Write a .csv file."""
 
     points = write_data_points
@@ -1361 +1336 @@
 # @return ??
 def _point_atts2array(point_atts):
-    point_atts['pointlist'] = num.array(point_atts['pointlist'], num.Float)
+    point_atts['pointlist'] = num.array(point_atts['pointlist'], num.float)
 
     for key in point_atts['attributelist'].keys():
         point_atts['attributelist'][key] = \
-                num.array(point_atts['attributelist'][key], num.Float)
+                num.array(point_atts['attributelist'][key], num.float)
 
     return point_atts
@@ -1375 +1350 @@
 # @return A points dictionary.
 def geospatial_data2points_dictionary(geospatial_data):
-    """Convert geospatial data to points_dictionary
-    """
+    """Convert geospatial data to points_dictionary"""
 
     points_dictionary = {}
@@ -1396 +1370 @@
 # @param points_dictionary A points dictionary to convert.
 def points_dictionary2geospatial_data(points_dictionary):
-    """Convert points_dictionary to geospatial data object
-    """
-
-    msg = 'Points dictionary must have key pointlist'
+    """Convert points_dictionary to geospatial data object"""
+
+    msg = "Points dictionary must have key 'pointlist'"
     assert points_dictionary.has_key('pointlist'), msg
 
-    msg = 'Points dictionary must have key attributelist'
+    msg = "Points dictionary must have key 'attributelist'"
     assert points_dictionary.has_key('attributelist'), msg
 
@@ -1412 +1385 @@
     return Geospatial_data(points_dictionary['pointlist'],
                            points_dictionary['attributelist'],
-                           geo_reference = geo)
-
-
-##
-# @brief Split a string into 'clean' fields.
-# @param line The string to process.
-# @param delimiter The delimiter string to split 'line' with.
-# @return A list of 'cleaned' field strings.
-# @note Any fields that were initially zero length will be removed.
-def clean_line(line,delimiter):
-    """Remove whitespace
-    """
-
-    line = line.strip()                 # probably unnecessary RW
-    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
-
-    return numbers
+                           geo_reference=geo)
 
 
@@ -1460 +1409 @@
     """
 
-    import copy
     # Input check
     if isinstance(points, basestring):
-        #It's a string - assume it is a point file
+        # It's a string - assume it is a point file
         points = Geospatial_data(file_name=points)
 
@@ -1471 +1419 @@
         assert geo_reference == None, msg
     else:
-        points = ensure_numeric(copy.copy(points), num.Float)
+        points = ensure_numeric(copy.copy(points), num.float)
 
     # Sort of geo_reference and convert points
     if geo_reference is None:
-        geo = None #Geo_reference()
+        geo = None    # Geo_reference()
     else:
         if isinstance(geo_reference, Geo_reference):
@@ -1494 +1442 @@
 # @return A geospatial object.
 def ensure_geospatial(points, geo_reference=None):
-    """
-    This function inputs several formats and
-    outputs one format. - a geospatial_data instance.
+    """Convert various data formats to a geospatial_data instance.
 
     Inputed formats are;
@@ -1508 +1454 @@
     """
 
-    import copy
     # Input check
     if isinstance(points, Geospatial_data):
@@ -1516 +1461 @@
     else:
         # List or numeric array of absolute points
-        points = ensure_numeric(points, num.Float)
+        points = ensure_numeric(points, num.float)
 
     # Sort out geo reference
@@ -1565 +1510 @@
                                      cache=False,
                                      verbose=False):
-    """
-    Removes a small random sample of points from 'data_file'. Then creates
-    models with different alpha values from 'alpha_list' and cross validates
-    the predicted value to the previously removed point data. Returns the
-    alpha value which has the smallest covariance.
+    """Removes a small random sample of points from 'data_file'.
+    Then creates models with different alpha values from 'alpha_list' and
+    cross validates the predicted value to the previously removed point data.
+    Returns the alpha value which has the smallest covariance.
 
     data_file: must not contain points outside the boundaries defined
-               and it either a pts, txt or csv file.
+               and it must be either a pts, txt or csv file.
 
     alpha_list: the alpha values to test in a single list
@@ -1615 +1559 @@
     from anuga.fit_interpolate.benchmark_least_squares import mem_usage
 
-
-    attribute_smoothed='elevation'
+    attribute_smoothed = 'elevation'
 
     if mesh_file is None:
@@ -1622 +1565 @@
         mesh_file = 'temp.msh'
 
-        if north_boundary is None or south_boundary is None \
-           or east_boundary is None or west_boundary is None:
+        if (north_boundary is None or south_boundary is None
+            or east_boundary is None or west_boundary is None):
             no_boundary = True
         else:
@@ -1630 +1573 @@
         if no_boundary is True:
             msg = 'All boundaries must be defined'
-            raise msg
-
-        poly_topo = [[east_boundary,south_boundary],
-                     [east_boundary,north_boundary],
-                     [west_boundary,north_boundary],
-                     [west_boundary,south_boundary]]
+            raise Expection, msg
+
+        poly_topo = [[east_boundary, south_boundary],
+                     [east_boundary, north_boundary],
+                     [west_boundary, north_boundary],
+                     [west_boundary, south_boundary]]
 
         create_mesh_from_regions(poly_topo,
@@ -1647 +1590 @@
 
     else: # if mesh file provided
-        #test mesh file exists?
+        # test mesh file exists?
         if verbose: "reading from file: %s" % mesh_file
         if access(mesh_file,F_OK) == 0:
@@ -1653 +1596 @@
             raise IOError, msg
 
-    #split topo data
+    # split topo data
     if verbose: print 'Reading elevation file: %s' % data_file
     G = Geospatial_data(file_name = data_file)
@@ -1665 +1608 @@
     if alpha_list == None:
         alphas = [0.001,0.01,100]
-        #alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01,\
-         #         0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
+        #alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01,
+        #          0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
     else:
         alphas = alpha_list
 
-    #creates array with columns 1 and 2 are x, y. column 3 is elevation
-    #4 onwards is the elevation_predicted using the alpha, which will
-    #be compared later against the real removed data
-    data = num.array([], typecode=num.Float)
-
-    data=num.resize(data, (len(points), 3+len(alphas)))
-
-    #gets relative point from sample
+    # creates array with columns 1 and 2 are x, y. column 3 is elevation
+    # 4 onwards is the elevation_predicted using the alpha, which will
+    # be compared later against the real removed data
+    data = num.array([], dtype=num.float)
+
+    data = num.resize(data, (len(points), 3+len(alphas)))
+
+    # gets relative point from sample
     data[:,0] = points[:,0]
     data[:,1] = points[:,1]
@@ -1683 +1626 @@
     data[:,2] = elevation_sample
 
-    normal_cov=num.array(num.zeros([len(alphas), 2]), typecode=num.Float)
+    normal_cov = num.array(num.zeros([len(alphas), 2]), dtype=num.float)
 
     if verbose: print 'Setup computational domains with different alphas'
 
     for i, alpha in enumerate(alphas):
-        #add G_other data to domains with different alphas
+        # add G_other data to domains with different alphas
         if verbose:
-            print '\n Calculating domain and mesh for Alpha = ', alpha, '\n'
+            print '\nCalculating domain and mesh for Alpha =', alpha, '\n'
         domain = Domain(mesh_file, use_cache=cache, verbose=verbose)
         if verbose: print domain.statistics()
@@ -1702 +1645 @@
         points_geo = Geospatial_data(points, domain.geo_reference)
 
-        #returns the predicted elevation of the points that were "split" out
-        #of the original data set for one particular alpha
+        # returns the predicted elevation of the points that were "split" out
+        # of the original data set for one particular alpha
         if verbose: print 'Get predicted elevation for location to be compared'
         elevation_predicted = \
@@ -1709 +1652 @@
                     get_values(interpolation_points=points_geo)
 
-        #add predicted elevation to array that starts with x, y, z...
+        # add predicted elevation to array that starts with x, y, z...
         data[:,i+3] = elevation_predicted
 
         sample_cov = cov(elevation_sample)
         ele_cov = cov(elevation_sample - elevation_predicted)
-        normal_cov[i,:] = [alpha,ele_cov / sample_cov]
+        normal_cov[i,:] = [alpha, ele_cov / sample_cov]
 
         if verbose:
@@ -1722 +1665 @@
 
     normal_cov0 = normal_cov[:,0]
-    normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0))
+    normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0), axis=0)
 
     if plot_name is not None:
@@ -1737 +1680 @@
         print 'Final results:'
         for i, alpha in enumerate(alphas):
-            print 'covariance for alpha %s = %s ' \
-                  % (normal_cov[i][0], normal_cov[i][1])
-        print '\n Optimal alpha is: %s ' \
-              % normal_cov_new[(num.argmin(normal_cov_new, axis=0))[1], 0]
+            print ('covariance for alpha %s = %s '
+                   % (normal_cov[i][0], normal_cov[i][1]))
+        print ('\nOptimal alpha is: %s '
+               % normal_cov_new[(num.argmin(normal_cov_new, axis=0))[1], 0])
 
     # covariance and optimal alpha
@@ -1822 +1765 @@
     from anuga.fit_interpolate.benchmark_least_squares import mem_usage
 
-
     attribute_smoothed = 'elevation'
 
@@ -1828 +1770 @@
         mesh_file = 'temp.msh'
 
-        if north_boundary is None or south_boundary is None \
-           or east_boundary is None or west_boundary is None:
+        if (north_boundary is None or south_boundary is None
+            or east_boundary is None or west_boundary is None):
             no_boundary = True
         else:
@@ -1836 +1778 @@
         if no_boundary is True:
             msg = 'All boundaries must be defined'
-            raise msg
-
-        poly_topo = [[east_boundary,south_boundary],
-                     [east_boundary,north_boundary],
-                     [west_boundary,north_boundary],
-                     [west_boundary,south_boundary]]
+            raise Expection, msg
+
+        poly_topo = [[east_boundary, south_boundary],
+                     [east_boundary, north_boundary],
+                     [west_boundary, north_boundary],
+                     [west_boundary, south_boundary]]
 
         create_mesh_from_regions(poly_topo,
@@ -1853 +1795 @@
 
     else: # if mesh file provided
-        #test mesh file exists?
+        # test mesh file exists?
         if access(mesh_file,F_OK) == 0:
             msg = "file %s doesn't exist!" % mesh_file
             raise IOError, msg
 
-    #split topo data
+    # split topo data
     G = Geospatial_data(file_name=data_file)
     if verbose: print 'start split'
@@ -1865 +1807 @@
     points = G_small.get_data_points()
 
-    #FIXME: Remove points outside boundary polygon
-#    print 'new point',len(points)
-#
-#    new_points=[]
-#    new_points=array([],typecode=Float)
-#    new_points=resize(new_points,(len(points),2))
-#    print "BOUNDARY", boundary_poly
-#    for i,point in enumerate(points):
-#        if is_inside_polygon(point,boundary_poly, verbose=True):
-#            new_points[i] = point
-#            print"WOW",i,new_points[i]
-#    points = new_points
-
     if verbose: print "Number of points in sample to compare: ", len(points)
 
     if alpha_list == None:
         alphas = [0.001,0.01,100]
-        #alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01,\
-         #         0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
+        #alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01,
+        #          0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
     else:
         alphas = alpha_list
@@ -1892 +1821 @@
 
     for alpha in alphas:
-        #add G_other data to domains with different alphas
+        # add G_other data to domains with different alphas
         if verbose:
-            print '\n Calculating domain and mesh for Alpha = ', alpha, '\n'
+            print '\nCalculating domain and mesh for Alpha =', alpha, '\n'
         domain = Domain(mesh_file, use_cache=cache, verbose=verbose)
         if verbose: print domain.statistics()
@@ -1904 +1833 @@
         domains[alpha] = domain
 
-    #creates array with columns 1 and 2 are x, y. column 3 is elevation
-    #4 onwards is the elevation_predicted using the alpha, which will
-    #be compared later against the real removed data
-    data = num.array([], typecode=num.Float)
+    # creates array with columns 1 and 2 are x, y. column 3 is elevation
+    # 4 onwards is the elevation_predicted using the alpha, which will
+    # be compared later against the real removed data
+    data = num.array([], dtype=num.float)
 
     data = num.resize(data, (len(points), 3+len(alphas)))
 
-    #gets relative point from sample
+    # gets relative point from sample
     data[:,0] = points[:,0]
     data[:,1] = points[:,1]
@@ -1917 +1846 @@
     data[:,2] = elevation_sample
 
-    normal_cov = num.array(num.zeros([len(alphas), 2]), typecode=num.Float)
+    normal_cov = num.array(num.zeros([len(alphas), 2]), dtype=num.float)
 
     if verbose:
         print 'Determine difference between predicted results and actual data'
+
     for i, alpha in enumerate(domains):
         if verbose: print'Alpha =', alpha
 
         points_geo = domains[alpha].geo_reference.change_points_geo_ref(points)
-        #returns the predicted elevation of the points that were "split" out
-        #of the original data set for one particular alpha
+        # returns the predicted elevation of the points that were "split" out
+        # of the original data set for one particular alpha
         elevation_predicted = \
                 domains[alpha].quantities[attribute_smoothed].\
                         get_values(interpolation_points=points_geo)
 
-        #add predicted elevation to array that starts with x, y, z...
+        # add predicted elevation to array that starts with x, y, z...
         data[:,i+3] = elevation_predicted
 
@@ -1941 +1871 @@
 
     normal_cov0 = normal_cov[:,0]
-    normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0))
+    normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0), axis=0)
 
     if plot_name is not None:

anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

@@ -1 +1 @@
 #!/usr/bin/env python
-
 
 import unittest
 import os
+import tempfile
 from math import sqrt, pi
-import tempfile
-
-import Numeric as num
+
+import numpy as num
 
 from anuga.geospatial_data.geospatial_data import *
@@ -15 +14 @@
 from anuga.utilities.system_tools import get_host_name
 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
+from anuga.config import netcdf_float
+
 
 class Test_Geospatial_data(unittest.TestCase):
@@ -23 +24 @@
         pass
 
-
     def test_0(self):
         #Basic points
         from anuga.coordinate_transforms.geo_reference import Geo_reference
-        
+
         points = [[1.0, 2.1], [3.0, 5.3]]
         G = Geospatial_data(points)
-
         assert num.allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]])
 
@@ -37 +36 @@
         rep = `G`
         ref = '[[ 1.   2.1]\n [ 3.   5.3]]'
-
-        msg = 'Representation %s is not equal to %s' %(rep, ref)
+        msg = 'Representation %s is not equal to %s' % (rep, ref)
         assert rep == ref, msg
 
         #Check getter
         assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3]])
-        
+
         #Check defaults
         assert G.attributes is None
-        
         assert G.geo_reference.zone == Geo_reference().zone
         assert G.geo_reference.xllcorner == Geo_reference().xllcorner
         assert G.geo_reference.yllcorner == Geo_reference().yllcorner
-        
 
     def test_1(self):
         points = [[1.0, 2.1], [3.0, 5.3]]
         attributes = [2, 4]
-        G = Geospatial_data(points, attributes)       
+        G = Geospatial_data(points, attributes)
         assert G.attributes.keys()[0] == DEFAULT_ATTRIBUTE
         assert num.allclose(G.attributes.values()[0], [2, 4])
-        
 
     def test_2(self):
         from anuga.coordinate_transforms.geo_reference import Geo_reference
+
         points = [[1.0, 2.1], [3.0, 5.3]]
         attributes = [2, 4]
@@ -71 +67 @@
         assert G.geo_reference.yllcorner == 200
 
-
     def test_get_attributes_1(self):
         from anuga.coordinate_transforms.geo_reference import Geo_reference
+
         points = [[1.0, 2.1], [3.0, 5.3]]
         attributes = [2, 4]
@@ -79 +75 @@
                             geo_reference=Geo_reference(56, 100, 200))
 
-
         P = G.get_data_points(absolute=False)
-        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])        
+        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])
 
         P = G.get_data_points(absolute=True)
-        assert num.allclose(P, [[101.0, 202.1], [103.0, 205.3]])        
+        assert num.allclose(P, [[101.0, 202.1], [103.0, 205.3]])
 
         V = G.get_attributes() #Simply get them
@@ -94 +89 @@
     def test_get_attributes_2(self):
         #Multiple attributes
-        
-        
         from anuga.coordinate_transforms.geo_reference import Geo_reference
+
         points = [[1.0, 2.1], [3.0, 5.3]]
         attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
@@ -103 +97 @@
                             default_attribute_name='a1')
 
-
         P = G.get_data_points(absolute=False)
-        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])        
-        
+        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])
+
         V = G.get_attributes() #Get default attribute
         assert num.allclose(V, [2, 4])
@@ -119 +112 @@
         assert num.allclose(V, [79.4, -7])
 
+                # FIXME: use failUnlessRaises()
         try:
             V = G.get_attributes('hdnoatedu') #Invalid
@@ -124 +118 @@
             pass
         else:
-            raise 'Should have raised exception' 
+            raise Exception, 'Should have raised exception'
 
     def test_get_data_points(self):
-        points_ab = [[12.5,34.7],[-4.5,-60.0]]
+        points_ab = [[12.5, 34.7], [-4.5, -60.0]]
         x_p = -10
         y_p = -40
         geo_ref = Geo_reference(56, x_p, y_p)
         points_rel = geo_ref.change_points_geo_ref(points_ab)
-        
+
         spatial = Geospatial_data(points_rel, geo_reference=geo_ref)
-
         results = spatial.get_data_points(absolute=False)
-        
         assert num.allclose(results, points_rel)
-        
+
         x_p = -1770
         y_p = 4.01
         geo_ref = Geo_reference(56, x_p, y_p)
         points_rel = geo_ref.change_points_geo_ref(points_ab)
-        results = spatial.get_data_points \
-                  ( geo_reference=geo_ref)
-        
+        results = spatial.get_data_points(geo_reference=geo_ref)
+
         assert num.allclose(results, points_rel)
 
-  
     def test_get_data_points_lat_long(self):
-        # lat long [-30.],[130]
-        #Zone:   52    
-        #Easting:  596450.153  Northing: 6680793.777 
-        # lat long [-32.],[131]
-        #Zone:   52    
-        #Easting:  688927.638  Northing: 6457816.509 
-        
-        points_Lat_long = [[-30.,130], [-32,131]]
-        
-        spatial = Geospatial_data(latitudes=[-30, -32.],
-                                  longitudes=[130, 131])
-
+        # lat long [-30.], [130]
+        # Zone:   52
+        # Easting:  596450.153  Northing: 6680793.777
+        # lat long [-32.], [131]
+        # Zone:   52
+        # Easting:  688927.638  Northing: 6457816.509
+
+        points_Lat_long = [[-30., 130], [-32, 131]]
+
+        spatial = Geospatial_data(latitudes=[-30, -32.], longitudes=[130, 131])
         results = spatial.get_data_points(as_lat_long=True)
-        #print "test_get_data_points_lat_long - results", results
-        #print "points_Lat_long",points_Lat_long 
         assert num.allclose(results, points_Lat_long)
-      
+
     def test_get_data_points_lat_longII(self):
         # x,y  North,east long,lat
         boundary_polygon = [[ 250000, 7630000]]
         zone = 50
-        
+
         geo_reference = Geo_reference(zone=zone)
-        geo = Geospatial_data(boundary_polygon,geo_reference=geo_reference)
+        geo = Geospatial_data(boundary_polygon ,geo_reference=geo_reference)
         seg_lat_long = geo.get_data_points(as_lat_long=True)
-        lat_result = degminsec2decimal_degrees(-21,24,54)
-        long_result = degminsec2decimal_degrees(114,35,17.89)
-        #print "seg_lat_long", seg_lat_long [0][0]
-        #print "lat_result",lat_result 
+        lat_result = degminsec2decimal_degrees(-21, 24, 54)
+        long_result = degminsec2decimal_degrees(114, 35, 17.89)
+        assert num.allclose(seg_lat_long[0][0], lat_result)     #lat
+        assert num.allclose(seg_lat_long[0][1], long_result)    #long
+
+    def test_get_data_points_lat_longIII(self):
+        # x,y  North,east long,lat
+        # for northern hemisphere
+        boundary_polygon = [[419944.8, 918642.4]]
+        zone = 47
+
+        geo_reference = Geo_reference(zone=zone)
+        geo = Geospatial_data(boundary_polygon, geo_reference=geo_reference)
+        seg_lat_long = geo.get_data_points(as_lat_long=True,
+                                           isSouthHemisphere=False)
+        lat_result = degminsec2decimal_degrees(8.31, 0, 0)
+        long_result = degminsec2decimal_degrees(98.273, 0, 0)
         assert num.allclose(seg_lat_long[0][0], lat_result)#lat
         assert num.allclose(seg_lat_long[0][1], long_result)#long
 
-
-    def test_get_data_points_lat_longIII(self):
-        # x,y  North,east long,lat
-        #for northern hemisphere
-        boundary_polygon = [[419944.8, 918642.4]]
-        zone = 47
-        
-        geo_reference = Geo_reference(zone=zone)
-        geo = Geospatial_data(boundary_polygon,
-                              geo_reference=geo_reference)
-                              
-        seg_lat_long = geo.get_data_points(as_lat_long=True,
-                                           isSouthHemisphere=False)
-                                           
-        lat_result = degminsec2decimal_degrees(8.31,0,0)
-        long_result = degminsec2decimal_degrees(98.273,0,0)
-        #print "seg_lat_long", seg_lat_long [0]
-        #print "lat_result",lat_result 
-        assert num.allclose(seg_lat_long[0][0], lat_result)#lat
-        assert num.allclose(seg_lat_long[0][1], long_result)#long
-
-
-              
     def test_set_geo_reference(self):
-        """test_set_georeference
-        
-        Test that georeference can be changed without changing the 
+        '''test_set_georeference
+
+        Test that georeference can be changed without changing the
         absolute values.
-        """
-            
-        points_ab = [[12.5,34.7],[-4.5,-60.0]]
+        '''
+
+        points_ab = [[12.5, 34.7], [-4.5, -60.0]]
         x_p = -10
         y_p = -40
         geo_ref = Geo_reference(56, x_p, y_p)
         points_rel = geo_ref.change_points_geo_ref(points_ab)
-        
+
         # Create without geo_ref properly set
-        G = Geospatial_data(points_rel)        
+        G = Geospatial_data(points_rel)
         assert not num.allclose(points_ab, G.get_data_points(absolute=True))
-        
+
         # Create the way it should be
         G = Geospatial_data(points_rel, geo_reference=geo_ref)
         assert num.allclose(points_ab, G.get_data_points(absolute=True))
-        
+
         # Change georeference and check that absolute values are unchanged.
         x_p = 10
@@ -231 +207 @@
         new_geo_ref = Geo_reference(56, x_p, y_p)
         G.set_geo_reference(new_geo_ref)
-        assert num.allclose(points_ab, G.get_data_points(absolute=True))
-        
-
-                
-        
+        msg = ('points_ab=\n%s\nbut G.get_data_points(absolute=True)=\n%s'
+               % (str(points_ab), str(G.get_data_points(absolute=True))))
+        assert num.allclose(points_ab, G.get_data_points(absolute=True)), msg
+
     def test_conversions_to_points_dict(self):
         #test conversions to points_dict
-        
-        
         from anuga.coordinate_transforms.geo_reference import Geo_reference
+
         points = [[1.0, 2.1], [3.0, 5.3]]
         attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
@@ -247 +221 @@
                             default_attribute_name='a1')
 
-
         points_dict = geospatial_data2points_dictionary(G)
 
         assert points_dict.has_key('pointlist')
-        assert points_dict.has_key('attributelist')        
+        assert points_dict.has_key('attributelist')
         assert points_dict.has_key('geo_reference')
 
@@ -259 +232 @@
         assert A.has_key('a0')
         assert A.has_key('a1')
-        assert A.has_key('a2')        
+        assert A.has_key('a2')
 
         assert num.allclose( A['a0'], [0, 0] )
-        assert num.allclose( A['a1'], [2, 4] )        
+        assert num.allclose( A['a1'], [2, 4] )
         assert num.allclose( A['a2'], [79.4, -7] )
-
 
         geo = points_dict['geo_reference']
         assert geo is G.geo_reference
 
-
     def test_conversions_from_points_dict(self):
-        """test conversions from points_dict
-        """
+        '''test conversions from points_dict'''
 
         from anuga.coordinate_transforms.geo_reference import Geo_reference
-        
+
         points = [[1.0, 2.1], [3.0, 5.3]]
         attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
@@ -283 +253 @@
         points_dict['attributelist'] = attributes
         points_dict['geo_reference'] = Geo_reference(56, 100, 200)
-        
 
         G = points_dictionary2geospatial_data(points_dict)
-
         P = G.get_data_points(absolute=False)
-        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])        
-        
-        #V = G.get_attribute_values() #Get default attribute
-        #assert allclose(V, [2, 4])
+        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])
 
         V = G.get_attributes('a0') #Get by name
@@ -303 +268 @@
 
     def test_add(self):
-        """ test the addition of two geospatical objects
-            no geo_reference see next test
-        """
+        '''test the addition of two geospatical objects
+        no geo_reference see next test
+        '''
+
         points = [[1.0, 2.1], [3.0, 5.3]]
-        attributes = {'depth':[2, 4], 'elevation':[6.1, 5]}
-        attributes1 = {'depth':[2, 4], 'elevation':[2.5, 1]}
-        G1 = Geospatial_data(points, attributes)        
-        G2 = Geospatial_data(points, attributes1) 
-        
-#        g3 = geospatial_data2points_dictionary(G1)
-#        print 'g3=', g3
-        
+        attributes = {'depth': [2, 4], 'elevation': [6.1, 5]}
+        attributes1 = {'depth': [2, 4], 'elevation': [2.5, 1]}
+        G1 = Geospatial_data(points, attributes)
+        G2 = Geospatial_data(points, attributes1)
+
         G = G1 + G2
 
@@ -325 +288 @@
 
     def test_addII(self):
-        """ test the addition of two geospatical objects
-            no geo_reference see next test
-        """
+        '''test the addition of two geospatical objects
+        no geo_reference see next test
+        '''
+
         points = [[1.0, 2.1], [3.0, 5.3]]
-        attributes = {'depth':[2, 4]}
-        G1 = Geospatial_data(points, attributes) 
-        
+        attributes = {'depth': [2, 4]}
+        G1 = Geospatial_data(points, attributes)
+
         points = [[5.0, 2.1], [3.0, 50.3]]
-        attributes = {'depth':[200, 400]}
+        attributes = {'depth': [200, 400]}
         G2 = Geospatial_data(points, attributes)
-        
-#        g3 = geospatial_data2points_dictionary(G1)
-#        print 'g3=', g3
-        
+
         G = G1 + G2
 
-        assert G.attributes.has_key('depth') 
+        assert G.attributes.has_key('depth')
         assert G.attributes.keys(), ['depth']
         assert num.allclose(G.attributes['depth'], [2, 4, 200, 400])
         assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3],
                                                   [5.0, 2.1], [3.0, 50.3]])
+
     def test_add_with_geo (self):
-        """
-        Difference in Geo_reference resolved
-        """
+        '''Difference in Geo_reference resolved'''
+
         points1 = [[1.0, 2.1], [3.0, 5.3]]
         points2 = [[5.0, 6.1], [6.0, 3.3]]
@@ -361 +322 @@
                                  projection='UTM',
                                  units='m')
-                                
+
         G1 = Geospatial_data(points1, attributes1, geo_ref1)
         G2 = Geospatial_data(points2, attributes2, geo_ref2)
@@ -370 +331 @@
 
         P2 = G2.get_data_points(absolute=True)
-        assert num.allclose(P2, [[5.1, 9.1], [6.1, 6.3]])        
-        
+        assert num.allclose(P2, [[5.1, 9.1], [6.1, 6.3]])
+
         G = G1 + G2
 
@@ -380 +341 @@
         P = G.get_data_points(absolute=True)
 
-        #P_relative = G.get_data_points(absolute=False)
-        #
-        #assert allclose(P_relative, P - [0.1, 2.0])
-
-        assert num.allclose(P, num.concatenate( (P1,P2) ))
+        assert num.allclose(P, num.concatenate((P1,P2), axis=0))    #??default#
+        msg = ('P=\n%s\nshould be close to\n'
+               '[[2.0, 4.1], [4.0, 7.3],\n'
+               ' [5.1, 9.1], [6.1, 6.3]]'
+               % str(P))
         assert num.allclose(P, [[2.0, 4.1], [4.0, 7.3],
-                                [5.1, 9.1], [6.1, 6.3]])
-        
-
-
-        
+                                [5.1, 9.1], [6.1, 6.3]]), msg
 
     def test_add_with_geo_absolute (self):
-        """
-        Difference in Geo_reference resolved
-        """
+        '''Difference in Geo_reference resolved'''
+
         points1 = num.array([[2.0, 4.1], [4.0, 7.3]])
-        points2 = num.array([[5.1, 9.1], [6.1, 6.3]])        
+        points2 = num.array([[5.1, 9.1], [6.1, 6.3]])
         attributes1 = [2, 4]
         attributes2 = [5, 76]
@@ -403 +359 @@
         geo_ref2 = Geo_reference(55, 2.0, 3.0)
 
-        
-                                
-        G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()],
+        G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(),
+                                        geo_ref1.get_yllcorner()],
                              attributes1, geo_ref1)
-        
-        G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()],
+
+        G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(),
+                                        geo_ref2.get_yllcorner()],
                              attributes2, geo_ref2)
 
@@ -416 +372 @@
 
         P1 = G1.get_data_points(absolute=False)
-        assert num.allclose(P1, points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()])        
+        msg = ('P1=\n%s\nbut points1 - <...>=\n%s'
+               % (str(P1),
+                  str(points1 - [geo_ref1.get_xllcorner(),
+                                 geo_ref1.get_yllcorner()])))
+        assert num.allclose(P1, points1 - [geo_ref1.get_xllcorner(),
+                                           geo_ref1.get_yllcorner()]), msg
 
         P2 = G2.get_data_points(absolute=True)
@@ -422 +383 @@
 
         P2 = G2.get_data_points(absolute=False)
-        assert num.allclose(P2, points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()])                
-        
+        assert num.allclose(P2, points2 - [geo_ref2.get_xllcorner(),
+                                           geo_ref2.get_yllcorner()])
+
         G = G1 + G2
-        
-        #assert allclose(G.get_geo_reference().get_xllcorner(), 1.0)
-        #assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
-
         P = G.get_data_points(absolute=True)
 
-        #P_relative = G.get_data_points(absolute=False)
-        #
-        #assert allclose(P_relative, [[1.0, 2.1], [3.0, 5.3], [4.1, 7.1], [5.1, 4.3]])
-
-        assert num.allclose(P, num.concatenate( (points1,points2) ))
-
+        assert num.allclose(P, num.concatenate((points1, points2), axis=0)) #??default#
 
     def test_add_with_None(self):
-        """ test that None can be added to a geospatical objects
-        """
-        
+        '''test that None can be added to a geospatical objects'''
+
         points1 = num.array([[2.0, 4.1], [4.0, 7.3]])
-        points2 = num.array([[5.1, 9.1], [6.1, 6.3]])        
+        points2 = num.array([[5.1, 9.1], [6.1, 6.3]])
 
         geo_ref1= Geo_reference(55, 1.0, 2.0)
@@ -452 +404 @@
                                  projection='UTM',
                                  units='m')
-        
-
-        attributes1 = {'depth':[2, 4.7], 'elevation':[6.1, 5]}
-        attributes2 = {'depth':[-2.3, 4], 'elevation':[2.5, 1]}
-
+
+        attributes1 = {'depth': [2, 4.7], 'elevation': [6.1, 5]}
+        attributes2 = {'depth': [-2.3, 4], 'elevation': [2.5, 1]}
 
         G1 = Geospatial_data(points1, attributes1, geo_ref1)
@@ -464 +414 @@
         assert G1.attributes.has_key('elevation')
         assert num.allclose(G1.attributes['depth'], [2, 4.7])
-        assert num.allclose(G1.attributes['elevation'], [6.1, 5])        
-        
+        assert num.allclose(G1.attributes['elevation'], [6.1, 5])
+
         G2 = Geospatial_data(points2, attributes2, geo_ref2)
         assert num.allclose(G2.get_geo_reference().get_xllcorner(), 0.1)
@@ -472 +422 @@
         assert G2.attributes.has_key('elevation')
         assert num.allclose(G2.attributes['depth'], [-2.3, 4])
-        assert num.allclose(G2.attributes['elevation'], [2.5, 1])        
+        assert num.allclose(G2.attributes['elevation'], [2.5, 1])
 
         #Check that absolute values are as expected
@@ -479 +429 @@
 
         P2 = G2.get_data_points(absolute=True)
-        assert num.allclose(P2, [[5.2, 12.1], [6.2, 9.3]])        
+        assert num.allclose(P2, [[5.2, 12.1], [6.2, 9.3]])
 
         # Normal add
         G = G1 + None
-
         assert G.attributes.has_key('depth')
         assert G.attributes.has_key('elevation')
         assert num.allclose(G.attributes['depth'], [2, 4.7])
-        assert num.allclose(G.attributes['elevation'], [6.1, 5])        
+        assert num.allclose(G.attributes['elevation'], [6.1, 5])
 
         # Points are now absolute.
         assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
         assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
-        P = G.get_data_points(absolute=True)        
-        assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]])
-
+        P = G.get_data_points(absolute=True)
+        msg = 'P=\n%s' % str(P)
+        assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]]), msg
 
         G = G2 + None
@@ -500 +449 @@
         assert G.attributes.has_key('elevation')
         assert num.allclose(G.attributes['depth'], [-2.3, 4])
-        assert num.allclose(G.attributes['elevation'], [2.5, 1])        
-
+        assert num.allclose(G.attributes['elevation'], [2.5, 1])
         assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
         assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
-        P = G.get_data_points(absolute=True)        
+
+        P = G.get_data_points(absolute=True)
         assert num.allclose(P, [[5.2, 12.1], [6.2, 9.3]])
-        
-
 
         # Reverse add
         G = None + G1
-
         assert G.attributes.has_key('depth')
         assert G.attributes.has_key('elevation')
         assert num.allclose(G.attributes['depth'], [2, 4.7])
-        assert num.allclose(G.attributes['elevation'], [6.1, 5])        
+        assert num.allclose(G.attributes['elevation'], [6.1, 5])
 
         # Points are now absolute.
         assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
         assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
-        P = G.get_data_points(absolute=True)        
-        assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]])        
-
+
+        P = G.get_data_points(absolute=True)
+        assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]])
 
         G = None + G2
@@ -528 +474 @@
         assert G.attributes.has_key('elevation')
         assert num.allclose(G.attributes['depth'], [-2.3, 4])
-        assert num.allclose(G.attributes['elevation'], [2.5, 1])        
+        assert num.allclose(G.attributes['elevation'], [2.5, 1])
 
         assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
         assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
-        P = G.get_data_points(absolute=True)        
+
+        P = G.get_data_points(absolute=True)
         assert num.allclose(P, [[5.2, 12.1], [6.2, 9.3]])
 
-        
-
-        
-                           
-        
     def test_clip0(self):
-        """test_clip0(self):
-        
+        '''test_clip0(self):
+
         Test that point sets can be clipped by a polygon
-        """
-        
+        '''
+
         from anuga.coordinate_transforms.geo_reference import Geo_reference
-        
+
         points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
                   [0, 0], [2.4, 3.3]]
         G = Geospatial_data(points)
 
-        # First try the unit square    
-        U = [[0,0], [1,0], [1,1], [0,1]] 
-        assert num.allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]])
+        # First try the unit square
+        U = [[0,0], [1,0], [1,1], [0,1]]
+        assert num.allclose(G.clip(U).get_data_points(),
+                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
 
         # Then a more complex polygon
         polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
-        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
+                  [0.5, 1.5], [0.5, -0.5]]
         G = Geospatial_data(points)
 
         assert num.allclose(G.clip(polygon).get_data_points(),
-                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
+                            [[0.5, 0.5], [1, -0.5], [1.5, 0]])
 
     def test_clip0_with_attributes(self):
-        """test_clip0_with_attributes(self):
-        
+        '''test_clip0_with_attributes(self):
+
         Test that point sets with attributes can be clipped by a polygon
-        """
-        
+        '''
+
         from anuga.coordinate_transforms.geo_reference import Geo_reference
-        
+
         points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
                   [0, 0], [2.4, 3.3]]
@@ -578 +522 @@
         att_dict = {'att1': attributes,
                     'att2': num.array(attributes)+1}
-        
+
         G = Geospatial_data(points, att_dict)
 
-        # First try the unit square    
-        U = [[0,0], [1,0], [1,1], [0,1]] 
-        assert num.allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]])
+        # First try the unit square
+        U = [[0,0], [1,0], [1,1], [0,1]]
+        assert num.allclose(G.clip(U).get_data_points(),
+                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
         assert num.allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
-        assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])                
+        assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])
 
         # Then a more complex polygon
         polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
-        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
+                  [0.5, 1.5], [0.5, -0.5]]
 
         # This time just one attribute
         attributes = [2, -4, 5, 76, -2, 0.1]
         G = Geospatial_data(points, attributes)
-
         assert num.allclose(G.clip(polygon).get_data_points(),
-                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
+                            [[0.5, 0.5], [1, -0.5], [1.5, 0]])
         assert num.allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
-        
 
     def test_clip1(self):
-        """test_clip1(self):
-        
+        '''test_clip1(self):
+
         Test that point sets can be clipped by a polygon given as
         another Geospatial dataset
-        """
-        
+        '''
+
         from anuga.coordinate_transforms.geo_reference import Geo_reference
-        
+
         points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
                   [0, 0], [2.4, 3.3]]
@@ -615 +559 @@
                     'att2': num.array(attributes)+1}
         G = Geospatial_data(points, att_dict)
-        
-        # First try the unit square    
-        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
+
+        # First try the unit square
+        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
         assert num.allclose(G.clip(U).get_data_points(),
-                        [[0.2, 0.5], [0.4, 0.3], [0, 0]])
-
+                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
         assert num.allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
-        assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])                        
-        
+        assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])
+
         # Then a more complex polygon
-        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
-        attributes = [2, -4, 5, 76, -2, 0.1]        
+        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
+                  [0.5, 1.5], [0.5, -0.5]]
+        attributes = [2, -4, 5, 76, -2, 0.1]
         G = Geospatial_data(points, attributes)
-        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
-        
+        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1],
+                                   [2,1], [0,1]])
 
         assert num.allclose(G.clip(polygon).get_data_points(),
                             [[0.5, 0.5], [1, -0.5], [1.5, 0]])
         assert num.allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
-        
 
     def test_clip0_outside(self):
-        """test_clip0_outside(self):
-        
+        '''test_clip0_outside(self):
+
         Test that point sets can be clipped outside of a polygon
-        """
-        
+        '''
+
         from anuga.coordinate_transforms.geo_reference import Geo_reference
-        
+
         points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
                   [0, 0], [2.4, 3.3]]
-        attributes = [2, -4, 5, 76, -2, 0.1, 3]        
+        attributes = [2, -4, 5, 76, -2, 0.1, 3]
         G = Geospatial_data(points, attributes)
 
-        # First try the unit square    
+        # First try the unit square
         U = [[0,0], [1,0], [1,1], [0,1]]
         assert num.allclose(G.clip_outside(U).get_data_points(),
                             [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
-        #print G.clip_outside(U).get_attributes()
-        assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])        
-        
+        assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])
 
         # Then a more complex polygon
         polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
-        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
-        attributes = [2, -4, 5, 76, -2, 0.1]        
+        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
+                  [0.5, 1.5], [0.5, -0.5]]
+        attributes = [2, -4, 5, 76, -2, 0.1]
         G = Geospatial_data(points, attributes)
-
         assert num.allclose(G.clip_outside(polygon).get_data_points(),
                             [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
-        assert num.allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])                
-
+        assert num.allclose(G.clip_outside(polygon).get_attributes(),
+                            [2, -2, 0.1])
 
     def test_clip1_outside(self):
-        """test_clip1_outside(self):
-        
+        '''test_clip1_outside(self):
+
         Test that point sets can be clipped outside of a polygon given as
         another Geospatial dataset
-        """
-        
+        '''
+
         from anuga.coordinate_transforms.geo_reference import Geo_reference
-        
+
         points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
                   [0, 0], [2.4, 3.3]]
-        attributes = [2, -4, 5, 76, -2, 0.1, 3]        
-        G = Geospatial_data(points, attributes)        
-
-        # First try the unit square    
-        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
+        attributes = [2, -4, 5, 76, -2, 0.1, 3]
+        G = Geospatial_data(points, attributes)
+
+        # First try the unit square
+        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
         assert num.allclose(G.clip_outside(U).get_data_points(),
                             [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
-        assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])        
+        assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])
 
         # Then a more complex polygon
-        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
-        attributes = [2, -4, 5, 76, -2, 0.1]        
+        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
+                  [0.5, 1.5], [0.5, -0.5]]
+        attributes = [2, -4, 5, 76, -2, 0.1]
         G = Geospatial_data(points, attributes)
-
-        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
-        
-
+        polygon = Geospatial_data([[0, 0], [1, 0], [0.5, -1], [2, -1],
+                                   [2, 1], [0, 1]])
         assert num.allclose(G.clip_outside(polygon).get_data_points(),
                             [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
-        assert num.allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])
-        
-
+        assert num.allclose(G.clip_outside(polygon).get_attributes(),
+                            [2, -2, 0.1])
 
     def test_clip1_inside_outside(self):
-        """test_clip1_inside_outside(self):
-        
+        '''test_clip1_inside_outside(self):
+
         Test that point sets can be clipped outside of a polygon given as
         another Geospatial dataset
-        """
-        
+        '''
+
         from anuga.coordinate_transforms.geo_reference import Geo_reference
-        
+
         points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
                   [0, 0], [2.4, 3.3]]
-        attributes = [2, -4, 5, 76, -2, 0.1, 3]        
+        attributes = [2, -4, 5, 76, -2, 0.1, 3]
         G = Geospatial_data(points, attributes)
 
-        # First try the unit square    
-        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
+        # First try the unit square
+        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
         G1 = G.clip(U)
-        assert num.allclose(G1.get_data_points(),[[0.2, 0.5], [0.4, 0.3], [0, 0]])
+        assert num.allclose(G1.get_data_points(),
+                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
         assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])
-        
         G2 = G.clip_outside(U)
         assert num.allclose(G2.get_data_points(),[[-1, 4], [1.0, 2.1],
                                                   [3.0, 5.3], [2.4, 3.3]])
-        assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])                
-
-        
+        assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])
+
         # New ordering
-        new_points = [[0.2, 0.5], [0.4, 0.3], [0, 0]] +\
-                     [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]
-
-        new_attributes = [-4, 76, 0.1, 2, 5, -2, 3]                 
-        
+        new_points = [[0.2, 0.5], [0.4, 0.3], [0, 0], [-1, 4],
+                      [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]
+        new_attributes = [-4, 76, 0.1, 2, 5, -2, 3]
+
         assert num.allclose((G1+G2).get_data_points(), new_points)
         assert num.allclose((G1+G2).get_attributes(), new_attributes)
@@ -741 +678 @@
         G.export_points_file(FN)
 
-
         # Read it back in
         G3 = Geospatial_data(FN)
 
-
         # Check result
-        assert num.allclose(G3.get_data_points(), new_points)        
-        assert num.allclose(G3.get_attributes(), new_attributes)        
-        
+        assert num.allclose(G3.get_data_points(), new_points)
+        assert num.allclose(G3.get_attributes(), new_attributes)
+
         os.remove(FN)
 
-        
     def test_load_csv(self):
-        
-        import os
-        import tempfile
-       
-        fileName = tempfile.mktemp(".csv")
-        file = open(fileName,"w")
-        file.write("x,y,elevation speed \n\
+        fileName = tempfile.mktemp('.csv')
+        file = open(fileName,'w')
+        file.write('x,y,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")
+1.0 0.0 10.4 40.0\n')
         file.close()
-        #print fileName
+
         results = Geospatial_data(fileName)
         os.remove(fileName)
-#        print 'data', results.get_data_points()
-        assert num.allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],
-                                                        [1.0, 0.0]])
+        assert num.allclose(results.get_data_points(),
+                            [[1.0, 0.0],[0.0, 1.0], [1.0, 0.0]])
         assert num.allclose(results.get_attributes(attribute_name='elevation'),
                             [10.0, 0.0, 10.4])
@@ -776 +705 @@
                             [0.0, 10.0, 40.0])
 
-
-  ###################### .CSV ##############################
+################################################################################
+# Test CSV files
+################################################################################
 
     def test_load_csv_lat_long_bad_blocking(self):
-        """
-        test_load_csv_lat_long_bad_blocking(self):
+        '''test_load_csv_lat_long_bad_blocking(self):
         Different zones in Geo references
-        """
-        fileName = tempfile.mktemp(".csv")
-        file = open(fileName,"w")
-        file.write("Lati,LONG,z \n\
+        '''
+
+        fileName = tempfile.mktemp('.csv')
+        file = open(fileName, 'w')
+        file.write('Lati,LONG,z \n\
 -25.0,180.0,452.688000\n\
--34,150.0,459.126000\n")
+-34,150.0,459.126000\n')
         file.close()
-        
+
         results = Geospatial_data(fileName, max_read_lines=1,
                                   load_file_now=False)
-        
-        #for i in results:
-        #    pass
+
         try:
             for i in results:
@@ -803 +731 @@
         else:
             msg = 'Different zones in Geo references not caught.'
-            raise msg        
-        
-        os.remove(fileName)
-        
+            raise Exception, msg
+
+        os.remove(fileName)
+
     def test_load_csv(self):
-        
-        fileName = tempfile.mktemp(".txt")
-        file = open(fileName,"w")
-        file.write(" x,y, 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")
+        fileName = tempfile.mktemp('.txt')
+        file = open(fileName, 'w')
+        file.write(' x,y, 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()
 
         results = Geospatial_data(fileName, max_read_lines=2)
 
-
-        assert num.allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
-        assert num.allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
-        assert num.allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
+        assert num.allclose(results.get_data_points(),
+                            [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        assert num.allclose(results.get_attributes(attribute_name='elevation'),
+                            [10.0, 0.0, 10.4])
+        assert num.allclose(results.get_attributes(attribute_name='speed'),
+                            [0.0, 10.0, 40.0])
 
         # Blocking
@@ -828 +757 @@
         for i in results:
             geo_list.append(i)
-            
+
         assert num.allclose(geo_list[0].get_data_points(),
                             [[1.0, 0.0],[0.0, 1.0]])
-
         assert num.allclose(geo_list[0].get_attributes(attribute_name='elevation'),
                             [10.0, 0.0])
         assert num.allclose(geo_list[1].get_data_points(),
-                            [[1.0, 0.0]])        
+                            [[1.0, 0.0]])
         assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'),
                             [10.4])
-           
-        os.remove(fileName)         
-        
+
+        os.remove(fileName)
+
     def test_load_csv_bad(self):
-        """ test_load_csv_bad(self):
+        '''test_load_csv_bad(self):
         header column, body column missmatch
         (Format error)
-        """
-        import os
-       
-        fileName = tempfile.mktemp(".txt")
-        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")
+        '''
+
+        fileName = tempfile.mktemp('.txt')
+        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()
 
@@ -861 +788 @@
         # Blocking
         geo_list = []
-        #for i in results:
-        #    geo_list.append(i)
         try:
             for i in results:
@@ -869 +794 @@
             pass
         else:
-            msg = 'bad file did not raise error!'
-            raise msg
+            msg = 'Bad file did not raise error!'
+            raise Exception, msg
         os.remove(fileName)
 
     def test_load_csv_badII(self):
-        """ test_load_csv_bad(self):
+        '''test_load_csv_bad(self):
         header column, body column missmatch
         (Format error)
-        """
-        import os
-       
-        fileName = tempfile.mktemp(".txt")
-        file = open(fileName,"w")
-        file.write(" x,y,elevation ,  speed \n\
+        '''
+
+        fileName = tempfile.mktemp('.txt')
+        file = open(fileName, 'w')
+        file.write(' x,y,elevation ,  speed \n\
 1.0, 0.0, 10.0, 0.0\n\
 0.0, 1.0, 0.0, 10\n\
-1.0, 0.0 ,10.4 yeah\n")
+1.0, 0.0 ,10.4 yeah\n')
         file.close()
 
@@ -893 +817 @@
         # Blocking
         geo_list = []
-        #for i in results:
-        #    geo_list.append(i)
         try:
             for i in results:
@@ -901 +823 @@
             pass
         else:
-            msg = 'bad file did not raise error!'
-            raise msg
+            msg = 'Bad file did not raise error!'
+            raise Exception, msg
         os.remove(fileName)
 
     def test_load_csv_badIII(self):
-        """ test_load_csv_bad(self):
+        '''test_load_csv_bad(self):
         space delimited
-        """
-        import os
-       
-        fileName = tempfile.mktemp(".txt")
-        file = open(fileName,"w")
-        file.write(" x y elevation   speed \n\
+        '''
+
+        fileName = tempfile.mktemp('.txt')
+        file = open(fileName, 'w')
+        file.write(' x y elevation   speed \n\
 1. 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")
+1.0 0.0 10.4 40.0\n')
         file.close()
 
@@ -925 +846 @@
             pass
         else:
-            msg = 'bad file did not raise error!'
-            raise msg
-        os.remove(fileName)
-        
+            msg = 'Bad file did not raise error!'
+            raise Exception, msg
+        os.remove(fileName)
+
     def test_load_csv_badIV(self):
-        """ test_load_csv_bad(self):
+        ''' test_load_csv_bad(self):
         header column, body column missmatch
         (Format error)
-        """
-        import os
-       
-        fileName = tempfile.mktemp(".txt")
-        file = open(fileName,"w")
-        file.write(" x,y,elevation ,  speed \n\
+        '''
+
+        fileName = tempfile.mktemp('.txt')
+        file = open(fileName, 'w')
+        file.write(' x,y,elevation ,  speed \n\
 1.0, 0.0, 10.0, wow\n\
 0.0, 1.0, 0.0, ha\n\
-1.0, 0.0 ,10.4, yeah\n")
+1.0, 0.0 ,10.4, yeah\n')
         file.close()
 
@@ -949 +869 @@
         # Blocking
         geo_list = []
-        #for i in results:
-         #   geo_list.append(i)
         try:
             for i in results:
@@ -957 +875 @@
             pass
         else:
-            msg = 'bad file did not raise error!'
-            raise msg
+            msg = 'Bad file did not raise error!'
+            raise Exception, msg
         os.remove(fileName)
 
     def test_load_pts_blocking(self):
         #This is pts!
-       
-        import os
-       
-        fileName = tempfile.mktemp(".txt")
-        file = open(fileName,"w")
-        file.write(" x,y, 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")
+        fileName = tempfile.mktemp('.txt')
+        file = open(fileName, 'w')
+        file.write(' x,y, 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()
 
-        pts_file = tempfile.mktemp(".pts")
-        
+        pts_file = tempfile.mktemp('.pts')
+
         convert = Geospatial_data(fileName)
         convert.export_points_file(pts_file)
         results = Geospatial_data(pts_file, max_read_lines=2)
 
-        assert num.allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],
-                                                        [1.0, 0.0]])
+        assert num.allclose(results.get_data_points(),
+                            [[1.0, 0.0],[0.0, 1.0], [1.0, 0.0]])
         assert num.allclose(results.get_attributes(attribute_name='elevation'),
                             [10.0, 0.0, 10.4])
@@ -990 +905 @@
         geo_list = []
         for i in results:
-            geo_list.append(i) 
+            geo_list.append(i)
         assert num.allclose(geo_list[0].get_data_points(),
                             [[1.0, 0.0],[0.0, 1.0]])
@@ -996 +911 @@
                             [10.0, 0.0])
         assert num.allclose(geo_list[1].get_data_points(),
-                            [[1.0, 0.0]])        
+                            [[1.0, 0.0]])
         assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'),
                             [10.4])
-           
-        os.remove(fileName)  
-        os.remove(pts_file)               
+
+        os.remove(fileName)
+        os.remove(pts_file)
 
     def verbose_test_load_pts_blocking(self):
-        
-        import os
-       
-        fileName = tempfile.mktemp(".txt")
-        file = open(fileName,"w")
-        file.write(" x,y, 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.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\n\
-1.0, 0.0, 10.0, 0.0\n\
-0.0, 1.0, 0.0, 10.0\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")
+        fileName = tempfile.mktemp('.txt')
+        file = open(fileName, 'w')
+        file.write(' x,y, 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.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\n\
+1.0, 0.0, 10.0, 0.0\n\
+0.0, 1.0, 0.0, 10.0\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()
 
-        pts_file = tempfile.mktemp(".pts")
-        
+        pts_file = tempfile.mktemp('.pts')
         convert = Geospatial_data(fileName)
         convert.export_points_file(pts_file)
@@ -1068 +979 @@
         geo_list = []
         for i in results:
-            geo_list.append(i) 
+            geo_list.append(i)
         assert num.allclose(geo_list[0].get_data_points(),
-                            [[1.0, 0.0],[0.0, 1.0]])
+                            [[1.0, 0.0], [0.0, 1.0]])
         assert num.allclose(geo_list[0].get_attributes(attribute_name='elevation'),
                             [10.0, 0.0])
         assert num.allclose(geo_list[1].get_data_points(),
-                            [[1.0, 0.0],[0.0, 1.0] ])        
+                            [[1.0, 0.0],[0.0, 1.0] ])
         assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'),
                             [10.0, 0.0])
-           
-        os.remove(fileName)  
+
+        os.remove(fileName)
         os.remove(pts_file)
-        
-        
 
     def test_new_export_pts_file(self):
@@ -1088 +997 @@
         att_dict['elevation'] = num.array([10.1, 0.0, 10.4])
         att_dict['brightness'] = num.array([10.0, 1.0, 10.4])
-        
-        fileName = tempfile.mktemp(".pts")
-        
+
+        fileName = tempfile.mktemp('.pts')
         G = Geospatial_data(pointlist, att_dict)
-        
         G.export_points_file(fileName)
-
         results = Geospatial_data(file_name = fileName)
-
-        os.remove(fileName)
-        
-        assert num.allclose(results.get_data_points(),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
-        assert num.allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
+        os.remove(fileName)
+
+        assert num.allclose(results.get_data_points(),
+                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
+        assert num.allclose(results.get_attributes(attribute_name='elevation'),
+                            [10.1, 0.0, 10.4])
         answer = [10.0, 1.0, 10.4]
-        assert num.allclose(results.get_attributes(attribute_name='brightness'), answer)
+        assert num.allclose(results.get_attributes(attribute_name='brightness'),
+                            answer)
 
     def test_new_export_absolute_pts_file(self):
         att_dict = {}
-        pointlist = num.array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
         att_dict['elevation'] = num.array([10.1, 0.0, 10.4])
         att_dict['brightness'] = num.array([10.0, 1.0, 10.4])
         geo_ref = Geo_reference(50, 25, 55)
-        
-        fileName = tempfile.mktemp(".pts")
-        
+
+        fileName = tempfile.mktemp('.pts')
         G = Geospatial_data(pointlist, att_dict, geo_ref)
-        
         G.export_points_file(fileName, absolute=True)
-
         results = Geospatial_data(file_name = fileName)
-
-        os.remove(fileName)
-        
-        assert num.allclose(results.get_data_points(), G.get_data_points(True))
-        assert num.allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4])
+        os.remove(fileName)
+
+        msg = ('results.get_data_points()=\n%s\nbut G.get_data_points(True)=\n%s'
+               % (str(results.get_data_points()), str(G.get_data_points(True))))
+        assert num.allclose(results.get_data_points(),
+                            G.get_data_points(True)), msg
+        msg = ("results.get_attributes(attribute_name='elevation')=%s"
+               % str(results.get_attributes(attribute_name='elevation')))
+        assert num.allclose(results.get_attributes(attribute_name='elevation'),
+                            [10.1, 0.0, 10.4]), msg
         answer = [10.0, 1.0, 10.4]
-        assert num.allclose(results.get_attributes(attribute_name='brightness'), answer)
+        msg = ("results.get_attributes(attribute_name='brightness')=%s, "
+               'answer=%s' %
+               (str(results.get_attributes(attribute_name='brightness')),
+                str(answer)))
+        assert num.allclose(results.get_attributes(attribute_name='brightness'),
+                            answer), msg
 
     def test_loadpts(self):
-        
         from Scientific.IO.NetCDF import NetCDFFile
 
-        fileName = tempfile.mktemp(".pts")
+        fileName = tempfile.mktemp('.pts')
         # NetCDF file definition
         outfile = NetCDFFile(fileName, netcdf_mode_w)
-        
+
         # dimension definitions
-        outfile.createDimension('number_of_points', 3)    
-        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
-    
+        outfile.createDimension('number_of_points', 3)
+        outfile.createDimension('number_of_dimensions', 2) # This is 2d data
+
         # variable definitions
-        outfile.createVariable('points', num.Float, ('number_of_points',
-                                                     'number_of_dimensions'))
-        outfile.createVariable('elevation', num.Float, ('number_of_points',))
-    
+        outfile.createVariable('points', netcdf_float, ('number_of_points',
+                                                        'number_of_dimensions'))
+        outfile.createVariable('elevation', netcdf_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 
+        elevation[0] = 10.0
         points[1, :] = [0.0,1.0]
-        elevation[1] = 0.0  
+        elevation[1] = 0.0
         points[2, :] = [1.0,0.0]
-        elevation[2] = 10.4    
+        elevation[2] = 10.4
 
         outfile.close()
-        
+
         results = Geospatial_data(file_name = fileName)
         os.remove(fileName)
-        answer =  [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
-        assert num.allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
-        assert num.allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
-        
+
+        answer = [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]
+        assert num.allclose(results.get_data_points(),
+                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
+        assert num.allclose(results.get_attributes(attribute_name='elevation'),
+                            [10.0, 0.0, 10.4])
+
     def test_writepts(self):
-        #test_writepts: Test that storage of x,y,attributes works
-        
+        '''Test that storage of x,y,attributes works'''
+
         att_dict = {}
-        pointlist = num.array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
         att_dict['elevation'] = num.array([10.0, 0.0, 10.4])
         att_dict['brightness'] = num.array([10.0, 0.0, 10.4])
-        geo_reference=Geo_reference(56,1.9,1.9)
+        geo_reference=Geo_reference(56, 1.9, 1.9)
 
         # Test pts format
-        fileName = tempfile.mktemp(".pts")
+        fileName = tempfile.mktemp('.pts')
         G = Geospatial_data(pointlist, att_dict, geo_reference)
         G.export_points_file(fileName, False)
@@ -1178 +1095 @@
         os.remove(fileName)
 
-        assert num.allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
-        assert num.allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
+        assert num.allclose(results.get_data_points(False),
+                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
+        assert num.allclose(results.get_attributes('elevation'),
+                            [10.0, 0.0, 10.4])
         answer = [10.0, 0.0, 10.4]
         assert num.allclose(results.get_attributes('brightness'), answer)
         self.failUnless(geo_reference == geo_reference,
-                         'test_writepts failed. Test geo_reference')
+                        'test_writepts failed. Test geo_reference')
 
     def test_write_csv_attributes(self):
-        #test_write : Test that storage of x,y,attributes works
-        
+        '''Test that storage of x,y,attributes works'''
+
         att_dict = {}
-        pointlist = num.array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
         att_dict['elevation'] = num.array([10.0, 0.0, 10.4])
         att_dict['brightness'] = num.array([10.0, 0.0, 10.4])
-        geo_reference=Geo_reference(56,0,0)
+        geo_reference=Geo_reference(56, 0, 0)
+
         # Test txt format
-        fileName = tempfile.mktemp(".txt")
+        fileName = tempfile.mktemp('.txt')
         G = Geospatial_data(pointlist, att_dict, geo_reference)
         G.export_points_file(fileName)
-        #print "fileName", fileName 
         results = Geospatial_data(file_name=fileName)
         os.remove(fileName)
-        assert num.allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
-        assert num.allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
+
+        assert num.allclose(results.get_data_points(False),
+                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
+        assert num.allclose(results.get_attributes('elevation'),
+                            [10.0, 0.0, 10.4])
         answer = [10.0, 0.0, 10.4]
         assert num.allclose(results.get_attributes('brightness'), answer)
-        
- 
+
     def test_write_csv_attributes_lat_long(self):
-        #test_write : Test that storage of x,y,attributes works
-        
+        '''Test that storage of x,y,attributes works'''
+
         att_dict = {}
-        pointlist = num.array([[-21.5,114.5],[-21.6,114.5],[-21.7,114.5]])
+        pointlist = num.array([[-21.5, 114.5], [-21.6, 114.5], [-21.7, 114.5]])
         att_dict['elevation'] = num.array([10.0, 0.0, 10.4])
         att_dict['brightness'] = num.array([10.0, 0.0, 10.4])
+
         # Test txt format
-        fileName = tempfile.mktemp(".txt")
+        fileName = tempfile.mktemp('.txt')
         G = Geospatial_data(pointlist, att_dict, points_are_lats_longs=True)
         G.export_points_file(fileName, as_lat_long=True)
-        #print "fileName", fileName 
         results = Geospatial_data(file_name=fileName)
         os.remove(fileName)
+
         assert num.allclose(results.get_data_points(False, as_lat_long=True),
-                        pointlist)
-        assert num.allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4])
+                            pointlist)
+        assert num.allclose(results.get_attributes('elevation'),
+                            [10.0, 0.0, 10.4])
         answer = [10.0, 0.0, 10.4]
         assert num.allclose(results.get_attributes('brightness'), answer)
-        
+
     def test_writepts_no_attributes(self):
-
-        #test_writepts_no_attributes: Test that storage of x,y alone works
-        
+        '''Test that storage of x,y alone works'''
+
         att_dict = {}
-        pointlist = num.array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
-        geo_reference=Geo_reference(56,1.9,1.9)
+        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
+        geo_reference=Geo_reference(56, 1.9, 1.9)
 
         # Test pts format
-        fileName = tempfile.mktemp(".pts")
+        fileName = tempfile.mktemp('.pts')
         G = Geospatial_data(pointlist, None, geo_reference)
         G.export_points_file(fileName, False)
@@ -1241 +1163 @@
         os.remove(fileName)
 
-        assert num.allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        assert num.allclose(results.get_data_points(False),
+                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
         self.failUnless(geo_reference == geo_reference,
-                         'test_writepts failed. Test geo_reference')
-        
-       
+                        'test_writepts failed. Test geo_reference')
+
     def test_write_csv_no_attributes(self):
-        #test_write txt _no_attributes: Test that storage of x,y alone works
-        
+        '''Test that storage of x,y alone works'''
+
         att_dict = {}
-        pointlist = num.array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
         geo_reference=Geo_reference(56,0,0)
+
         # Test format
-        fileName = tempfile.mktemp(".txt")
+        fileName = tempfile.mktemp('.txt')
         G = Geospatial_data(pointlist, None, geo_reference)
         G.export_points_file(fileName)
         results = Geospatial_data(file_name=fileName)
         os.remove(fileName)
-        assert num.allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
-
-         
-        
- ########################## BAD .PTS ##########################          
+
+        assert num.allclose(results.get_data_points(False),
+                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
+
+################################################################################
+# Check bad PTS files.
+################################################################################
 
     def test_load_bad_no_file_pts(self):
-        import os
-        import tempfile
-       
-        fileName = tempfile.mktemp(".pts")
-        #print fileName
+        fileName = tempfile.mktemp('.pts')
         try:
-            results = Geospatial_data(file_name = fileName)
-#            dict = import_points_file(fileName)
+            results = Geospatial_data(file_name=fileName)
         except IOError:
             pass
         else:
             msg = 'imaginary file did not raise error!'
-            raise msg
-#            self.failUnless(0 == 1,
-#                        'imaginary file did not raise error!')
-
+            raise Exception, msg
 
     def test_create_from_pts_file(self):
-        
         from Scientific.IO.NetCDF import NetCDFFile
 
-#        fileName = tempfile.mktemp(".pts")
+        # NetCDF file definition
         FN = 'test_points.pts'
-        # NetCDF file definition
         outfile = NetCDFFile(FN, netcdf_mode_w)
-        
+
         # dimension definitions
-        outfile.createDimension('number_of_points', 3)    
-        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
-    
+        outfile.createDimension('number_of_points', 3)
+        outfile.createDimension('number_of_dimensions', 2) # This is 2d data
+
         # variable definitions
-        outfile.createVariable('points', num.Float, ('number_of_points',
-                                                     'number_of_dimensions'))
-        outfile.createVariable('elevation', num.Float, ('number_of_points',))
-    
+        outfile.createVariable('points', netcdf_float, ('number_of_points',
+                                                        'number_of_dimensions'))
+        outfile.createVariable('elevation', netcdf_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 
+        elevation[0] = 10.0
         points[1, :] = [0.0,1.0]
-        elevation[1] = 0.0  
+        elevation[1] = 0.0
         points[2, :] = [1.0,0.0]
-        elevation[2] = 10.4    
+        elevation[2] = 10.4
 
         outfile.close()
@@ -1317 +1232 @@
         assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
         assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
-
-        assert num.allclose(G.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        assert num.allclose(G.get_data_points(),
+                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
         assert num.allclose(G.get_attributes(), [10.0, 0.0, 10.4])
         os.remove(FN)
 
     def test_create_from_pts_file_with_geo(self):
-        """This test reveals if Geospatial data is correctly instantiated from a pts file.
-        """
-        
+        '''Test if Geospatial data is correctly instantiated from a pts file.'''
+
         from Scientific.IO.NetCDF import NetCDFFile
 
+        # NetCDF file definition
         FN = 'test_points.pts'
-        # NetCDF file definition
         outfile = NetCDFFile(FN, netcdf_mode_w)
 
@@ -1339 +1253 @@
 
         # dimension definitions
-        outfile.createDimension('number_of_points', 3)    
-        outfile.createDimension('number_of_dimensions', 2) #This is 2d data
-    
+        outfile.createDimension('number_of_points', 3)
+        outfile.createDimension('number_of_dimensions', 2) # This is 2d data
+
         # variable definitions
-        outfile.createVariable('points', num.Float, ('number_of_points',
-                                                     'number_of_dimensions'))
-        outfile.createVariable('elevation', num.Float, ('number_of_points',))
-    
+        outfile.createVariable('points', netcdf_float, ('number_of_points',
+                                                        'number_of_dimensions'))
+        outfile.createVariable('elevation', netcdf_float, ('number_of_points',))
+
         # Get handles to the variables
         points = outfile.variables['points']
@@ -1352 +1266 @@
 
         points[0, :] = [1.0,0.0]
-        elevation[0] = 10.0 
+        elevation[0] = 10.0
         points[1, :] = [0.0,1.0]
-        elevation[1] = 0.0  
+        elevation[1] = 0.0
         points[2, :] = [1.0,0.0]
-        elevation[2] = 10.4    
+        elevation[2] = 10.4
 
         outfile.close()
@@ -1364 +1278 @@
         assert num.allclose(G.get_geo_reference().get_xllcorner(), xll)
         assert num.allclose(G.get_geo_reference().get_yllcorner(), yll)
-
         assert num.allclose(G.get_data_points(), [[1.0+xll, 0.0+yll],
                                                   [0.0+xll, 1.0+yll],
                                                   [1.0+xll, 0.0+yll]])
-        
         assert num.allclose(G.get_attributes(), [10.0, 0.0, 10.4])
+
         os.remove(FN)
 
-        
     def test_add_(self):
         '''test_add_(self):
         adds an txt and pts files, reads the files and adds them
-           checking results are correct
+        checking results are correct
         '''
+
         # create files
         att_dict1 = {}
-        pointlist1 = num.array([[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+        pointlist1 = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
         att_dict1['elevation'] = num.array([-10.0, 0.0, 10.4])
         att_dict1['brightness'] = num.array([10.0, 0.0, 10.4])
         geo_reference1 = Geo_reference(56, 2.0, 1.0)
-        
+
         att_dict2 = {}
-        pointlist2 = num.array([[2.0, 1.0],[1.0, 2.0],[2.0, 1.0]])
+        pointlist2 = num.array([[2.0, 1.0], [1.0, 2.0], [2.0, 1.0]])
         att_dict2['elevation'] = num.array([1.0, 15.0, 1.4])
         att_dict2['brightness'] = num.array([14.0, 1.0, -12.4])
-        geo_reference2 = Geo_reference(56, 1.0, 2.0) 
+        geo_reference2 = Geo_reference(56, 1.0, 2.0)
 
         G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1)
         G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2)
-        
-        fileName1 = tempfile.mktemp(".txt")
-        fileName2 = tempfile.mktemp(".pts")
-
-        #makes files
+
+        fileName1 = tempfile.mktemp('.txt')
+        fileName2 = tempfile.mktemp('.pts')
+
+        # makes files
         G1.export_points_file(fileName1)
         G2.export_points_file(fileName2)
-        
+
         # add files
-        
-        G3 = Geospatial_data(file_name = fileName1)
-        G4 = Geospatial_data(file_name = fileName2)
-        
+        G3 = Geospatial_data(file_name=fileName1)
+        G4 = Geospatial_data(file_name=fileName2)
         G = G3 + G4
 
-        
         #read results
-#        print'res', G.get_data_points()
-#        print'res1', G.get_data_points(False)
         assert num.allclose(G.get_data_points(),
                             [[ 3.0, 1.0], [ 2.0, 2.0],
                              [ 3.0, 1.0], [ 3.0, 3.0],
                              [ 2.0, 4.0], [ 3.0, 3.0]])
-                         
         assert num.allclose(G.get_attributes(attribute_name='elevation'),
                             [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4])
-        
         answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4]
-        assert num.allclose(G.get_attributes(attribute_name='brightness'), answer)
-        
+        assert num.allclose(G.get_attributes(attribute_name='brightness'),
+                            answer)
         self.failUnless(G.get_geo_reference() == geo_reference1,
-                         'test_writepts failed. Test geo_reference')
-                         
+                        'test_writepts failed. Test geo_reference')
+
         os.remove(fileName1)
         os.remove(fileName2)
-        
+
     def test_ensure_absolute(self):
-        points = [[2.0, 0.0],[1.0, 1.0],
-                         [2.0, 0.0],[2.0, 2.0],
-                         [1.0, 3.0],[2.0, 2.0]]
+        points = [[2.0, 0.0], [1.0, 1.0],
+                  [2.0, 0.0], [2.0, 2.0],
+                  [1.0, 3.0], [2.0, 2.0]]
         new_points = ensure_absolute(points)
-        
+
         assert num.allclose(new_points, points)
 
-        points = num.array([[2.0, 0.0],[1.0, 1.0],
-                            [2.0, 0.0],[2.0, 2.0],
-                            [1.0, 3.0],[2.0, 2.0]])
+        points = num.array([[2.0, 0.0], [1.0, 1.0],
+                            [2.0, 0.0], [2.0, 2.0],
+                            [1.0, 3.0], [2.0, 2.0]])
         new_points = ensure_absolute(points)
-        
+
         assert num.allclose(new_points, points)
-        
-        ab_points = num.array([[2.0, 0.0],[1.0, 1.0],
-                               [2.0, 0.0],[2.0, 2.0],
-                               [1.0, 3.0],[2.0, 2.0]])
-        
-        mesh_origin = (56, 290000, 618000) #zone, easting, northing
-
-        data_points = num.zeros((ab_points.shape), num.Float)
+
+        ab_points = num.array([[2.0, 0.0], [1.0, 1.0],
+                               [2.0, 0.0], [2.0, 2.0],
+                               [1.0, 3.0], [2.0, 2.0]])
+
+        mesh_origin = (56, 290000, 618000) # zone, easting, northing
+        data_points = num.zeros((ab_points.shape), num.float)
+
         #Shift datapoints according to new origins
         for k in range(len(ab_points)):
             data_points[k][0] = ab_points[k][0] - mesh_origin[1]
             data_points[k][1] = ab_points[k][1] - mesh_origin[2]
-        #print "data_points",data_points     
-        new_points = ensure_absolute(data_points,
-                                     geo_reference=mesh_origin)
-        #print "new_points",new_points
-        #print "ab_points",ab_points
-           
+        new_points = ensure_absolute(data_points, geo_reference=mesh_origin)
+
         assert num.allclose(new_points, ab_points)
 
         geo = Geo_reference(56,67,-56)
-
-        data_points = geo.change_points_geo_ref(ab_points)   
-        new_points = ensure_absolute(data_points,
-                                     geo_reference=geo)
-        #print "new_points",new_points
-        #print "ab_points",ab_points
-           
+        data_points = geo.change_points_geo_ref(ab_points)
+        new_points = ensure_absolute(data_points, geo_reference=geo)
+
         assert num.allclose(new_points, ab_points)
-
 
         geo_reference = Geo_reference(56, 100, 200)
@@ -1478 +1375 @@
         points = geo_reference.change_points_geo_ref(ab_points)
         attributes = [2, 4]
-        #print "geo in points", points
-        G = Geospatial_data(points, attributes,
-                            geo_reference=geo_reference)
-          
+        G = Geospatial_data(points, attributes, geo_reference=geo_reference)
         new_points = ensure_absolute(G)
-        #print "new_points",new_points
-        #print "ab_points",ab_points
-           
+
         assert num.allclose(new_points, ab_points)
 
-
-        
     def test_ensure_geospatial(self):
-        points = [[2.0, 0.0],[1.0, 1.0],
-                         [2.0, 0.0],[2.0, 2.0],
-                         [1.0, 3.0],[2.0, 2.0]]
+        points = [[2.0, 0.0], [1.0, 1.0],
+                  [2.0, 0.0], [2.0, 2.0],
+                  [1.0, 3.0], [2.0, 2.0]]
         new_points = ensure_geospatial(points)
-        
-        assert num.allclose(new_points.get_data_points(absolute = True), points)
-
-        points = num.array([[2.0, 0.0],[1.0, 1.0],
-                            [2.0, 0.0],[2.0, 2.0],
-                            [1.0, 3.0],[2.0, 2.0]])
+        assert num.allclose(new_points.get_data_points(absolute=True), points)
+
+        points = num.array([[2.0, 0.0], [1.0, 1.0],
+                            [2.0, 0.0], [2.0, 2.0],
+                            [1.0, 3.0], [2.0, 2.0]])
         new_points = ensure_geospatial(points)
-        
-        assert num.allclose(new_points.get_data_points(absolute = True), points)
-        
+        assert num.allclose(new_points.get_data_points(absolute=True), points)
+
         ab_points = num.array([[2.0, 0.0],[1.0, 1.0],
                                [2.0, 0.0],[2.0, 2.0],
                                [1.0, 3.0],[2.0, 2.0]])
-        
-        mesh_origin = (56, 290000, 618000) #zone, easting, northing
-
-        data_points = num.zeros((ab_points.shape), num.Float)
+        mesh_origin = (56, 290000, 618000)      # zone, easting, northing
+        data_points = num.zeros((ab_points.shape), num.float)
+
         #Shift datapoints according to new origins
         for k in range(len(ab_points)):
             data_points[k][0] = ab_points[k][0] - mesh_origin[1]
             data_points[k][1] = ab_points[k][1] - mesh_origin[2]
-        #print "data_points",data_points     
         new_geospatial = ensure_geospatial(data_points,
-                                             geo_reference=mesh_origin)
+                                           geo_reference=mesh_origin)
         new_points = new_geospatial.get_data_points(absolute=True)
-        #print "new_points",new_points
-        #print "ab_points",ab_points
-           
         assert num.allclose(new_points, ab_points)
 
-        geo = Geo_reference(56,67,-56)
-
-        data_points = geo.change_points_geo_ref(ab_points)   
-        new_geospatial = ensure_geospatial(data_points,
-                                             geo_reference=geo)
+        geo = Geo_reference(56, 67, -56)
+        data_points = geo.change_points_geo_ref(ab_points)
+        new_geospatial = ensure_geospatial(data_points, geo_reference=geo)
         new_points = new_geospatial.get_data_points(absolute=True)
-        #print "new_points",new_points
-        #print "ab_points",ab_points
-           
-        assert num.allclose(new_points, ab_points)
-
+        msg = ('new_points=\n%s\nab_points=\n%s'
+               % (str(new_points), str(ab_points)))
+        assert num.allclose(new_points, ab_points), msg
 
         geo_reference = Geo_reference(56, 100, 200)
@@ -1541 +1420 @@
         points = geo_reference.change_points_geo_ref(ab_points)
         attributes = [2, 4]
-        #print "geo in points", points
-        G = Geospatial_data(points, attributes,
-                            geo_reference=geo_reference)
-          
-        new_geospatial  = ensure_geospatial(G)
+        G = Geospatial_data(points, attributes, geo_reference=geo_reference)
+        new_geospatial = ensure_geospatial(G)
         new_points = new_geospatial.get_data_points(absolute=True)
-        #print "new_points",new_points
-        #print "ab_points",ab_points
-           
         assert num.allclose(new_points, ab_points)
-        
+
     def test_isinstance(self):
-
-        import os
-       
-        fileName = tempfile.mktemp(".csv")
-        file = open(fileName,"w")
-        file.write("x,y,  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")
+        fileName = tempfile.mktemp('.csv')
+        file = open(fileName, 'w')
+        file.write('x,y,  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()
 
         results = Geospatial_data(fileName)
         assert num.allclose(results.get_data_points(absolute=True),
-                            [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
+                            [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]])
         assert num.allclose(results.get_attributes(attribute_name='elevation'),
                             [10.0, 0.0, 10.4])
@@ -1573 +1443 @@
 
         os.remove(fileName)
-        
 
     def test_no_constructors(self):
-        
         try:
             G = Geospatial_data()
-#            results = Geospatial_data(file_name = fileName)
-#            dict = import_points_file(fileName)
         except ValueError:
             pass
         else:
             msg = 'Instance must have a filename or data points'
-            raise msg        
+            raise Exception, msg
 
     def test_load_csv_lat_long(self):
-        """ 
-        comma delimited
-
-        """
-        fileName = tempfile.mktemp(".csv")
-        file = open(fileName,"w")
-        file.write("long,lat, elevation, yeah \n\
+        '''comma delimited'''
+
+        fileName = tempfile.mktemp('.csv')
+        file = open(fileName, 'w')
+        file.write('long,lat, elevation, yeah \n\
 150.916666667,-34.50,452.688000, 10\n\
-150.0,-34,459.126000, 10\n")
+150.0,-34,459.126000, 10\n')
         file.close()
+
         results = Geospatial_data(fileName)
         os.remove(fileName)
         points = results.get_data_points()
-        
+
         assert num.allclose(points[0][0], 308728.009)
         assert num.allclose(points[0][1], 6180432.601)
         assert num.allclose(points[1][0],  222908.705)
         assert num.allclose(points[1][1], 6233785.284)
-        
-      
+
     def test_load_csv_lat_longII(self):
-        """ 
-        comma delimited
-
-        """
-        fileName = tempfile.mktemp(".csv")
-        file = open(fileName,"w")
-        file.write("Lati,LONG,z \n\
+        '''comma delimited'''
+
+        fileName = tempfile.mktemp('.csv')
+        file = open(fileName, 'w')
+        file.write('Lati,LONG,z \n\
 -34.50,150.916666667,452.688000\n\
--34,150.0,459.126000\n")
+-34,150.0,459.126000\n')
         file.close()
+
         results = Geospatial_data(fileName)
         os.remove(fileName)
         points = results.get_data_points()
-        
+
         assert num.allclose(points[0][0], 308728.009)
         assert num.allclose(points[0][1], 6180432.601)
-        assert num.allclose(points[1][0],  222908.705)
+        assert num.allclose(points[1][0], 222908.705)
         assert num.allclose(points[1][1], 6233785.284)
 
-          
     def test_load_csv_lat_long_bad(self):
-        """ 
-        comma delimited
-
-        """
-        fileName = tempfile.mktemp(".csv")
-        file = open(fileName,"w")
-        file.write("Lati,LONG,z \n\
+        '''comma delimited'''
+
+        fileName = tempfile.mktemp('.csv')
+        file = open(fileName, 'w')
+        file.write('Lati,LONG,z \n\
 -25.0,180.0,452.688000\n\
--34,150.0,459.126000\n")
+-34,150.0,459.126000\n')
         file.close()
+
         try:
             results = Geospatial_data(fileName)
@@ -1646 +1507 @@
         else:
             msg = 'Different zones in Geo references not caught.'
-            raise msg        
-        
-        os.remove(fileName)
-        
+            raise Exception, msg
+
+        os.remove(fileName)
+
     def test_lat_long(self):
-        lat_gong = degminsec2decimal_degrees(-34,30,0.)
-        lon_gong = degminsec2decimal_degrees(150,55,0.)
-        
-        lat_2 = degminsec2decimal_degrees(-34,00,0.)
-        lon_2 = degminsec2decimal_degrees(150,00,0.)
-        
+        lat_gong = degminsec2decimal_degrees(-34, 30, 0.)
+        lon_gong = degminsec2decimal_degrees(150, 55, 0.)
+
+        lat_2 = degminsec2decimal_degrees(-34, 00, 0.)
+        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
+
         lats = [lat_gong, lat_2]
         longs = [lon_gong, lon_2]
@@ -1662 +1523 @@
 
         points = gsd.get_data_points(absolute=True)
-        
+
         assert num.allclose(points[0][0], 308728.009)
         assert num.allclose(points[0][1], 6180432.601)
-        assert num.allclose(points[1][0],  222908.705)
+        assert num.allclose(points[1][0], 222908.705)
         assert num.allclose(points[1][1], 6233785.284)
         self.failUnless(gsd.get_geo_reference().get_zone() == 56,
                         'Bad zone error!')
-        
+
+        # use self.failUnlessRaises(ValueError, Geospatial_data(latitudes=lats))
         try:
             results = Geospatial_data(latitudes=lats)
@@ -1675 +1537 @@
             pass
         else:
-            self.failUnless(0 ==1,  'Error not thrown error!')
+            self.fail('Error not thrown error!')
         try:
             results = Geospatial_data(latitudes=lats)
@@ -1681 +1543 @@
             pass
         else:
-            self.failUnless(0 ==1,  'Error not thrown error!')
+            self.fail('Error not thrown error!')
         try:
             results = Geospatial_data(longitudes=lats)
@@ -1687 +1549 @@
             pass
         else:
-            self.failUnless(0 ==1, 'Error not thrown error!')
+            self.fail('Error not thrown error!')
         try:
             results = Geospatial_data(latitudes=lats, longitudes=longs,
@@ -1694 +1556 @@
             pass
         else:
-            self.failUnless(0 ==1,  'Error not thrown error!')
-            
+            self.fail('Error not thrown error!')
+
         try:
             results = Geospatial_data(latitudes=lats, longitudes=longs,
@@ -1702 +1564 @@
             pass
         else:
-            self.failUnless(0 ==1,  'Error not thrown error!')
+            self.fail('Error not thrown error!')
 
     def test_lat_long2(self):
-        lat_gong = degminsec2decimal_degrees(-34,30,0.)
-        lon_gong = degminsec2decimal_degrees(150,55,0.)
-        
-        lat_2 = degminsec2decimal_degrees(-34,00,0.)
-        lon_2 = degminsec2decimal_degrees(150,00,0.)
-        
+        lat_gong = degminsec2decimal_degrees(-34, 30, 0.)
+        lon_gong = degminsec2decimal_degrees(150, 55, 0.)
+
+        lat_2 = degminsec2decimal_degrees(-34, 00, 0.)
+        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
+
         points = [[lat_gong, lon_gong], [lat_2, lon_2]]
         gsd = Geospatial_data(data_points=points, points_are_lats_longs=True)
 
         points = gsd.get_data_points(absolute=True)
-        
+
         assert num.allclose(points[0][0], 308728.009)
         assert num.allclose(points[0][1], 6180432.601)
-        assert num.allclose(points[1][0],  222908.705)
+        assert num.allclose(points[1][0], 222908.705)
         assert num.allclose(points[1][1], 6233785.284)
         self.failUnless(gsd.get_geo_reference().get_zone() == 56,
@@ -1728 +1590 @@
             pass
         else:
-            self.failUnless(0 ==1,  'Error not thrown error!')
-
+            self.fail('Error not thrown error!')
 
     def test_write_urs_file(self):
-        lat_gong = degminsec2decimal_degrees(-34,00,0)
-        lon_gong = degminsec2decimal_degrees(150,30,0.)
-        
-        lat_2 = degminsec2decimal_degrees(-34,00,1)
-        lon_2 = degminsec2decimal_degrees(150,00,0.)
+        lat_gong = degminsec2decimal_degrees(-34, 00, 0)
+        lon_gong = degminsec2decimal_degrees(150, 30, 0.)
+
+        lat_2 = degminsec2decimal_degrees(-34, 00, 1)
+        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
         p1 = (lat_gong, lon_gong)
         p2 = (lat_2, lon_2)
@@ -1742 +1603 @@
         gsd = Geospatial_data(data_points=list(points),
                               points_are_lats_longs=True)
-        
-        fn = tempfile.mktemp(".urs")
+
+        fn = tempfile.mktemp('.urs')
         gsd.export_points_file(fn)
-        #print "fn", fn 
         handle = open(fn)
         lines = handle.readlines()
-        assert lines[0],'2'
-        assert lines[1],'-34.0002778 150.0 0'
-        assert lines[2],'-34.0 150.5 1'
+        assert lines[0], '2'
+        assert lines[1], '-34.0002778 150.0 0'
+        assert lines[2], '-34.0 150.5 1'
         handle.close()
         os.remove(fn)
-        
+
     def test_lat_long_set(self):
-        lat_gong = degminsec2decimal_degrees(-34,30,0.)
-        lon_gong = degminsec2decimal_degrees(150,55,0.)
-        
-        lat_2 = degminsec2decimal_degrees(-34,00,0.)
-        lon_2 = degminsec2decimal_degrees(150,00,0.)
+        lat_gong = degminsec2decimal_degrees(-34, 30, 0.)
+        lon_gong = degminsec2decimal_degrees(150, 55, 0.)
+
+        lat_2 = degminsec2decimal_degrees(-34, 00, 0.)
+        lon_2 = degminsec2decimal_degrees(150, 00, 0.)
         p1 = (lat_gong, lon_gong)
         p2 = (lat_2, lon_2)
@@ -1767 +1627 @@
 
         points = gsd.get_data_points(absolute=True)
-        #print "points[0][0]", points[0][0]
         #Note the order is unknown, due to using sets
         # and it changes from windows to linux
@@ -1773 +1632 @@
             assert num.allclose(points[1][0], 308728.009)
             assert num.allclose(points[1][1], 6180432.601)
-            assert num.allclose(points[0][0],  222908.705)
+            assert num.allclose(points[0][0], 222908.705)
             assert num.allclose(points[0][1], 6233785.284)
         except AssertionError:
             assert num.allclose(points[0][0], 308728.009)
             assert num.allclose(points[0][1], 6180432.601)
-            assert num.allclose(points[1][0],  222908.705)
+            assert num.allclose(points[1][0], 222908.705)
             assert num.allclose(points[1][1], 6233785.284)
-            
+
         self.failUnless(gsd.get_geo_reference().get_zone() == 56,
                         'Bad zone error!')
         points = gsd.get_data_points(as_lat_long=True)
-        #print "test_lat_long_set points", points
         try:
             assert num.allclose(points[0][0], -34)
@@ -1793 +1651 @@
 
     def test_len(self):
-        
         points = [[1.0, 2.1], [3.0, 5.3]]
         G = Geospatial_data(points)
-        self.failUnless(2 ==len(G),  'Len error!')
-        
+        self.failUnless(2 == len(G), 'Len error!')
+
         points = [[1.0, 2.1]]
         G = Geospatial_data(points)
-        self.failUnless(1 ==len(G),  'Len error!')
+        self.failUnless(1 == len(G), 'Len error!')
 
         points = [[1.0, 2.1], [3.0, 5.3], [3.0, 5.3], [3.0, 5.3]]
         G = Geospatial_data(points)
-        self.failUnless(4 ==len(G),  'Len error!')
-        
+        self.failUnless(4 == len(G), 'Len error!')
+
     def test_split(self):
         """test if the results from spilt are disjoin sets"""
-        
-        #below is a work around until the randint works on cyclones compute nodes
-        if get_host_name()[8:9]!='0':
-                
-            
+
+        # below is a workaround until randint works on cyclones compute nodes
+        if get_host_name()[8:9] != '0':
             points = [[1.0, 1.0], [1.0, 2.0],[1.0, 3.0], [1.0, 4.0], [1.0, 5.0],
                       [2.0, 1.0], [2.0, 2.0],[2.0, 3.0], [2.0, 4.0], [2.0, 5.0],
@@ -1818 +1673 @@
                       [4.0, 1.0], [4.0, 2.0],[4.0, 3.0], [4.0, 4.0], [4.0, 5.0],
                       [5.0, 1.0], [5.0, 2.0],[5.0, 3.0], [5.0, 4.0], [5.0, 5.0]]
-            attributes = {'depth':[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 
-                          14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25],
-                          'speed':[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 
-                          14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]}
+            attributes = {'depth': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
+                                    11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
+                                    21, 22, 23, 24, 25],
+                          'speed': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
+                                    11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
+                                    21, 22, 23, 24, 25]}
             G = Geospatial_data(points, attributes)
-    
+
             factor = 0.21
-    
-            #will return G1 with 10% of points and G2 with 90%
-            G1, G2  = G.split(factor,100) 
-            
+
+            # will return G1 with 10% of points and G2 with 90%
+            G1, G2  = G.split(factor, 100)
             assert num.allclose(len(G), len(G1)+len(G2))
             assert num.allclose(round(len(G)*factor), len(G1))
-    
+
             P = G1.get_data_points(absolute=False)
-            assert num.allclose(P, [[5.0,4.0],[4.0,3.0],[4.0,2.0],[3.0,1.0],[2.0,3.0]])
-    
+            expected = [[5.,4.], [4.,1.], [2.,4.], [2.,3.], [1.,4.]]
+            #            expected = [[5.0, 4.0], [4.0, 3.0], [4.0, 2.0],
+                                     #            [3.0, 1.0], [2.0, 3.0]]
+            msg = 'Expected %s, but\nP=%s' % (str(expected), str(P))
+            assert num.allclose(P, expected), msg
+
             A = G1.get_attributes()
-            assert num.allclose(A,[24, 18, 17, 11, 8])
-        
+            expected = [24, 16, 9, 8, 4]
+            msg = 'expected=%s, but A=%s' % (str(expected), str(A))
+            assert num.allclose(A, expected), msg
+
     def test_split1(self):
-        """test if the results from spilt are disjoin sets"""
-        #below is a work around until the randint works on cyclones compute nodes
-        if get_host_name()[8:9]!='0':
-
-            from RandomArray import randint,seed
-            seed(100,100)
-            a_points = randint(0,999999,(10,2))
+        """test if the results from split are disjoint sets"""
+
+        # below is a workaround until randint works on cyclones compute nodes
+        if get_host_name()[8:9] != '0':
+            from numpy.random import randint, seed
+
+            seed((100, 100))
+            a_points = randint(0, 999999, (10,2))
             points = a_points.tolist()
-    #        print points
-    
+
             G = Geospatial_data(points)
-    
+
             factor = 0.1
-    
-            #will return G1 with 10% of points and G2 with 90%
-            G1, G2  = G.split(factor,100) 
-            
-    #        print 'G1',G1
+
+            # will return G1 with 10% of points and G2 with 90%
+            G1, G2  = G.split(factor, 100)
             assert num.allclose(len(G), len(G1)+len(G2))
             assert num.allclose(round(len(G)*factor), len(G1))
-    
+
             P = G1.get_data_points(absolute=False)
-            assert num.allclose(P, [[982420.,28233.]])
-
- 
+            expected = [[425835., 137518.]]
+            msg = 'expected=%s, but\nP=%s' % (str(expected), str(P))
+            assert num.allclose(P, expected), msg
+
     def test_find_optimal_smoothing_parameter(self):
         """
-        Creates a elevation file represting hill (sort of) and runs 
+        Creates a elevation file representing a hill (sort of) and runs
         find_optimal_smoothing_parameter for 3 different alphas,
-        
+
         NOTE the random number seed is provided to control the results
         """
+
         from cmath import cos
 
-        #below is a work around until the randint works on cyclones compute nodes
+        # below is a workaround until randint works on cyclones compute nodes
         if get_host_name()[8:9]!='0':
-
-            filename = tempfile.mktemp(".csv")
-            file = open(filename,"w")
-            file.write("x,y,elevation \n")
-    
-            for i in range(-5,6):
-                for j in range(-5,6):
-                    #this equation made surface like a circle ripple
+            filename = tempfile.mktemp('.csv')
+            file = open(filename, 'w')
+            file.write('x,y,elevation \n')
+
+            for i in range(-5, 6):
+                for j in range(-5, 6):
+                    # this equation makes a surface like a circle ripple
                     z = abs(cos(((i*i) + (j*j))*.1)*2)
-    #                print 'x,y,f',i,j,z
-                    file.write("%s, %s, %s\n" %(i, j, z))
-                    
+                    file.write("%s, %s, %s\n" % (i, j, z))
+
             file.close()
-     
-            value, alpha = find_optimal_smoothing_parameter(data_file=filename, 
+
+            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
                                                  alpha_list=[0.0001, 0.01, 1],
                                                  mesh_file=None,
@@ -1900 +1760 @@
                                                  seed_num=100000,
                                                  verbose=False)
-    
+
             os.remove(filename)
-            
+
             # print value, alpha
-            assert (alpha==0.01)
+            assert(alpha == 0.01)
 
     def test_find_optimal_smoothing_parameter1(self):
         """
-        Creates a elevation file represting hill (sort of) and
+        Creates a elevation file representing  a hill (sort of) and
         Then creates a mesh file and passes the mesh file and the elevation
         file to find_optimal_smoothing_parameter for 3 different alphas,
-        
+
         NOTE the random number seed is provided to control the results
         """
-        #below is a work around until the randint works on cyclones compute nodes
+
+        # below is a workaround until randint works on cyclones compute nodes
         if get_host_name()[8:9]!='0':
-
             from cmath import cos
             from anuga.pmesh.mesh_interface import create_mesh_from_regions
-            
-            filename = tempfile.mktemp(".csv")
-            file = open(filename,"w")
-            file.write("x,y,elevation \n")
-    
-            for i in range(-5,6):
-                for j in range(-5,6):
-                    #this equation made surface like a circle ripple
+
+            filename = tempfile.mktemp('.csv')
+            file = open(filename, 'w')
+            file.write('x,y,elevation \n')
+
+            for i in range(-5 ,6):
+                for j in range(-5, 6):
+                    # this equation makes a surface like a circle ripple
                     z = abs(cos(((i*i) + (j*j))*.1)*2)
-    #                print 'x,y,f',i,j,z
-                    file.write("%s, %s, %s\n" %(i, j, z))
-                    
+                    file.write('%s, %s, %s\n' % (i, j, z))
+
             file.close()
-            poly=[[5,5],[5,-5],[-5,-5],[-5,5]]
-            internal_poly=[[[[1,1],[1,-1],[-1,-1],[-1,1]],.5]]
-            mesh_filename= tempfile.mktemp(".msh")
-            
+
+            poly=[[5,5], [5,-5], [-5,-5], [-5,5]]
+            internal_poly=[[[[1,1], [1,-1], [-1,-1], [-1,1]], .5]]
+            mesh_filename= tempfile.mktemp('.msh')
+
             create_mesh_from_regions(poly,
-                                 boundary_tags={'back': [2],
-                                                'side': [1,3],
-                                                'ocean': [0]},
-                             maximum_triangle_area=3,
-                             interior_regions=internal_poly,
-                             filename=mesh_filename,
-                             use_cache=False,
-                             verbose=False)
-     
-            value, alpha = find_optimal_smoothing_parameter(data_file=filename, 
+                                     boundary_tags={'back': [2],
+                                                    'side': [1,3],
+                                                    'ocean': [0]},
+                                     maximum_triangle_area=3,
+                                     interior_regions=internal_poly,
+                                     filename=mesh_filename,
+                                     use_cache=False,
+                                     verbose=False)
+
+            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
                                                  alpha_list=[0.0001, 0.01, 1],
                                                  mesh_file=mesh_filename,
@@ -1952 +1812 @@
                                                  seed_num=174,
                                                  verbose=False)
-    
+
             os.remove(filename)
             os.remove(mesh_filename)
-            
-    #        print value, alpha
-            assert (alpha==0.01)
+
+            msg = 'alpha=%s' % str(alpha)
+            # 0.01 was expected with Numeric.RandomArray RNG
+            assert alpha==1.0, msg
 
     def test_find_optimal_smoothing_parameter2(self):
-        """
-        Tests requirement that mesh file must exist or IOError is thrown
-        
+        '''Tests requirement that mesh file must exist or IOError is thrown
+
         NOTE the random number seed is provided to control the results
-        """
+        '''
+
         from cmath import cos
         from anuga.pmesh.mesh_interface import create_mesh_from_regions
-        
-        filename = tempfile.mktemp(".csv")
-        mesh_filename= tempfile.mktemp(".msh")
-        
+
+        filename = tempfile.mktemp('.csv')
+        mesh_filename= tempfile.mktemp('.msh')
+
         try:
-            value, alpha = find_optimal_smoothing_parameter(data_file=filename, 
+            value, alpha = find_optimal_smoothing_parameter(data_file=filename,
                                              alpha_list=[0.0001, 0.01, 1],
                                              mesh_file=mesh_filename,
@@ -1981 +1842 @@
             pass
         else:
-            self.failUnless(0 ==1,  'Error not thrown error!')
-        
-         
+            self.fail('Error not thrown error!')
+
+################################################################################
+
 if __name__ == "__main__":
-
-    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_write_csv_attributes_lat_long')
-    suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter')
-    #suite = unittest.makeSuite(Test_Geospatial_data, 'test_split1')
-    #suite = unittest.makeSuite(Test_Geospatial_data, 'test')
+    suite = unittest.makeSuite(Test_Geospatial_data, 'test')
     runner = unittest.TextTestRunner() #verbosity=2)
     runner.run(suite)
 
-    
+