Ignore:
Timestamp:
Feb 27, 2009, 11:54:09 AM (14 years ago)
Author:
rwilson
Message:

numpy changes

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/numpy/anuga/geospatial_data/geospatial_data.py

    r6410 r6428  
    1 """Class Geospatial_data - Manipulation of locations on the planet and
    2 associated attributes.
    3 
     1"""Class Geospatial_data
     2
     3Manipulation of locations on the planet and associated attributes.
    44"""
    55
    … …  
    1111from RandomArray import randint, seed, get_seed
    1212from copy import deepcopy
     13
    1314from Scientific.IO.NetCDF import NetCDFFile
    14 
    1515import numpy as num
    1616
    … …  
    2525from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    2626from anuga.config import netcdf_float
     27
    2728
    2829DEFAULT_ATTRIBUTE = 'elevation'
    … …  
    139140
    140141        self.set_verbose(verbose)
    141         self.geo_reference = None # create the attribute
     142        self.geo_reference = None
    142143        self.file_name = file_name
    143144
    … …  
    148149
    149150        if file_name is None:
    150             if latitudes is not None \
    151                or longitudes is not None \
    152                or points_are_lats_longs:
     151            if (latitudes is not None or longitudes is not None
     152                    or points_are_lats_longs):
    153153                data_points, geo_reference = \
    154154                    _set_using_lat_long(latitudes=latitudes,
    … …  
    156156                                        geo_reference=geo_reference,
    157157                                        data_points=data_points,
    158                                         points_are_lats_longs=\
    159                                         points_are_lats_longs)
     158                                        points_are_lats_longs=
     159                                            points_are_lats_longs)
    160160            self.check_data_points(data_points)
    161161            self.set_attributes(attributes)
    … …  
    213213            msg = 'There is no data or file provided!'
    214214            raise ValueError, msg
    215 
    216215        else:
    217216            self.data_points = ensure_numeric(data_points)
    … …  
    257256        """
    258257
    259         from anuga.coordinate_transforms.geo_reference import Geo_reference
    260 
    261258        if geo_reference is None:
    262259            # Use default - points are in absolute coordinates
    … …  
    269266            # FIXME (Ole): This exception will be raised even
    270267            # if geo_reference is None. Is that the intent Duncan?
    271             msg = 'Argument geo_reference must be a valid Geo_reference '
    272             msg += 'object or None.'
     268            msg = ('Argument geo_reference must be a valid Geo_reference '
     269                   'object or None.')
    273270            raise Expection, msg
    274271
    … …  
    378375    # @param isSouthHemisphere If True, return lat/lon points in S.Hemi.
    379376    # @return A set of data points, in appropriate form.
    380     def get_data_points(self, absolute=True, geo_reference=None,
    381                         as_lat_long=False, isSouthHemisphere=True):
     377    def get_data_points(self,
     378                        absolute=True,
     379                        geo_reference=None,
     380                        as_lat_long=False,
     381                        isSouthHemisphere=True):
    382382        """Get coordinates for all data points as an Nx2 array
    383383
    … …  
    460460    # @return The new geospatial object.
    461461    def __add__(self, other):
    462         """Returns the addition of 2 geospatical objects,
     462        """Returns the addition of 2 geospatial objects,
    463463        objects are concatencated to the end of each other
    464464
    … …  
    488488            if self.attributes is None:
    489489                if other.attributes is not None:
    490                     msg = 'Geospatial data must have the same '
    491                     msg += 'attributes to allow addition.'
     490                    msg = ('Geospatial data must have the same '
     491                           'attributes to allow addition.')
    492492                    raise Exception, msg
    493493
    … …  
    499499                        attrib1 = self.attributes[x]
    500500                        attrib2 = other.attributes[x]
    501                         new_attributes[x] = num.concatenate((attrib1, attrib2))
     501                        new_attributes[x] = num.concatenate((attrib1, attrib2),
     502                                                            axis=0) #??default#
    502503                    else:
    503                         msg = 'Geospatial data must have the same \n'
    504                         msg += 'attributes to allow addition.'
     504                        msg = ('Geospatial data must have the same '
     505                               'attributes to allow addition.')
    505506                        raise Exception, msg
    506507        else:
    … …  
    523524        return self + other
    524525
    525     ###
    526     #  IMPORT/EXPORT POINTS FILES
    527     ###
     526################################################################################
     527#  IMPORT/EXPORT POINTS FILES
     528################################################################################
    528529
    529530    ##
    … …  
    560561            except IOError, e:
    561562                # This should only be if a file is not found
    562                 msg = 'Could not open file %s. ' % file_name
    563                 msg += 'Check the file location.'
     563                msg = ('Could not open file %s. Check the file location.'
     564                       % file_name)
    564565                raise IOError, msg
    565566            except SyntaxError, e:
    566567                # This should only be if there is a format error
    567                 msg = 'Problem with format of file %s. \n' %file_name
    568                 msg += Error_message['IOError']
     568                msg = ('Problem with format of file %s.\n%s'
     569                       % (file_name, Error_message['IOError']))
    569570                raise SyntaxError, msg
    570571        else:
    571             msg = 'Extension %s is unknown' %file_name[-4:]
     572            msg = 'Extension %s is unknown' % file_name[-4:]
    572573            raise IOError, msg
    573574
    … …  
    614615                                self.get_all_attributes(),
    615616                                self.get_geo_reference())
    616 
    617617        elif file_name[-4:] == ".txt" or file_name[-4:] == ".csv":
    618618            msg = "ERROR: trying to write a .txt file with relative data."
    … …  
    621621                            self.get_data_points(absolute=True,
    622622                                                 as_lat_long=as_lat_long,
    623                                   isSouthHemisphere=isSouthHemisphere),
     623                                           isSouthHemisphere=isSouthHemisphere),
    624624                            self.get_all_attributes(),
    625625                            as_lat_long=as_lat_long)
    626 
    627626        elif file_name[-4:] == ".urs" :
    628627            msg = "ERROR: Can not write a .urs file as a relative file."
    … …  
    630629            _write_urs_file(file_name,
    631630                            self.get_data_points(as_lat_long=True,
    632                                   isSouthHemisphere=isSouthHemisphere))
    633 
     631                                           isSouthHemisphere=isSouthHemisphere))
    634632        else:
    635633            msg = 'Unknown file type %s ' %file_name
    … …  
    694692        random_list = []
    695693        remainder_list = []
    696         new_size = round(factor*self_size)
     694        new_size = round(factor * self_size)
    697695
    698696        # Find unique random numbers
    … …  
    713711        # Set seed if provided, mainly important for unit test!
    714712        # plus recalcule seed when no seed provided.
    715         if seed_num != None:
     713        if seed_num is not None:
    716714            seed(seed_num, seed_num)
    717715        else:
    … …  
    734732
    735733        # pops array index (random_num) from remainder_list
    736         # (which starts as the
    737         # total_list and appends to random_list
     734        # (which starts as the total_list and appends to random_list)
    738735        random_num_len = len(random_num)
    739736        for i in random_num:
    … …  
    746743
    747744        # FIXME: move to tests, it might take a long time
    748         # then create an array of random lenght between 500 and 1000,
     745        # then create an array of random length between 500 and 1000,
    749746        # and use a random factor between 0 and 1
    750747        # setup for assertion
    … …  
    752749        test_total.extend(remainder_list)
    753750        test_total.sort()
    754         msg = 'The two random lists made from the original list when added ' \
    755               'together DO NOT equal the original list'
     751        msg = ('The two random lists made from the original list when added '
     752               'together DO NOT equal the original list')
    756753        assert total_list == test_total, msg
    757754
    … …  
    770767        file pointer position
    771768        """
    772         from Scientific.IO.NetCDF import NetCDFFile
     769
    773770        # FIXME - what to do if the file isn't there
    774771
    … …  
    799796            self.number_of_blocks = self.number_of_points/self.max_read_lines
    800797            # This computes the number of full blocks. The last block may be
    801             # smaller and won't be ircluded in this estimate.
     798            # smaller and won't be included in this estimate.
    802799
    803800            if self.verbose is True:
    804                 print 'Reading %d points (in ~%d blocks) from file %s. ' \
    805                       % (self.number_of_points,
    806                          self.number_of_blocks,
    807                          self.file_name),
    808                 print 'Each block consists of %d data points' \
    809                       % self.max_read_lines
    810 
     801                print ('Reading %d points (in ~%d blocks) from file %s. '
     802                       % (self.number_of_points, self.number_of_blocks,
     803                          self.file_name)),
     804                print ('Each block consists of %d data points'
     805                       % self.max_read_lines)
    811806        else:
    812807            # Assume the file is a csv file
    … …  
    839834
    840835            if self.verbose is True:
    841                 if self.show_verbose >= self.start_row \
    842                    and self.show_verbose < fin_row:
    843                     print 'Reading block %d (points %d to %d) out of %d'\
    844                           %(self.block_number,
    845                             self.start_row,
    846                             fin_row,
    847                             self.number_of_blocks)
     836                if (self.show_verbose >= self.start_row
     837                    and self.show_verbose < fin_row):
     838                    print ('Reading block %d (points %d to %d) out of %d'
     839                           % (self.block_number, self.start_row,
     840                              fin_row, self.number_of_blocks))
    848841
    849842                    self.show_verbose += max(self.max_read_lines,
    850843                                             self.verbose_block_size)
    851 
    852844
    853845            # Read next block
    … …  
    861853
    862854            self.block_number += 1
    863 
    864855        else:
    865856            # Assume the file is a csv file
    … …  
    868859                 att_dict,
    869860                 geo_ref,
    870                  self.file_pointer) = \
    871                         _read_csv_file_blocking(self.file_pointer,
    872                                                 self.header[:],
    873                                                 max_read_lines=\
    874                                                     self.max_read_lines,
    875                                                 verbose=self.verbose)
     861                 self.file_pointer) = _read_csv_file_blocking(self.file_pointer,
     862                                                              self.header[:],
     863                                                              max_read_lines= \
     864                                                           self.max_read_lines,
     865                                                              verbose= \
     866                                                                  self.verbose)
    876867
    877868                # Check that the zones haven't changed.
    … …  
    880871                    self.blocking_georef = geo_ref
    881872                elif self.blocking_georef is not None:
    882                     msg = 'Geo reference given, then not given.'
    883                     msg += ' This should not happen.'
     873                    msg = ('Geo reference given, then not given.'
     874                           ' This should not happen.')
    884875                    raise ValueError, msg
    885876                geo = Geospatial_data(pointlist, att_dict, geo_ref)
    … …  
    899890                del self.file_pointer
    900891                # This should only be if there is a format error
    901                 msg = 'Could not open file %s. \n' % self.file_name
    902                 msg += Error_message['IOError']
     892                msg = ('Could not open file %s.\n%s'
     893                       % (self.file_name, Error_message['IOError']))
    903894                raise SyntaxError, msg
    904895        return geo
    905 
    906896
    907897##################### Error messages ###########
    908898Error_message = {}
    909899Em = Error_message
    910 Em['IOError'] = "NOTE: The format for a comma separated .txt/.csv file is:\n"
    911 Em['IOError'] += "        1st line:     [column names]\n"
    912 Em['IOError'] += "        other lines:  [x value], [y value], [attributes]\n"
    913 Em['IOError'] += "\n"
    914 Em['IOError'] += "           for example:\n"
    915 Em['IOError'] += "           x, y, elevation, friction\n"
    916 Em['IOError'] += "           0.6, 0.7, 4.9, 0.3\n"
    917 Em['IOError'] += "           1.9, 2.8, 5, 0.3\n"
    918 Em['IOError'] += "           2.7, 2.4, 5.2, 0.3\n"
    919 Em['IOError'] += "\n"
    920 Em['IOError'] += "The first two columns are assumed to be x, y coordinates.\n"
    921 Em['IOError'] += "The attribute values must be numeric.\n"
     900Em['IOError'] = ('NOTE: The format for a comma separated .txt/.csv file is:\n'
     901                 '        1st line:     [column names]\n'
     902                 '        other lines:  [x value], [y value], [attributes]\n'
     903                 '\n'
     904                 '           for example:\n'
     905                 '           x, y, elevation, friction\n'
     906                 '           0.6, 0.7, 4.9, 0.3\n'
     907                 '           1.9, 2.8, 5, 0.3\n'
     908                 '           2.7, 2.4, 5.2, 0.3\n'
     909                 '\n'
     910                 'The first two columns are assumed to be x, y coordinates.\n'
     911                 'The attribute values must be numeric.\n')
    922912
    923913##
    … …  
    936926
    937927    if geo_reference is not None:
    938         msg = "A georeference is specified yet latitude and longitude " \
    939               "are also specified!"
     928        msg = ('A georeference is specified yet latitude and longitude '
     929               'are also specified!')
    940930        raise ValueError, msg
    941931
    942932    if data_points is not None and not points_are_lats_longs:
    943         msg = "Data points are specified yet latitude and longitude are " \
    944               "also specified."
     933        msg = ('Data points are specified yet latitude and longitude are '
     934               'also specified.')
    945935        raise ValueError, msg
    946936
    … …  
    983973    dic['attributelist']['elevation'] = [[7.0,5.0]]
    984974    """
    985 
    986     from Scientific.IO.NetCDF import NetCDFFile
    987975
    988976    if verbose: print 'Reading ', file_name
    … …  
    11271115            if len(header) != len(numbers):
    11281116                file_pointer.close()
    1129                 msg = "File load error. " \
    1130                       "There might be a problem with the file header."
     1117                msg = ('File load error. '
     1118                       'There might be a problem with the file header.')
    11311119                raise SyntaxError, msg
    11321120            for i,n in enumerate(numbers):
    11331121                n.strip()
    11341122                if n != '\n' and n != '':
    1135                     #attributes.append(float(n))
    11361123                    att_dict.setdefault(header[i],[]).append(float(n))
    11371124        except ValueError:
    … …  
    11741161# @note Will throw IOError and AttributeError exceptions.
    11751162def _read_pts_file_header(fid, verbose=False):
    1176     """Read the geo_reference and number_of_points from a .pts file"""
     1163    '''Read the geo_reference and number_of_points from a .pts file'''
    11771164
    11781165    keys = fid.variables.keys()
    … …  
    12021189# @return Tuple of (pointlist, attributes).
    12031190def _read_pts_file_blocking(fid, start_row, fin_row, keys):
    1204     """Read the body of a .csv file."""
     1191    '''Read the body of a .csv file.'''
    12051192
    12061193    pointlist = num.array(fid.variables['points'][start_row:fin_row])
    … …  
    12401227    """
    12411228
    1242     from Scientific.IO.NetCDF import NetCDFFile
    1243 
    12441229    # NetCDF file definition
    12451230    outfile = NetCDFFile(file_name, netcdf_mode_w)
    … …  
    12561241
    12571242    # Variable definition
    1258     outfile.createVariable('points', netcdf_float, ('number_of_points',
    1259                                                    'number_of_dimensions'))
     1243    outfile.createVariable('points', netcdf_float,
     1244                           ('number_of_points', 'number_of_dimensions'))
    12601245
    12611246    # create variables
    … …  
    13851370    """Convert points_dictionary to geospatial data object"""
    13861371
    1387     msg = 'Points dictionary must have key pointlist'
     1372    msg = "Points dictionary must have key 'pointlist'"
    13881373    assert points_dictionary.has_key('pointlist'), msg
    13891374
    1390     msg = 'Points dictionary must have key attributelist'
     1375    msg = "Points dictionary must have key 'attributelist'"
    13911376    assert points_dictionary.has_key('attributelist'), msg
    13921377
    … …  
    13981383    return Geospatial_data(points_dictionary['pointlist'],
    13991384                           points_dictionary['attributelist'],
    1400                            geo_reference = geo)
     1385                           geo_reference=geo)
    14011386
    14021387
    … …  
    14361421    # Sort of geo_reference and convert points
    14371422    if geo_reference is None:
    1438         geo = None # Geo_reference()
     1423        geo = None    # Geo_reference()
    14391424    else:
    14401425        if isinstance(geo_reference, Geo_reference):
    … …  
    15721557    from anuga.fit_interpolate.benchmark_least_squares import mem_usage
    15731558
    1574 
    15751559    attribute_smoothed = 'elevation'
    15761560
    … …  
    15791563        mesh_file = 'temp.msh'
    15801564
    1581         if north_boundary is None or south_boundary is None \
    1582            or east_boundary is None or west_boundary is None:
     1565        if (north_boundary is None or south_boundary is None
     1566            or east_boundary is None or west_boundary is None):
    15831567            no_boundary = True
    15841568        else:
    … …  
    15891573            raise Expection, msg
    15901574
    1591         poly_topo = [[east_boundary,south_boundary],
    1592                      [east_boundary,north_boundary],
    1593                      [west_boundary,north_boundary],
    1594                      [west_boundary,south_boundary]]
     1575        poly_topo = [[east_boundary, south_boundary],
     1576                     [east_boundary, north_boundary],
     1577                     [west_boundary, north_boundary],
     1578                     [west_boundary, south_boundary]]
    15951579
    15961580        create_mesh_from_regions(poly_topo,
    … …  
    16221606    if alpha_list == None:
    16231607        alphas = [0.001,0.01,100]
    1624         #alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01,\
    1625          #         0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
     1608        #alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01,
     1609        #          0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
    16261610    else:
    16271611        alphas = alpha_list
    … …  
    16471631        # add G_other data to domains with different alphas
    16481632        if verbose:
    1649             print '\n Calculating domain and mesh for Alpha = ', alpha, '\n'
     1633            print '\nCalculating domain and mesh for Alpha =', alpha, '\n'
    16501634        domain = Domain(mesh_file, use_cache=cache, verbose=verbose)
    16511635        if verbose: print domain.statistics()
    … …  
    16711655        sample_cov = cov(elevation_sample)
    16721656        ele_cov = cov(elevation_sample - elevation_predicted)
    1673         normal_cov[i,:] = [alpha,ele_cov / sample_cov]
     1657        normal_cov[i,:] = [alpha, ele_cov / sample_cov]
    16741658
    16751659        if verbose:
    … …  
    16941678        print 'Final results:'
    16951679        for i, alpha in enumerate(alphas):
    1696             print 'covariance for alpha %s = %s ' \
    1697                   % (normal_cov[i][0], normal_cov[i][1])
    1698         print '\n Optimal alpha is: %s ' \
    1699               % normal_cov_new[(num.argmin(normal_cov_new, axis=0))[1], 0]
     1680            print ('covariance for alpha %s = %s '
     1681                   % (normal_cov[i][0], normal_cov[i][1]))
     1682        print ('\nOptimal alpha is: %s '
     1683               % normal_cov_new[(num.argmin(normal_cov_new, axis=0))[1], 0])
    17001684
    17011685    # covariance and optimal alpha
    … …  
    17791763    from anuga.fit_interpolate.benchmark_least_squares import mem_usage
    17801764
    1781 
    17821765    attribute_smoothed = 'elevation'
    17831766
    … …  
    17851768        mesh_file = 'temp.msh'
    17861769
    1787         if north_boundary is None or south_boundary is None \
    1788            or east_boundary is None or west_boundary is None:
     1770        if (north_boundary is None or south_boundary is None
     1771            or east_boundary is None or west_boundary is None):
    17891772            no_boundary = True
    17901773        else:
    … …  
    17951778            raise Expection, msg
    17961779
    1797         poly_topo = [[east_boundary,south_boundary],
    1798                      [east_boundary,north_boundary],
    1799                      [west_boundary,north_boundary],
    1800                      [west_boundary,south_boundary]]
     1780        poly_topo = [[east_boundary, south_boundary],
     1781                     [east_boundary, north_boundary],
     1782                     [west_boundary, north_boundary],
     1783                     [west_boundary, south_boundary]]
    18011784
    18021785        create_mesh_from_regions(poly_topo,
    … …  
    18221805    points = G_small.get_data_points()
    18231806
    1824     # FIXME: Remove points outside boundary polygon
    1825 #    print 'new point',len(points)
    1826 #
    1827 #    new_points=[]
    1828 #    new_points=array([],dtype=float)
    1829 #    new_points=resize(new_points,(len(points),2))
    1830 #    print "BOUNDARY", boundary_poly
    1831 #    for i,point in enumerate(points):
    1832 #        if is_inside_polygon(point,boundary_poly, verbose=True):
    1833 #            new_points[i] = point
    1834 #            print"WOW",i,new_points[i]
    1835 #    points = new_points
    1836 
    18371807    if verbose: print "Number of points in sample to compare: ", len(points)
    18381808
    18391809    if alpha_list == None:
    18401810        alphas = [0.001,0.01,100]
    1841         #alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01,\
    1842          #         0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
     1811        #alphas = [0.000001, 0.00001, 0.0001, 0.001, 0.01,
     1812        #          0.1, 1.0, 10.0, 100.0,1000.0,10000.0]
    18431813    else:
    18441814        alphas = alpha_list
    … …  
    18511821        # add G_other data to domains with different alphas
    18521822        if verbose:
    1853             print '\n Calculating domain and mesh for Alpha = ', alpha, '\n'
     1823            print '\nCalculating domain and mesh for Alpha =', alpha, '\n'
    18541824        domain = Domain(mesh_file, use_cache=cache, verbose=verbose)
    18551825        if verbose: print domain.statistics()
    … …  
    18781848    if verbose:
    18791849        print 'Determine difference between predicted results and actual data'
     1850
    18801851    for i, alpha in enumerate(domains):
    18811852        if verbose: print'Alpha =', alpha
Note: See TracChangeset for help on using the changeset viewer.