Context Navigation

← Previous Change
Next Change →

geospatial_data.py

Timestamp:

Feb 10, 2009, 11:11:04 AM (15 years ago)

Author:

rwilson

Message:

Initial commit of numpy changes. Still a long way to go.

Location:

branches/numpy

Files:

: 1 edited
: 1 copied

. (copied) (copied from anuga_core/source)
anuga/geospatial_data/geospatial_data.py (modified) (55 diffs)

Legend:

: Unmodified
: Added
: Removed

branches/numpy/anuga/geospatial_data/geospatial_data.py

-                      r6166
+                      r6304
 from Scientific.IO.NetCDF import NetCDFFile
 import Numeric as num
+import numpy as num
 from anuga.coordinate_transforms.lat_long_UTM_conversion import UTMtoLL
 …
 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'
 …
                  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.
 …
         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 # create the attribute
         self.file_name = file_name
 …
                                         data_points=data_points,
                                         points_are_lats_longs=\
                                             points_are_lats_longs)
+                                        points_are_lats_longs)
             self.check_data_points(data_points)
             self.set_attributes(attributes)
 …
                         # 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'
     ##
 …
     ##
     # @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:
 …
         else:
             self.data_points = ensure_numeric(data_points)
+            return
+            print 'self.data_points=%s' % str(self.data_points)
+            print 'self.data_points.shape=%s' % str(self.data_points.shape)
             if not (0,) == self.data_points.shape:
                 assert len(self.data_points.shape) == 2
 …
     # @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:
 …
         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):
 …
             msg = 'Argument geo_reference must be a valid Geo_reference '
             msg += 'object or None.'
             raise msg
+            raise Expection, msg
         # If a geo_reference already exists, change the point data according to
 …
                         raise Exception, msg
         else:
             #other is None:
+            # other is None:
             new_points = self.get_data_points(absolute=True)
             new_attributes = self.attributes
 …
     # @return The new geospatial object.
     def __radd__(self, other):
+        """Handle cases like None + Geospatial_data(...)
+        """
+        """Handle cases like None + Geospatial_data(...)"""
         return self + other
 …
     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!
 …
             except SyntaxError, e:
                 # This should only be if there is a format error
                 msg = 'Could not open file %s. \n' %file_name
+                msg = 'Problem with format of file %s. \n' %file_name
                 msg += Error_message['IOError']
                 raise SyntaxError, msg
 …
                            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.
 …
         """
         #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)
 …
     # @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.
 …
         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))
         total_list = total.tolist()
 …
         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()
 …
             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)
 …
         """
         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?
 …
                         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:
 …
                                                  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:
 …
     """
     line = file_pointer.readline()
+    line = file_pointer.readline().strip()
     header = clean_line(line, delimiter)
 …
 # @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,
 …
         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
 …
 # @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()
 …
 # @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 .csv file."""
     pointlist = num.array(fid.variables['points'][start_row:fin_row])
 …
     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
 …
     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:
 …
                     as_lat_long=False,
                     delimiter=','):
+    """Write a .csv file.
+    """
+    """Write a .csv file."""
     points = write_data_points
 …
 # @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
 …
 # @return A points dictionary.
 def geospatial_data2points_dictionary(geospatial_data):
+    """Convert geospatial data to points_dictionary
+    """
+    """Convert geospatial data to points_dictionary"""
     points_dictionary = {}
 …
 # @param points_dictionary A points dictionary to convert.
 def points_dictionary2geospatial_data(points_dictionary):
+    """Convert points_dictionary to geospatial data object
+    """
+    """Convert points_dictionary to geospatial data object"""
     msg = 'Points dictionary must have key pointlist'
 …
 ##
 # @brief Split a string into 'clean' fields.
 # @param line The string to process.
+# @param str 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
+# @note If a field contains '\n' it isn't zero length.
+def clean_line(str, delimiter):
+    """Split string on given delimiter, remove whitespace from each field."""
+    return [x.strip() for x in str.split(delimiter) if x != '']
 …
     # 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)
 …
         assert geo_reference == None, msg
     else:
         points = ensure_numeric(points, num.Float)
+        points = ensure_numeric(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):
 …
 # @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;
 …
     else:
         # List or numeric array of absolute points
         points = ensure_numeric(points, num.Float)
+        points = ensure_numeric(points, num.float)
     # Sort out geo reference
 …
                                      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
 …
         if no_boundary is True:
             msg = 'All boundaries must be defined'
             raise msg
+            raise Expection, msg
         poly_topo = [[east_boundary,south_boundary],
 …
     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:
 …
             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)
 …
         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)
+    # 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]
 …
     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'
 …
         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 = \
 …
                     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
 …
         if no_boundary is True:
             msg = 'All boundaries must be defined'
             raise msg
+            raise Expection, msg
         poly_topo = [[east_boundary,south_boundary],
 …
     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'
 …
     points = G_small.get_data_points()
     #FIXME: Remove points outside boundary polygon
+    # FIXME: Remove points outside boundary polygon
 #    print 'new point',len(points)
+#
 #    new_points=[]
 #    new_points=array([],typecode=Float)
+#    new_points=array([],dtype=float)
 #    new_points=resize(new_points,(len(points),2))
 #    print "BOUNDARY", boundary_poly
 …
     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'
 …
         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]
 …
     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:
 …
         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

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 6304 for branches/numpy/anuga/geospatial_data/geospatial_data.py

Legend:

branches/numpy/anuga/geospatial_data/geospatial_data.py

Download in other formats: