source: branches/numpy/anuga/geospatial_data/geospatial_data.py @ 6517

"""Class Geospatial_data

Manipulation of locations on the planet and associated attributes.
"""

from sys import maxint
from os import access, F_OK, R_OK,remove
from types import DictType
from warnings import warn
from string import lower
from RandomArray import randint, seed, get_seed
from copy import deepcopy

from Scientific.IO.NetCDF import NetCDFFile
import numpy as num

from anuga.coordinate_transforms.lat_long_UTM_conversion import UTMtoLL
from anuga.utilities.numerical_tools import ensure_numeric
from anuga.coordinate_transforms.geo_reference import Geo_reference, \
     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'


##
# @brief ??
class Geospatial_data:

    ##
    # @brief 
    # @param data_points Mx2 iterable of tuples or array of x,y coordinates.
    # @param attributes Associated values for each data point.
    # @param geo_reference ??
    # @param default_attribute_name ??
    # @param file_name 
    # @param latitudes ??
    # @param longitudes ??
    # @param points_are_lats_longs True if points are lat/long, not UTM.
    # @param max_read_lines Size of block to read, if blocking.
    # @param load_file_now True if blocking but we eant to read file now.
    # @param verbose True if this class instance is verbose.
    def __init__(self,
                 data_points=None, # this can also be a points file name
                 attributes=None,
                 geo_reference=None,
                 default_attribute_name=None,
                 file_name=None,
                 latitudes=None,
                 longitudes=None,
                 points_are_lats_longs=False,
                 max_read_lines=None,
                 load_file_now=True,
                 verbose=False):
        """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
        with extension .txt, .cvs or .pts can also be passed in here.

        attributes: Associated values for each data point. The type
        must be either a list or an array of length M or a dictionary
        of lists (or arrays) of length M. In the latter case the keys
        in the dictionary represent the attribute names, in the former
        the attribute will get the default name "elevation".

        geo_reference: Object representing the origin of the data
        points. It contains UTM zone, easting and northing and data
        points are assumed to be relative to this origin.
        If geo_reference is None, the default geo ref object is used.

        default_attribute_name: Name of default attribute to be used with
        get_attribute_values. The idea is that the dataset can be
        equipped with information about which attribute to return.
        If None, the default is the "first"

        latitudes, longitudes: Vectors of latitudes and longitudes,
        used to specify location instead of points.

        points_are_lats_longs: Set this as true if the points are actually
        lats and longs, not UTM

        max_read_lines: The number of rows read into memory when using
        blocking to read a file.

        load_file_now:  If true the file is automatically loaded
        into the geospatial instance. Used when blocking.

        file_name: Name of input netCDF file or .txt file. netCDF file must
        have dimensions "points" etc.
        .txt file is a comma seperated file with x, y and attribute
        data.

        The first line has the titles of the columns.  The first two
        column titles are checked to see if they start with lat or
        long (not case sensitive).  If so the data is assumed to be
        latitude and longitude, in decimal format and converted to
        UTM.  Otherwise the first two columns are assumed to be the x
        and y, and the title names acually used are ignored.


        The format for a .txt file is:
            1st line:     [column names]
            other lines:  x y [attributes]

            for example:
            x, y, elevation, friction
            0.6, 0.7, 4.9, 0.3
            1.9, 2.8, 5, 0.3
            2.7, 2.4, 5.2, 0.3

        The first two columns have to be x, y or lat, long
        coordinates.


        The format for a Points dictionary is:
          ['pointlist'] a 2 column array describing points. 1st column x,
          2nd column y.
          ['attributelist'], a dictionary of 1D arrays, representing
          attribute values at the point.  The dictionary key is the attribute
          header.
          ['geo_reference'] a Geo_refernece object. Use if the point
          information is relative. This is optional.
            eg
            dic['pointlist'] = [[1.0,2.0],[3.0,5.0]]
            dic['attributelist']['elevation'] = [[7.0,5.0]

        verbose:
        """

        if isinstance(data_points, basestring):
            # assume data_points is really a file name
            file_name = data_points

        self.set_verbose(verbose)
        self.geo_reference = None
        self.file_name = file_name

        if max_read_lines is None:
            self.max_read_lines = MAX_READ_LINES
        else:
            self.max_read_lines = max_read_lines

        if file_name is None:
            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,
                                        longitudes=longitudes,
                                        geo_reference=geo_reference,
                                        data_points=data_points,
                                        points_are_lats_longs=
                                            points_are_lats_longs)
            self.check_data_points(data_points)
            self.set_attributes(attributes)
            self.set_geo_reference(geo_reference)
            self.set_default_attribute_name(default_attribute_name)
        elif load_file_now is True:
            # watch for case where file name and points,
            # attributes etc are provided!!
            # if file name then all provided info will be removed!

            if verbose is True:
                if file_name is not None:
                    print 'Loading Geospatial data from file: %s' % file_name

            self.import_points_file(file_name, verbose=verbose)

            self.check_data_points(self.data_points)
            self.set_attributes(self.attributes)
            self.set_geo_reference(self.geo_reference)
            self.set_default_attribute_name(default_attribute_name)

        if verbose is True:
            if file_name is not None:
                print 'Geospatial data created from file: %s' % file_name
                if load_file_now is False:
                    print 'Data will be loaded blockwise on demand'

                    if file_name.endswith('csv') or file_name.endswith('txt'):
                        pass
                        # 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'

    ##
    # @brief Return length of the points set.
    def __len__(self):
        return len(self.data_points)

    ##
    # @brief Return a string representation of the points set.
    def __repr__(self):
        return str(self.get_data_points(absolute=True))

    ##
    # @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"""

        if data_points is None:
            self.data_points = None
            msg = 'There is no data or file provided!'
            raise ValueError, msg
        else:
            self.data_points = ensure_numeric(data_points)
            if not (0,) == self.data_points.shape:
                assert len(self.data_points.shape) == 2
                assert self.data_points.shape[1] == 2

    ##
    # @brief Check and assign attributes data.
    # @param attributes Dictionary or scalar to save as .attributes.
    # @note Throws exception if unable to convert dict keys to numeric.
    def set_attributes(self, attributes):
        """Check and assign attributes dictionary"""

        if attributes is None:
            self.attributes = None
            return

        if not isinstance(attributes, DictType):
            # Convert single attribute into dictionary
            attributes = {DEFAULT_ATTRIBUTE: attributes}

        # Check input attributes
        for key in attributes.keys():
            try:
                attributes[key] = ensure_numeric(attributes[key])
            except:
                msg = ("Attribute '%s' (%s) could not be converted to a"
                       "numeric vector" % (str(key), str(attributes[key])))
                raise Exception, msg

        self.attributes = attributes

    ##
    # @brief Set the georeference of geospatial data.
    # @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):
        """Set the georeference of geospatial data.

        It can also be used to change the georeference and will ensure that
        the absolute coordinate values are unchanged.
        """

        if geo_reference is None:
            # Use default - points are in absolute coordinates
            geo_reference = Geo_reference()

        # Allow for tuple (zone, xllcorner, yllcorner)
        geo_reference = ensure_geo_reference(geo_reference)

        if not isinstance(geo_reference, Geo_reference):
            # 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 '
                   'object or None.')
            raise Expection, msg

        # If a geo_reference already exists, change the point data according to
        # the new geo reference
        if  self.geo_reference is not None:
            self.data_points = self.get_data_points(geo_reference=geo_reference)

        self.geo_reference = geo_reference

    ##
    # @brief Set default attribute name.
    # @param default_attribute_name The default to save.
    def set_default_attribute_name(self, default_attribute_name):
        self.default_attribute_name = default_attribute_name

    ##
    # @brief Set the instance verbose flag.
    # @param verbose The value to save.
    # @note Will raise exception if param is not True or False.
    def set_verbose(self, verbose=False):
        if verbose in [False, True]:
            self.verbose = verbose
        else:
            msg = 'Illegal value: %s' % str(verbose)
            raise Exception, msg

    ##
    # @brief Clip geospatial data by a given polygon.
    # @param polygon The polygon to clip with.
    # @param closed True if points on clip boundary are not included in result.
    # @param verbose True if this function is verbose.
    def clip(self, polygon, closed=True, verbose=False):
        """Clip geospatial data by a polygon

        Input
          polygon - Either a list of points, an Nx2 array or
                    a Geospatial data object.
          closed - (optional) determine whether points on boundary should be
          regarded as belonging to the polygon (closed = True)
          or not (closed = False). Default is True.

        Output
          New geospatial data object representing points inside
          specified polygon.


        Note - this method is non-destructive and leaves the data in 'self'
        unchanged
        """

        from anuga.utilities.polygon import inside_polygon

        if isinstance(polygon, Geospatial_data):
            # Polygon is an object - extract points
            polygon = polygon.get_data_points()

        points = self.get_data_points()
        inside_indices = inside_polygon(points, polygon, closed, verbose)

        clipped_G = self.get_sample(inside_indices)

        return clipped_G

    ##
    # @brief Clip points data by polygon, return points outside polygon.
    # @param polygon The polygon to clip with.
    # @param closed True if points on clip boundary are not included in result.
    # @param verbose True if this function is verbose.
    def clip_outside(self, polygon, closed=True, verbose=False):
        """Clip geospatial date by a polygon, keeping data OUTSIDE of polygon

        Input
          polygon - Either a list of points, an Nx2 array or
                    a Geospatial data object.
          closed - (optional) determine whether points on boundary should be
          regarded as belonging to the polygon (closed = True)
          or not (closed = False). Default is True.

        Output
          Geospatial data object representing point OUTSIDE specified polygon
        """

        from anuga.utilities.polygon import outside_polygon

        if isinstance(polygon, Geospatial_data):
            # Polygon is an object - extract points
            polygon = polygon.get_data_points()

        points = self.get_data_points()
        outside_indices = outside_polygon(points, polygon, closed,verbose)

        clipped_G = self.get_sample(outside_indices)

        return clipped_G

    ##
    # @brief Get instance geo_reference data.
    def get_geo_reference(self):
        return self.geo_reference

    ##
    # @brief Get coordinates for all data points as an Nx2 array.
    # @param absolute If True, return UTM, else relative to xll/yll corners.
    # @param geo_reference If supplied, points are relative to it.
    # @param as_lat_long If True, return points as lat/lon.
    # @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):
        """Get coordinates for all data points as an Nx2 array

        If absolute is False returned coordinates are relative to the
        internal georeference's xll and yll corners, otherwise
        absolute UTM coordinates are returned.

        If a geo_reference is passed the points are returned relative
        to that geo_reference.

        isSH (isSouthHemisphere) is only used when getting data
        points "as_lat_long" is True and if FALSE will return lats and
        longs valid for the Northern Hemisphere.

        Default: absolute is True.
        """

        if as_lat_long is True:
            msg = "Points need a zone to be converted into lats and longs"
            assert self.geo_reference is not None, msg
            zone = self.geo_reference.get_zone()
            assert self.geo_reference.get_zone() is not DEFAULT_ZONE, msg
            lats_longs = []
            for point in self.get_data_points(True):
                # UTMtoLL(northing, easting, zone,
                lat_calced, long_calced = UTMtoLL(point[1], point[0],
                                                  zone, isSouthHemisphere)
                lats_longs.append((lat_calced, long_calced)) # to hash
            return lats_longs

        if absolute is True and geo_reference is None:
            return self.geo_reference.get_absolute(self.data_points)
        elif geo_reference is not None:
            return geo_reference.change_points_geo_ref(self.data_points,
                                                       self.geo_reference)
        else:
            # If absolute is False
            return self.data_points

    ##
    # @brief Get value for attribute name.
    # @param attribute_name Name to get value for.
    # @note If name passed is None, return default attribute value.
    def get_attributes(self, attribute_name=None):
        """Return values for one named attribute.

        If attribute_name is None, default_attribute_name is used
        """

        if attribute_name is None:
            if self.default_attribute_name is not None:
                attribute_name = self.default_attribute_name
            else:
                attribute_name = self.attributes.keys()[0]
                # above line takes the first one from keys

        if self.verbose is True:
            print 'Using attribute %s' %attribute_name
            print 'Available attributes: %s' %(self.attributes.keys())

        msg = 'Attribute name %s does not exist in data set' % attribute_name
        assert self.attributes.has_key(attribute_name), msg

        return self.attributes[attribute_name]

    ##
    # @brief Get all instance attributes.
    # @return The instance attribute dictionary, or None if no attributes.
    def get_all_attributes(self):
        """Return values for all attributes.
        The return value is either None or a dictionary (possibly empty).
        """

        return self.attributes

    ##
    # @brief Override __add__() to allow addition of geospatial objects.
    # @param self This object.
    # @param other The second object.
    # @return The new geospatial object.
    def __add__(self, other):
        """Returns the addition of 2 geospatial objects,
        objects are concatencated to the end of each other

        NOTE: doesn't add if objects contain different
        attributes

        Always return absolute points!
        This also means, that if you add None to the object,
        it will be turned into absolute coordinates

        other can be None in which case nothing is added to self.
        """

        # find objects zone and checks if the same
        geo_ref1 = self.get_geo_reference()
        zone1 = geo_ref1.get_zone()

        if other is not None:
            geo_ref2 = other.get_geo_reference()
            zone2 = geo_ref2.get_zone()
            geo_ref1.reconcile_zones(geo_ref2)
            new_points = num.concatenate((self.get_data_points(absolute=True),
                                          other.get_data_points(absolute=True)),
                                         axis = 0)

            # Concatenate attributes if any
            if self.attributes is None:
                if other.attributes is not None:
                    msg = ('Geospatial data must have the same '
                           'attributes to allow addition.')
                    raise Exception, msg

                new_attributes = None
            else:
                new_attributes = {}
                for x in self.attributes.keys():
                    if other.attributes.has_key(x):
                        attrib1 = self.attributes[x]
                        attrib2 = other.attributes[x]
                        new_attributes[x] = num.concatenate((attrib1, attrib2),
                                                            axis=0) #??default#
                    else:
                        msg = ('Geospatial data must have the same '
                               'attributes to allow addition.')
                        raise Exception, msg
        else:
            # other is None:
            new_points = self.get_data_points(absolute=True)
            new_attributes = self.attributes

        # Instantiate new data object and return absolute coordinates
        new_geo_ref = Geo_reference(geo_ref1.get_zone(), 0.0, 0.0)
        return Geospatial_data(new_points, new_attributes, new_geo_ref)

    ##
    # @brief Override the addition case where LHS isn't geospatial object.
    # @param self This object.
    # @param other The second object.
    # @return The new geospatial object.
    def __radd__(self, other):
        """Handle cases like None + Geospatial_data(...)"""

        return self + other

################################################################################
#  IMPORT/EXPORT POINTS FILES
################################################################################

    ##
    # @brief Import a .txt, .csv or .pts points data file.
    # @param file_name 
    # @param delimiter 
    # @param verbose True if this function is to be verbose.
    # @note Will throw IOError or SyntaxError if there is a problem.
    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!

        Post condition: self.attributes dictionary has been set
        """

        if access(file_name, F_OK) == 0 :
            msg = 'File %s does not exist or is not accessible' % file_name
            raise IOError, msg

        attributes = {}
        if file_name[-4:] == ".pts":
            try:
                data_points, attributes, geo_reference = \
                             _read_pts_file(file_name, verbose)
            except IOError, e:
                msg = 'Could not open file %s ' % file_name
                raise IOError, msg
        elif file_name[-4:] == ".txt" or file_name[-4:]== ".csv":
            try:
                data_points, attributes, geo_reference = \
                             _read_csv_file(file_name, verbose)
            except IOError, e:
                # This should only be if a file is not found
                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 = ('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:]
            raise IOError, msg

        self.data_points = data_points
        self.attributes = attributes
        self.geo_reference = geo_reference

    ##
    # @brief Write points data to a file (.csv or .pts).
    # @param file_name Path to file to write.
    # @param absolute ??
    # @param as_lat_long ??
    # @param isSouthHemisphere ??
    def export_points_file(self, file_name, absolute=True,
                           as_lat_long=False, isSouthHemisphere=True):
        """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.

        If absolute is True data the xll and yll are added to the points value
        and the xll and yll of the geo_reference are set to 0.

        If absolute is False data points at returned as relative to the xll
        and yll and geo_reference remains uneffected

        isSouthHemisphere: is only used when getting data
        points "as_lat_long" is True and if FALSE will return lats and
        longs valid for the Northern Hemisphere.
        """

        if (file_name[-4:] == ".pts"):
            if absolute is True:
                geo_ref = deepcopy(self.geo_reference)
                geo_ref.xllcorner = 0
                geo_ref.yllcorner = 0
                _write_pts_file(file_name,
                                self.get_data_points(absolute),
                                self.get_all_attributes(),
                                geo_ref)
            else:
                _write_pts_file(file_name,
                                self.get_data_points(absolute),
                                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."
            assert absolute, msg
            _write_csv_file(file_name,
                            self.get_data_points(absolute=True,
                                                 as_lat_long=as_lat_long,
                                           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."
            assert absolute, msg
            _write_urs_file(file_name,
                            self.get_data_points(as_lat_long=True,
                                           isSouthHemisphere=isSouthHemisphere))
        else:
            msg = 'Unknown file type %s ' %file_name
            raise IOError, msg

    ##
    # @brief Get a subset of data that is referred to by 'indices'.
    # @param indices A list of indices to select data subset with.
    # @return A geospatial object containing data subset.
    def get_sample(self, indices):
        """ Returns a object which is a subset of the original
        and the data points and attributes in this new object refer to
        the indices provided

        Input
            indices- a list of integers that represent the new object
        Output
            New geospatial data object representing points specified by
            the indices
        """

        # FIXME: add the geo_reference to this
        points = self.get_data_points()
        sampled_points = num.take(points, indices, axis=0)

        attributes = self.get_all_attributes()

        sampled_attributes = {}
        if attributes is not None:
            for key, att in attributes.items():
                sampled_attributes[key] = num.take(att, indices, axis=0)

        return Geospatial_data(sampled_points, sampled_attributes)

    ##
    # @brief Split one geospatial object into two.
    # @param factor Relative size to make first result object.
    # @param seed_num Random 'seed' - used only for unit test.
    # @param verbose True if this function is to be verbose.
    # @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'
        smaller the original and the second is the remainder. The two
        new objects are disjoint sets of each other.

        Points of the two new object have selected RANDOMLY.

        This method create two lists of indices which are passed into
        get_sample.  The lists are created using random numbers, and
        they are unique sets eg.  total_list(1,2,3,4,5,6,7,8,9)
        random_list(1,3,6,7,9) and remainder_list(0,2,4,5,8)

        Input -  the factor which to split the object, if 0.1 then 10% of the
                 together object will be returned

        Output - two geospatial_data objects that are disjoint sets of the
                 original
        """

        i = 0
        self_size = len(self)
        random_list = []
        remainder_list = []
        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_list = total.tolist()

        if verbose: print "total list len", len(total_list)

        # There will be repeated random numbers however will not be a
        # problem as they are being 'pop'ed out of array so if there
        # are two numbers the same they will pop different indicies,
        # still basically random
        ## create list of non-unquie random numbers
        if verbose: print "create random numbers list %s long" %new_size

        # Set seed if provided, mainly important for unit test!
        # plus recalcule seed when no seed provided.
        if seed_num is not None:
            seed(seed_num, seed_num)
        else:
            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
        random_num.sort()
        random_num.reverse()

        if verbose: print "make random number list and get indices"

        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)
        random_num_len = len(random_num)
        for i in random_num:
            random_list.append(remainder_list.pop(i))
            j += 1
            # prints progress
            if verbose and round(random_num_len/10*k) == j:
                print '(%s/%s)' % (j, random_num_len)
                k += 1

        # FIXME: move to tests, it might take a long time
        # then create an array of random length between 500 and 1000,
        # and use a random factor between 0 and 1
        # setup for assertion
        test_total = random_list[:]
        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')
        assert total_list == test_total, msg

        # Get new samples
        if verbose: print "get values of indices for random list"
        G1 = self.get_sample(random_list)
        if verbose: print "get values of indices for opposite of random list"
        G2 = self.get_sample(remainder_list)

        return G1, G2

    ##
    # @brief Allow iteration over this object.
    def __iter__(self):
        """Read in the header, number_of_points and save the
        file pointer position
        """

        # FIXME - what to do if the file isn't there

        # FIXME (Ole): Shouldn't this go into the constructor?
        # This method acts like the constructor when blocking.
        # ... and shouldn't it be called block_size?
        #
        if self.max_read_lines is None:
            self.max_read_lines = MAX_READ_LINES

        if self.file_name[-4:] == ".pts":
            # See if the file is there.  Throw a QUIET IO error if it isn't
            fd = open(self.file_name,'r')
            fd.close()

            # Throws prints to screen if file not present
            self.fid = NetCDFFile(self.file_name, netcdf_mode_r)

            (self.blocking_georef,
             self.blocking_keys,
             self.number_of_points) = _read_pts_file_header(self.fid,
                                                            self.verbose)
            self.start_row = 0
            self.last_row = self.number_of_points
            self.show_verbose = 0
            self.verbose_block_size = (self.last_row + 10)/10
            self.block_number = 0
            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 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)
        else:
            # Assume the file is a csv file
            file_pointer = open(self.file_name)
            self.header, self.file_pointer = _read_csv_file_header(file_pointer)
            self.blocking_georef = None # Used for reconciling zones

        return self

    ##
    # @brief Read another block into the instance.
    def next(self):
        """read a block, instanciate a new geospatial and return it"""

        if self.file_name[-4:] == ".pts":
            if self.start_row == self.last_row:
                # Read the end of the file last iteration
                # Remove blocking attributes
                self.fid.close()
                del self.max_read_lines
                del self.blocking_georef
                del self.last_row
                del self.start_row
                del self.blocking_keys
                del self.fid
                raise StopIteration
            fin_row = self.start_row + self.max_read_lines
            if fin_row > self.last_row:
                fin_row = self.last_row

            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))

                    self.show_verbose += max(self.max_read_lines,
                                             self.verbose_block_size)

            # Read next block
            pointlist, att_dict, = _read_pts_file_blocking(self.fid,
                                                           self.start_row,
                                                           fin_row,
                                                           self.blocking_keys)

            geo = Geospatial_data(pointlist, att_dict, self.blocking_georef)
            self.start_row = fin_row

            self.block_number += 1
        else:
            # Assume the file is a csv file
            try:
                (pointlist,
                 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)

                # Check that the zones haven't changed.
                if geo_ref is not None:
                    geo_ref.reconcile_zones(self.blocking_georef)
                    self.blocking_georef = geo_ref
                elif self.blocking_georef is not None:
                    msg = ('Geo reference given, then not given.'
                           ' This should not happen.')
                    raise ValueError, msg
                geo = Geospatial_data(pointlist, att_dict, geo_ref)
            except StopIteration:
                self.file_pointer.close()
                del self.header
                del self.file_pointer
                raise StopIteration
            except ANUGAError:
                self.file_pointer.close()
                del self.header
                del self.file_pointer
                raise
            except SyntaxError:
                self.file_pointer.close()
                del self.header
                del self.file_pointer
                # This should only be if there is a format error
                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'
                 '        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')

##
# @brief ??
# @param latitudes ??
# @param longitudes ??
# @param geo_reference ??
# @param data_points ??
# @param points_are_lats_longs ??
# @note IS THIS USED???
def _set_using_lat_long(latitudes,
                        longitudes,
                        geo_reference,
                        data_points,
                        points_are_lats_longs):
    """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!')
        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.')
        raise ValueError, msg

    if points_are_lats_longs:
        if data_points is None:
            msg = "Data points are not specified."
            raise ValueError, msg
        lats_longs = ensure_numeric(data_points)
        latitudes = num.ravel(lats_longs[:,0:1])
        longitudes = num.ravel(lats_longs[:,1:])

    if latitudes is None and longitudes is None:
        msg = "Latitudes and Longitudes are not specified."
        raise ValueError, msg

    if latitudes is None:
        msg = "Longitudes are specified yet latitudes aren't."
        raise ValueError, msg

    if longitudes is None:
        msg = "Latitudes are specified yet longitudes aren't."
        raise ValueError, msg

    data_points, zone  = convert_from_latlon_to_utm(latitudes=latitudes,
                                                    longitudes=longitudes)
    return data_points, Geo_reference(zone=zone)


##
# @brief Read a .pts data file.
# @param file_name Path to file to read.
# @param verbose True if this function is to be verbose.
# @return (pointlist, attributes, geo_reference)
def _read_pts_file(file_name, verbose=False):
    """Read .pts NetCDF file

    Return a (dict_points, dict_attribute, geo_ref)
    eg
    dict['points'] = [[1.0,2.0],[3.0,5.0]]
    dict['attributelist']['elevation'] = [[7.0,5.0]]
    """

    if verbose: print 'Reading ', file_name

    # See if the file is there.  Throw a QUIET IO error if it isn't
    fd = open(file_name,'r')
    fd.close()

    # Throws prints to screen if file not present
    fid = NetCDFFile(file_name, netcdf_mode_r)

    pointlist = num.array(fid.variables['points'])
    keys = fid.variables.keys()

    if verbose: print 'Got %d variables: %s' % (len(keys), keys)

    try:
        keys.remove('points')
    except IOError, e:
        fid.close()
        msg = "Expected keyword 'points' but could not find it"
        raise IOError, msg

    attributes = {}
    for key in keys:
        if verbose: print "reading attribute '%s'" % key

        attributes[key] = num.array(fid.variables[key])

    try:
        geo_reference = Geo_reference(NetCDFObject=fid)
    except AttributeError, e:
        geo_reference = None

    fid.close()

    return pointlist, attributes, geo_reference


##
# @brief Read a .csv data file.
# @param file_name Path to the .csv file to read.
# @param verbose True if this function is to be verbose.
def _read_csv_file(file_name, verbose=False):
    """Read .csv file

    Return a dic of array of points, and dic of array of attribute
    eg
    dic['points'] = [[1.0,2.0],[3.0,5.0]]
    dic['attributelist']['elevation'] = [[7.0,5.0]]
    """

    file_pointer = open(file_name)
    header, file_pointer = _read_csv_file_header(file_pointer)
    try:
        (pointlist,
         att_dict,
         geo_ref,
         file_pointer) = _read_csv_file_blocking(file_pointer,
                                                 header,
                                                 max_read_lines=1e30)
                                    # If the file is bigger that this, block..
                                    # FIXME (Ole) What's up here?
    except ANUGAError:
        file_pointer.close()
        raise

    file_pointer.close()

    return pointlist, att_dict, geo_ref


##
# @brief Read a .csv file header.
# @param file_pointer Open descriptor of the file to read.
# @param delimiter Header line delimiter string, split on this string.
# @param verbose True if this function is to be verbose.
# @return A tuple of (<cleaned header string>, <input file_pointer>)

CSV_DELIMITER = ','

def _read_csv_file_header(file_pointer,
                          delimiter=CSV_DELIMITER,
                          verbose=False):
    """Read the header of a .csv file
    Return a list of the header names
    """

    line = file_pointer.readline()
    header = clean_line(line, delimiter)

    return header, file_pointer

##
# @brief Read a .csv file, with blocking.
# @param file_pointer Open descriptor of the file to read.
# @param header List of already read .csv header fields.
# @param delimiter Delimiter string header was split on.
# @param max_read_lines The max number of lines to read before blocking.
# @param verbose True if this function is to be verbose.
# @note Will throw IndexError, SyntaxError exceptions.
def _read_csv_file_blocking(file_pointer,
                            header,
                            delimiter=CSV_DELIMITER,
                            max_read_lines=MAX_READ_LINES,
                            verbose=False):
    """Read the body of a .csv file.
    header: The list header of the csv file, with the x and y labels.
    """

    points = []
    pointattributes = []
    att_dict = {}

    # This is to remove the x and y headers.
    header = header[:]
    try:
        x_header = header.pop(0)
        y_header = header.pop(0)
    except IndexError:
        # if there are not two columns this will occur.
        # eg if it is a space seperated file
        raise SyntaxError

    read_lines = 0
    while read_lines < max_read_lines:
        line = file_pointer.readline()
        numbers = clean_line(line, delimiter)
        if len(numbers) <= 1:
            break
        if line[0] == '#':
            continue

        read_lines += 1

        try:
            x = float(numbers[0])
            y = float(numbers[1])
            points.append([x,y])
            numbers.pop(0)
            numbers.pop(0)
            if len(header) != len(numbers):
                file_pointer.close()
                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 != '':
                    att_dict.setdefault(header[i],[]).append(float(n))
        except ValueError:
            raise SyntaxError

    if points == []:
        raise StopIteration

    pointlist = num.array(points, num.float)
    for key in att_dict.keys():
        att_dict[key] = num.array(att_dict[key], num.float)

    # Do stuff here so the info is in lat's and longs
    geo_ref = None
    x_header = lower(x_header[:3])
    y_header = lower(y_header[:3])
    if (x_header == 'lon' or  x_header == 'lat') \
       and (y_header == 'lon' or  y_header == 'lat'):
        if x_header == 'lon':
            longitudes = num.ravel(pointlist[:,0:1])
            latitudes = num.ravel(pointlist[:,1:])
        else:
            latitudes = num.ravel(pointlist[:,0:1])
            longitudes = num.ravel(pointlist[:,1:])

        pointlist, geo_ref = _set_using_lat_long(latitudes,
                                                 longitudes,
                                                 geo_reference=None,
                                                 data_points=None,
                                                 points_are_lats_longs=False)

    return pointlist, att_dict, geo_ref, file_pointer


##
# @brief Read a .pts file header.
# @param fid Handle to the open .pts file.
# @param verbose True if the function is to be verbose.
# @return (geo_reference, keys, fid.dimensions['number_of_points'])
# @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'''

    keys = fid.variables.keys()
    try:
        keys.remove('points')
    except IOError, e:
        fid.close()
        msg = "Expected keyword 'points' but could not find it."
        raise IOError, msg

    if verbose: print 'Got %d variables: %s' % (len(keys), keys)

    try:
        geo_reference = Geo_reference(NetCDFObject=fid)
    except AttributeError, e:
        geo_reference = None

    return geo_reference, keys, fid.dimensions['number_of_points']


##
# @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.
# @param fin_row End row index of points to return.
# @param keys Iterable of keys to return.
# @return Tuple of (pointlist, attributes).
def _read_pts_file_blocking(fid, start_row, fin_row, keys):
    '''Read the body of a .pts file.'''

    pointlist = num.array(fid.variables['points'][start_row:fin_row])

    attributes = {}
    for key in keys:
        attributes[key] = num.array(fid.variables[key][start_row:fin_row])

    return pointlist, attributes


##
# @brief Write a .pts data file.
# @param file_name Path to the file to write.
# @param write_data_points Data points to write.
# @param write_attributes Attributes to write.
# @param write_geo_reference Georef to write.
def _write_pts_file(file_name,
                    write_data_points,
                    write_attributes=None,
                    write_geo_reference=None):
    """Write .pts NetCDF file

    NOTE: Below might not be valid ask Duncan : NB 5/2006

    WARNING: This function mangles the point_atts data structure
    # F??ME: (DSG)This format has issues.
    # There can't be an attribute called points
    # consider format change
    # method changed by NB not sure if above statement is correct

    should create new test for this
    legal_keys = ['pointlist', 'attributelist', 'geo_reference']
    for key in point_atts.keys():
        msg = 'Key %s is illegal. Valid keys are %s' %(key, legal_keys)
        assert key in legal_keys, msg
    """

    # NetCDF file definition
    outfile = NetCDFFile(file_name, netcdf_mode_w)

    # Create new file
    outfile.institution = 'Geoscience Australia'
    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

    # Variable definition
    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, netcdf_float, ('number_of_points',))
            outfile.variables[key][:] = write_attributes[key]

    if write_geo_reference is not None:
        write_NetCDF_georeference(write_geo_reference, outfile)

    outfile.close()


##
# @brief Write a .csv data file.
# @param file_name Path to the file to write.
# @param write_data_points Data points to write.
# @param write_attributes Attributes to write.
# @param as_lat_long True if points are lat/lon, else x/y.
# @param delimiter The CSV delimiter to use.
def _write_csv_file(file_name,
                    write_data_points,
                    write_attributes=None,
                    as_lat_long=False,
                    delimiter=','):
    """Write a .csv file."""

    points = write_data_points
    pointattributes = write_attributes

    fd = open(file_name, 'w')

    if as_lat_long:
        titlelist = "latitude" + delimiter + "longitude"  + delimiter
    else:
        titlelist = "x" + delimiter + "y"  + delimiter

    if pointattributes is not None:
        for title in pointattributes.keys():
            titlelist = titlelist + title + delimiter
        titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter

    fd.write(titlelist + "\n")

    # <x/lat> <y/long> [attributes]
    for i, vert in enumerate( points):
        if pointattributes is not None:
            attlist = ","
            for att in pointattributes.keys():
                attlist = attlist + str(pointattributes[att][i]) + delimiter
            attlist = attlist[0:-len(delimiter)] # remove the last delimiter
            attlist.strip()
        else:
            attlist = ''

        fd.write(str(vert[0]) + delimiter + str(vert[1]) + attlist + "\n")

    fd.close()


##
# @brief Write a URS file.
# @param file_name The path of the file to write.
# @param points 
# @param delimiter 
def _write_urs_file(file_name, points, delimiter=' '):
    """Write a URS format file.
    export a file, file_name, with the urs format
    the data points are in lats and longs
    """

    fd = open(file_name, 'w')

    # first line is # points
    fd.write(str(len(points)) + "\n")

    # <lat> <long> <id#>
    for i, vert in enumerate( points):
        fd.write(str(round(vert[0],7)) + delimiter +
                 str(round(vert[1],7)) + delimiter + str(i) + "\n")

    fd.close()


##
# @brief ??
# @param point_atts ??
# @return ??
def _point_atts2array(point_atts):
    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)

    return point_atts


##
# @brief Convert geospatial object to a points dictionary.
# @param geospatial_data The geospatial object to convert.
# @return A points dictionary.
def geospatial_data2points_dictionary(geospatial_data):
    """Convert geospatial data to points_dictionary"""

    points_dictionary = {}
    points_dictionary['pointlist'] = geospatial_data.data_points

    points_dictionary['attributelist'] = {}

    for attribute_name in geospatial_data.attributes.keys():
        val = geospatial_data.attributes[attribute_name]
        points_dictionary['attributelist'][attribute_name] = val

    points_dictionary['geo_reference'] = geospatial_data.geo_reference

    return points_dictionary


##
# @brief Convert a points dictionary to a geospatial object.
# @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'"
    assert points_dictionary.has_key('pointlist'), msg

    msg = "Points dictionary must have key 'attributelist'"
    assert points_dictionary.has_key('attributelist'), msg

    if points_dictionary.has_key('geo_reference'):
        geo = points_dictionary['geo_reference']
    else:
        geo = None

    return Geospatial_data(points_dictionary['pointlist'],
                           points_dictionary['attributelist'],
                           geo_reference=geo)


##
# @brief Ensure that points are in absolute coordinates.
# @param points A list or array of points to check, or geospatial object.
# @param geo_reference If supplied, 
# @return ??
def ensure_absolute(points, geo_reference=None):
    """Ensure that points are in absolute coordinates.

    This function inputs several formats and
    outputs one format. - a numeric array of absolute points.

    Input formats are;
      points: List or numeric array of coordinate pairs [xi, eta] of
              points or geospatial object or points file name

    mesh_origin: A geo_reference object or 3-tuples consisting of
                 UTM zone, easting and northing.
                 If specified vertex coordinates are assumed to be
                 relative to their respective origins.
    """

    # Input check
    if isinstance(points, basestring):
        # It's a string - assume it is a point file
        points = Geospatial_data(file_name=points)

    if isinstance(points, Geospatial_data):
        points = points.get_data_points(absolute=True)
        msg = 'Use a Geospatial_data object or a mesh origin, not both.'
        assert geo_reference == None, msg
    else:
        points = ensure_numeric(points, num.float)

    # Sort of geo_reference and convert points
    if geo_reference is None:
        geo = None    # Geo_reference()
    else:
        if isinstance(geo_reference, Geo_reference):
            geo = geo_reference
        else:
            geo = Geo_reference(geo_reference[0],
                                geo_reference[1],
                                geo_reference[2])
        points = geo.get_absolute(points)

    return points


##
# @brief 
# @param points 
# @param geo_reference 
# @return A geospatial object.
def ensure_geospatial(points, geo_reference=None):
    """Convert various data formats to a geospatial_data instance.

    Inputed formats are;
    points:      List or numeric array of coordinate pairs [xi, eta] of
                 points or geospatial object

    mesh_origin: A geo_reference object or 3-tuples consisting of
                 UTM zone, easting and northing.
                 If specified vertex coordinates are assumed to be
                 relative to their respective origins.
    """

    # Input check
    if isinstance(points, Geospatial_data):
        msg = "Use a Geospatial_data object or a mesh origin, not both."
        assert geo_reference is None, msg
        return points
    else:
        # List or numeric array of absolute points
        points = ensure_numeric(points, num.float)

    # Sort out geo reference
    if geo_reference is None:
        geo = None
    else:
        if isinstance(geo_reference, Geo_reference):
            geo = geo_reference
        else:
            geo = Geo_reference(geo_reference[0],
                                geo_reference[1],
                                geo_reference[2])

    # Create Geospatial_data object with appropriate geo reference and return
    points = Geospatial_data(data_points=points, geo_reference=geo)

    return points


##
# @brief 
# @param data_file 
# @param alpha_list 
# @param mesh_file 
# @param boundary_poly 
# @param mesh_resolution 
# @param north_boundary 
# @param south_boundary 
# @param east_boundary 
# @param west_boundary 
# @param plot_name 
# @param split_factor 
# @param seed_num 
# @param cache 
# @param verbose 
def find_optimal_smoothing_parameter(data_file,
                                     alpha_list=None,
                                     mesh_file=None,
                                     boundary_poly=None,
                                     mesh_resolution=100000,
                                     north_boundary=None,
                                     south_boundary=None,
                                     east_boundary=None,
                                     west_boundary=None,
                                     plot_name='all_alphas',
                                     split_factor=0.1,
                                     seed_num=None,
                                     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.

    data_file: must not contain points outside the boundaries defined
               and it must be either a pts, txt or csv file.

    alpha_list: the alpha values to test in a single list

    mesh_file: name of the created mesh file or if passed in will read it.
               NOTE, if there is a mesh file mesh_resolution,
               north_boundary, south... etc will be ignored.

    mesh_resolution: the maximum area size for a triangle

    north_boundary... west_boundary: the value of the boundary

    plot_name: the name for the plot contain the results

    seed_num: the seed to the random number generator

    USAGE:
        value, alpha = find_optimal_smoothing_parameter(data_file=fileName,
                                             alpha_list=[0.0001, 0.01, 1],
                                             mesh_file=None,
                                             mesh_resolution=3,
                                             north_boundary=5,
                                             south_boundary=-5,
                                             east_boundary=5,
                                             west_boundary=-5,
                                             plot_name='all_alphas',
                                             seed_num=100000,
                                             verbose=False)

    OUTPUT: returns the minumum normalised covalance calculate AND the
           alpha that created it. PLUS writes a plot of the results

    NOTE: code will not work if the data_file extent is greater than the
    boundary_polygon or any of the boundaries, eg north_boundary...west_boundary
    """

    from anuga.shallow_water import Domain
    from anuga.geospatial_data.geospatial_data import Geospatial_data
    from anuga.pmesh.mesh_interface import create_mesh_from_regions
    from anuga.utilities.numerical_tools import cov
    from anuga.utilities.polygon import is_inside_polygon
    from anuga.fit_interpolate.benchmark_least_squares import mem_usage

    attribute_smoothed = 'elevation'

    if mesh_file is None:
        if verbose: print "building mesh"
        mesh_file = 'temp.msh'

        if (north_boundary is None or south_boundary is None
            or east_boundary is None or west_boundary is None):
            no_boundary = True
        else:
            no_boundary = False

        if no_boundary is True:
            msg = 'All boundaries must be defined'
            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,
                                 boundary_tags={'back': [2],
                                                'side': [1,3],
                                                'ocean': [0]},
                                 maximum_triangle_area=mesh_resolution,
                                 filename=mesh_file,
                                 use_cache=cache,
                                 verbose=verbose)

    else: # if mesh file provided
        # test mesh file exists?
        if verbose: "reading from file: %s" % mesh_file
        if access(mesh_file,F_OK) == 0:
            msg = "file %s doesn't exist!" % mesh_file
            raise IOError, msg

    # split topo data
    if verbose: print 'Reading elevation file: %s' % data_file
    G = Geospatial_data(file_name = data_file)
    if verbose: print 'Start split'
    G_small, G_other = G.split(split_factor, seed_num, verbose=verbose)
    if verbose: print 'Finish split'
    points = G_small.get_data_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]
    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([], 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]
    elevation_sample = G_small.get_attributes(attribute_name=attribute_smoothed)
    data[:,2] = elevation_sample

    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
        if verbose:
            print '\nCalculating domain and mesh for Alpha =', alpha, '\n'
        domain = Domain(mesh_file, use_cache=cache, verbose=verbose)
        if verbose: print domain.statistics()
        domain.set_quantity(attribute_smoothed,
                            geospatial_data=G_other,
                            use_cache=cache,
                            verbose=verbose,
                            alpha=alpha)

        # Convert points to geospatial data for use with get_values below
        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
        if verbose: print 'Get predicted elevation for location to be compared'
        elevation_predicted = \
                domain.quantities[attribute_smoothed].\
                    get_values(interpolation_points=points_geo)

        # 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]

        if verbose:
            print 'Covariance for alpha ', normal_cov[i][0], '= ', \
                      normal_cov[i][1]
            print '-------------------------------------------- \n'

    normal_cov0 = normal_cov[:,0]
    normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0), axis=0)

    if plot_name is not None:
        from pylab import savefig, semilogx, loglog

        semilogx(normal_cov_new[:,0], normal_cov_new[:,1])
        loglog(normal_cov_new[:,0], normal_cov_new[:,1])
        savefig(plot_name, dpi=300)

    if mesh_file == 'temp.msh':
        remove(mesh_file)

    if verbose:
        print 'Final results:'
        for i, alpha in enumerate(alphas):
            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
    return (min(normal_cov_new[:,1]),
            normal_cov_new[(num.argmin(normal_cov_new,axis=0))[1],0])


##
# @brief 
# @param data_file 
# @param alpha_list 
# @param mesh_file 
# @param boundary_poly 
# @param mesh_resolution 
# @param north_boundary 
# @param south_boundary 
# @param east_boundary 
# @param west_boundary 
# @param plot_name 
# @param split_factor 
# @param seed_num 
# @param cache 
# @param verbose 
def old_find_optimal_smoothing_parameter(data_file,
                                         alpha_list=None,
                                         mesh_file=None,
                                         boundary_poly=None,
                                         mesh_resolution=100000,
                                         north_boundary=None,
                                         south_boundary=None,
                                         east_boundary=None,
                                         west_boundary=None,
                                         plot_name='all_alphas',
                                         split_factor=0.1,
                                         seed_num=None,
                                         cache=False,
                                         verbose=False):
    """
    data_file: must not contain points outside the boundaries defined
               and it either a pts, txt or csv file.

    alpha_list: the alpha values to test in a single list

    mesh_file: name of the created mesh file or if passed in will read it.
               NOTE, if there is a mesh file mesh_resolution,
               north_boundary, south... etc will be ignored.

    mesh_resolution: the maximum area size for a triangle

    north_boundary... west_boundary: the value of the boundary

    plot_name: the name for the plot contain the results

    seed_num: the seed to the random number generator

    USAGE:
        value, alpha = find_optimal_smoothing_parameter(data_file=fileName,
                                             alpha_list=[0.0001, 0.01, 1],
                                             mesh_file=None,
                                             mesh_resolution=3,
                                             north_boundary=5,
                                             south_boundary=-5,
                                             east_boundary=5,
                                             west_boundary=-5,
                                             plot_name='all_alphas',
                                             seed_num=100000,
                                             verbose=False)

    OUTPUT: returns the minumum normalised covalance calculate AND the
            alpha that created it. PLUS writes a plot of the results

    NOTE: code will not work if the data_file extend is greater than the
          boundary_polygon or the north_boundary...west_boundary
    """

    from anuga.shallow_water import Domain
    from anuga.geospatial_data.geospatial_data import Geospatial_data
    from anuga.pmesh.mesh_interface import create_mesh_from_regions
    from anuga.utilities.numerical_tools import cov
    from anuga.utilities.polygon import is_inside_polygon
    from anuga.fit_interpolate.benchmark_least_squares import mem_usage

    attribute_smoothed = 'elevation'

    if mesh_file is None:
        mesh_file = 'temp.msh'

        if (north_boundary is None or south_boundary is None
            or east_boundary is None or west_boundary is None):
            no_boundary = True
        else:
            no_boundary = False

        if no_boundary is True:
            msg = 'All boundaries must be defined'
            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,
                                 boundary_tags={'back': [2],
                                                'side': [1,3],
                                                'ocean': [0]},
                                 maximum_triangle_area=mesh_resolution,
                                 filename=mesh_file,
                                 use_cache=cache,
                                 verbose=verbose)

    else: # if mesh file provided
        # 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
    G = Geospatial_data(file_name=data_file)
    if verbose: print 'start split'
    G_small, G_other = G.split(split_factor, seed_num, verbose=verbose)
    if verbose: print 'finish split'
    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]
    else:
        alphas = alpha_list

    domains = {}

    if verbose: print 'Setup computational domains with different alphas'

    for alpha in alphas:
        # add G_other data to domains with different alphas
        if verbose:
            print '\nCalculating domain and mesh for Alpha =', alpha, '\n'
        domain = Domain(mesh_file, use_cache=cache, verbose=verbose)
        if verbose: print domain.statistics()
        domain.set_quantity(attribute_smoothed,
                            geospatial_data=G_other,
                            use_cache=cache,
                            verbose=verbose,
                            alpha=alpha)
        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([], 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]
    elevation_sample = G_small.get_attributes(attribute_name=attribute_smoothed)
    data[:,2] = elevation_sample

    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
        elevation_predicted = \
                domains[alpha].quantities[attribute_smoothed].\
                        get_values(interpolation_points=points_geo)

        # 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]
        print 'memory usage during compare', mem_usage()
        if verbose: print 'cov', normal_cov[i][0], '= ', normal_cov[i][1]

    normal_cov0 = normal_cov[:,0]
    normal_cov_new = num.take(normal_cov, num.argsort(normal_cov0), axis=0)

    if plot_name is not None:
        from pylab import savefig,semilogx,loglog

        semilogx(normal_cov_new[:,0], normal_cov_new[:,1])
        loglog(normal_cov_new[:,0], normal_cov_new[:,1])
        savefig(plot_name, dpi=300)
    if mesh_file == 'temp.msh':
        remove(mesh_file)

    return (min(normal_cov_new[:,1]),
            normal_cov_new[(num.argmin(normal_cov_new, axis=0))[1],0])


if __name__ == "__main__":
    pass