Changeset 6304


Ignore:
Timestamp:
Feb 10, 2009, 11:11:04 AM (16 years ago)
Author:
rwilson
Message:

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

Location:
branches/numpy
Files:
1 deleted
93 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/numpy/anuga/abstract_2d_finite_volumes/domain.py

    r6246 r6304  
    3131from quantity import Quantity
    3232
    33 import Numeric as num
     33import numpy as num
    3434
    3535
    … …  
    185185            buffer_shape = self.full_send_dict[key][0].shape[0]
    186186            self.full_send_dict[key].append(num.zeros((buffer_shape, self.nsys),
    187                                                       num.Float))
     187                                                      num.float))
    188188
    189189        for key in self.ghost_recv_dict:
    190190            buffer_shape = self.ghost_recv_dict[key][0].shape[0]
    191191            self.ghost_recv_dict[key].append(num.zeros((buffer_shape, self.nsys),
    192                                              num.Float))
     192                                             num.float))
    193193
    194194        # Setup cell full flag
    … …  
    197197        N = len(self) #number_of_elements
    198198        self.number_of_elements = N
    199         self.tri_full_flag = num.ones(N, num.Int)
     199        self.tri_full_flag = num.ones(N, num.int)
    200200        for i in self.ghost_recv_dict.keys():
    201201            for id in self.ghost_recv_dict[i][0]:
    … …  
    258258        # (boolean) array, to be used during the flux calculation.
    259259        N = len(self) # Number_of_triangles
    260         self.already_computed_flux = num.zeros((N, 3), num.Int)
     260        self.already_computed_flux = num.zeros((N, 3), num.int)
    261261
    262262        # Storage for maximal speeds computed for each triangle by
    263263        # compute_fluxes.
    264264        # This is used for diagnostics only (reset at every yieldstep)
    265         self.max_speed = num.zeros(N, num.Float)
     265        self.max_speed = num.zeros(N, num.float)
    266266
    267267        if mesh_filename is not None:
    … …  
    375375            raise Exception, msg
    376376
    377         q = num.zeros(len(self.conserved_quantities), num.Float)
     377        q = num.zeros(len(self.conserved_quantities), num.float)
    378378
    379379        for i, name in enumerate(self.conserved_quantities):
    … …  
    447447
    448448        name:  Name of quantity
    449         value: Compatible list, Numeric array, const or function (see below)
     449        value: Compatible list, numeric array, const or function (see below)
    450450
    451451        The values will be stored in elements following their internal ordering.
    … …  
    879879
    880880            msg += '------------------------------------------------\n'
    881             msg += '  Speeds in [%f, %f]\n' % (min(self.max_speed),
    882                                                max(self.max_speed))
     881            msg += '  Speeds in [%f, %f]\n' % (min(self.max_speed.flat),
     882                                               max(self.max_speed.flat))
    883883            msg += '  Histogram:\n'
    884884
    … …  
    892892                else:
    893893                    # Closed upper interval
    894                     hi = max(self.max_speed)
     894                    hi = max(self.max_speed.flat)
    895895                    msg += '    [%f, %f]: %d\n' % (lo, hi, count)
    896896
    897             N = len(self.max_speed)
     897            N = len(self.max_speed.flat)
    898898            if N > 10:
    899899                msg += '  Percentiles (10%):\n'
    … …  
    13571357                self.number_of_steps = 0
    13581358                self.number_of_first_order_steps = 0
    1359                 self.max_speed = num.zeros(N, num.Float)
     1359                self.max_speed = num.zeros(N, num.float)
    13601360
    13611361    ##
    … …  
    16261626        # momentum
    16271627        if self.protect_against_isolated_degenerate_timesteps is True and \
    1628                self.max_speed > 10.0: # FIXME (Ole): Make this configurable
     1628               num.max(self.max_speed) > 10.0: # FIXME (Ole): Make this configurable
    16291629
    16301630            # Setup 10 bins for speed histogram
    … …  
    16411641
    16421642                    # Find triangles in last bin
    1643                     # FIXME - speed up using Numeric
     1643                    # FIXME - speed up using numeric package
    16441644                    d = 0
    16451645                    for i in range(self.number_of_full_triangles):

  • branches/numpy/anuga/abstract_2d_finite_volumes/ermapper_grids.py

    r6161 r6304  
    22
    33# from os import open, write, read
    4 import Numeric as num
    5 
    6 celltype_map = {'IEEE4ByteReal': num.Float32, 'IEEE8ByteReal': num.Float64}
     4import numpy as num
     5
     6celltype_map = {'IEEE4ByteReal': num.float32, 'IEEE8ByteReal': num.float64}
    77
    88
    … …  
    1111    write_ermapper_grid(ofile, data, header = {}):
    1212   
    13     Function to write a 2D Numeric array to an ERMapper grid.  There are a series of conventions adopted within
     13    Function to write a 2D numeric array to an ERMapper grid.  There are a series of conventions adopted within
    1414    this code, specifically:
    1515    1)  The registration coordinate for the data is the SW (or lower-left) corner of the data
    1616    2)  The registration coordinates refer to cell centres
    17     3)  The data is a 2D Numeric array with the NW-most data in element (0,0) and the SE-most data in element (N,M)
     17    3)  The data is a 2D numeric array with the NW-most data in element (0,0) and the SE-most data in element (N,M)
    1818        where N is the last line and M is the last column
    1919    4)  There has been no testng of the use of a rotated grid.  Best to keep data in an NS orientation
    … …  
    163163    return header                     
    164164
    165 def write_ermapper_data(grid, ofile, data_format = num.Float32):
     165def write_ermapper_data(grid, ofile, data_format=num.float32):
    166166
    167167
    … …  
    193193
    194194
    195 def read_ermapper_data(ifile, data_format = num.Float32):
     195def read_ermapper_data(ifile, data_format = num.float32):
    196196    # open input file in a binary format and read the input string
    197197    fid = open(ifile,'rb')
    … …  
    199199    fid.close()
    200200
    201     # convert input string to required format (Note default format is num.Float32)
     201    # convert input string to required format (Note default format is num.float32)
    202202    grid_as_float = num.fromstring(input_string,data_format)
    203203    return grid_as_float

  • branches/numpy/anuga/abstract_2d_finite_volumes/general_mesh.py

    r6191 r6304  
    1 import Numeric as num
     1import numpy as num
    22
    33from anuga.coordinate_transforms.geo_reference import Geo_reference
    … …  
    1010    which is represented as a coordinate set (x,y).
    1111    Vertices from different triangles can point to the same node.
    12     The nodes are implemented as an Nx2 Numeric array containing the
     12    The nodes are implemented as an Nx2 numeric array containing the
    1313    x and y coordinates.
    1414
    … …  
    1919    where
    2020
    21       nodes is either a list of 2-tuples or an Nx2 Numeric array of
     21      nodes is either a list of 2-tuples or an Nx2 numeric array of
    2222      floats representing all x, y coordinates in the mesh.
    2323
    24       triangles is either a list of 3-tuples or an Mx3 Numeric array of
     24      triangles is either a list of 3-tuples or an Mx3 numeric array of
    2525      integers representing indices of all vertices in the mesh.
    2626      Each vertex is identified by its index i in [0, N-1].
    … …  
    6868       
    6969          nodes: x,y coordinates represented as a sequence of 2-tuples or
    70                  a Nx2 Numeric array of floats.
     70                 a Nx2 numeric array of floats.
    7171                 
    72           triangles: sequence of 3-tuples or Mx3 Numeric array of
     72          triangles: sequence of 3-tuples or Mx3 numeric array of
    7373                     non-negative integers representing indices into
    7474                     the nodes array.
    … …  
    8888        if verbose: print 'General_mesh: Building basic mesh structure in ANUGA domain'
    8989
    90         self.triangles = num.array(triangles, num.Int)
    91         self.nodes = num.array(nodes, num.Float)
     90        self.triangles = num.array(triangles, num.int)
     91        self.nodes = num.array(nodes, num.float)
    9292
    9393
    … …  
    125125
    126126        # Input checks
    127         msg = 'Triangles must an Mx3 Numeric array or a sequence of 3-tuples. '
     127        msg = 'Triangles must an Mx3 numeric array or a sequence of 3-tuples. '
    128128        msg += 'The supplied array has the shape: %s'\
    129129               %str(self.triangles.shape)
    130130        assert len(self.triangles.shape) == 2, msg
    131131
    132         msg = 'Nodes must an Nx2 Numeric array or a sequence of 2-tuples'
     132        msg = 'Nodes must an Nx2 numeric array or a sequence of 2-tuples'
    133133        msg += 'The supplied array has the shape: %s'\
    134134               %str(self.nodes.shape)
    … …  
    144144                      max(self.nodes[:,0]), max(self.nodes[:,1]) ]
    145145
    146         self.xy_extent = num.array(xy_extent, num.Float)
     146        self.xy_extent = num.array(xy_extent, num.float)
    147147
    148148
    149149        # Allocate space for geometric quantities
    150         self.normals = num.zeros((N, 6), num.Float)
    151         self.areas = num.zeros(N, num.Float)
    152         self.edgelengths = num.zeros((N, 3), num.Float)
     150        self.normals = num.zeros((N, 6), num.float)
     151        self.areas = num.zeros(N, num.float)
     152        self.edgelengths = num.zeros((N, 3), num.float)
    153153
    154154        # Get x,y coordinates for all triangles and store
    … …  
    185185            #   - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle
    186186
    187             n0 = num.array([x2 - x1, y2 - y1], num.Float)
     187            n0 = num.array([x2 - x1, y2 - y1], num.float)
    188188            l0 = num.sqrt(num.sum(n0**2))
    189189
    190             n1 = num.array([x0 - x2, y0 - y2], num.Float)
     190            n1 = num.array([x0 - x2, y0 - y2], num.float)
    191191            l1 = num.sqrt(num.sum(n1**2))
    192192
    193             n2 = num.array([x1 - x0, y1 - y0], num.Float)
     193            n2 = num.array([x1 - x0, y1 - y0], num.float)
    194194            l2 = num.sqrt(num.sum(n2**2))
    195195
    … …  
    275275            if not self.geo_reference.is_absolute():
    276276                return V + num.array([self.geo_reference.get_xllcorner(),
    277                                       self.geo_reference.get_yllcorner()], num.Float)
     277                                      self.geo_reference.get_yllcorner()], num.float)
    278278            else:
    279279                return V
    … …  
    317317            if absolute is True and not self.geo_reference.is_absolute():
    318318                offset=num.array([self.geo_reference.get_xllcorner(),
    319                                   self.geo_reference.get_yllcorner()], num.Float)
     319                                  self.geo_reference.get_yllcorner()], num.float)
    320320                return num.array([V[i3,:]+offset,
    321321                                  V[i3+1,:]+offset,
    322                                   V[i3+2,:]+offset], num.Float)
     322                                  V[i3+2,:]+offset], num.float)
    323323            else:
    324                 return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.Float)
     324                return num.array([V[i3,:], V[i3+1,:], V[i3+2,:]], num.float)
    325325               
    326326
    … …  
    349349
    350350        M = self.number_of_triangles
    351         vertex_coordinates = num.zeros((3*M, 2), num.Float)
     351        vertex_coordinates = num.zeros((3*M, 2), num.float)
    352352
    353353        for i in range(M):
    … …  
    409409        K = 3*M       # Total number of unique vertices
    410410       
    411         T = num.reshape(num.arange(K, typecode=num.Int), (M,3))
     411        T = num.reshape(num.arange(K, dtype=num.int), (M,3))
    412412       
    413413        return T     
    … …  
    416416
    417417    def get_unique_vertices(self,  indices=None):
    418         """FIXME(Ole): This function needs a docstring
    419         """
     418        """FIXME(Ole): This function needs a docstring"""
     419
    420420        triangles = self.get_triangles(indices=indices)
    421421        unique_verts = {}
    422422        for triangle in triangles:
    423            unique_verts[triangle[0]] = 0
    424            unique_verts[triangle[1]] = 0
    425            unique_verts[triangle[2]] = 0
     423            #print 'triangle(%s)=%s' % (type(triangle), str(triangle))
     424            unique_verts[triangle[0]] = 0
     425            unique_verts[triangle[1]] = 0
     426            unique_verts[triangle[2]] = 0
    426427        return unique_verts.keys()
    427428
    … …  
    452453                triangle_list.append( (volume_id, vertex_id) )
    453454
    454             triangle_list = num.array(triangle_list, num.Int)    #array default#
     455            triangle_list = num.array(triangle_list, num.int)    #array default#
    455456        else:
    456457            # Get info for all nodes recursively.
    … …  
    529530
    530531        # Count number of triangles per node
    531         number_of_triangles_per_node = num.zeros(self.number_of_full_nodes)
     532        number_of_triangles_per_node = num.zeros(self.number_of_full_nodes,
     533                                                 num.int)       #array default#
    532534        for volume_id, triangle in enumerate(self.get_triangles()):
    533535            for vertex_id in triangle:
    … …  
    536538        # Allocate space for inverted structure
    537539        number_of_entries = num.sum(number_of_triangles_per_node)
    538         vertex_value_indices = num.zeros(number_of_entries)
     540        vertex_value_indices = num.zeros(number_of_entries, num.int) #array default#
    539541
    540542        # Register (triangle, vertex) indices for each node

  • branches/numpy/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py

    r6174 r6304  
    88from anuga.fit_interpolate.interpolate import Modeltime_too_early
    99
    10 import Numeric as num
     10import numpy as num
    1111
    1212
    … …  
    6969            raise Exception, msg
    7070
    71         self.conserved_quantities=num.array(conserved_quantities, num.Float)
     71        self.conserved_quantities=num.array(conserved_quantities, num.float)
    7272
    7373    def __repr__(self):
    … …  
    9999
    100100        try:
    101             q = num.array(q, num.Float)
     101            q = num.array(q, num.float)
    102102        except:
    103103            msg = 'Return value from time boundary function could '
    104             msg += 'not be converted into a Numeric array of floats.\n'
     104            msg += 'not be converted into a numeric array of floats.\n'
    105105            msg += 'Specified function should return either list or array.\n'
    106106            msg += 'I got %s' %str(q)
    … …  
    180180
    181181        if verbose: print 'Find midpoint coordinates of entire boundary'
    182         self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.Float)
     182        self.midpoint_coordinates = num.zeros((len(domain.boundary), 2), num.float)
    183183        boundary_keys = domain.boundary.keys()
    184184
    … …  
    200200           
    201201            # Compute midpoints
    202             if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2], num.Float)
    203             if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2], num.Float)
    204             if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2], num.Float)
     202            if edge_id == 0: m = num.array([(x1 + x2)/2, (y1 + y2)/2], num.float)
     203            if edge_id == 1: m = num.array([(x0 + x2)/2, (y0 + y2)/2], num.float)
     204            if edge_id == 2: m = num.array([(x1 + x0)/2, (y1 + y0)/2], num.float)
    205205
    206206            # Convert to absolute UTM coordinates
    … …  
    309309                        # See http://docs.python.org/lib/warning-filter.html
    310310                        self.default_boundary_invoked = True
    311                    
    312 
    313             if res == NAN:
     311           
     312            if num.any(res == NAN):
    314313                x,y=self.midpoint_coordinates[i,:]
    315314                msg = 'NAN value found in file_boundary at '

  • branches/numpy/anuga/abstract_2d_finite_volumes/mesh_factory.py

    r6145 r6304  
    33"""
    44
    5 import Numeric as num
     5import numpy as num
    66
    77
    … …  
    9494    index = Index(n,m)
    9595
    96     points = num.zeros((Np, 2), num.Float)
     96    points = num.zeros((Np, 2), num.float)
    9797
    9898    for i in range(m+1):
    … …  
    106106
    107107
    108     elements = num.zeros((Nt, 3), num.Int)
     108    elements = num.zeros((Nt, 3), num.int)
    109109    boundary = {}
    110110    nt = -1

  • branches/numpy/anuga/abstract_2d_finite_volumes/neighbour_mesh.py

    r6191 r6304  
    1212from math import pi, sqrt
    1313
    14 import Numeric as num
     14import numpy as num
    1515       
    1616
    … …  
    2424    coordinate set.
    2525
    26     Coordinate sets are implemented as an N x 2 Numeric array containing
     26    Coordinate sets are implemented as an N x 2 numeric array containing
    2727    x and y coordinates.
    2828
    … …  
    3333    where
    3434
    35       coordinates is either a list of 2-tuples or an Mx2 Numeric array of
     35      coordinates is either a list of 2-tuples or an Mx2 numeric array of
    3636      floats representing all x, y coordinates in the mesh.
    3737
    38       triangles is either a list of 3-tuples or an Nx3 Numeric array of
     38      triangles is either a list of 3-tuples or an Nx3 numeric array of
    3939      integers representing indices of all vertices in the mesh.
    4040      Each vertex is identified by its index i in [0, M-1].
    … …  
    7676        """
    7777        Build triangles from x,y coordinates (sequence of 2-tuples or
    78         Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
    79         or Nx3 Numeric array of non-negative integers).
     78        Mx2 numeric array of floats) and triangles (sequence of 3-tuples
     79        or Nx3 numeric array of non-negative integers).
    8080        """
    8181
    … …  
    9898
    9999        #Allocate space for geometric quantities
    100         self.centroid_coordinates = num.zeros((N, 2), num.Float)
    101 
    102         self.radii = num.zeros(N, num.Float)
    103 
    104         self.neighbours = num.zeros((N, 3), num.Int)
    105         self.neighbour_edges = num.zeros((N, 3), num.Int)
    106         self.number_of_boundaries = num.zeros(N, num.Int)
    107         self.surrogate_neighbours = num.zeros((N, 3), num.Int)
     100        self.centroid_coordinates = num.zeros((N, 2), num.float)
     101
     102        self.radii = num.zeros(N, num.float)
     103
     104        self.neighbours = num.zeros((N, 3), num.int)
     105        self.neighbour_edges = num.zeros((N, 3), num.int)
     106        self.number_of_boundaries = num.zeros(N, num.int)
     107        self.surrogate_neighbours = num.zeros((N, 3), num.int)
    108108
    109109        #Get x,y coordinates for all triangles and store
    … …  
    124124
    125125            #Compute centroid
    126             centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3], num.Float)
     126            centroid = num.array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3], num.float)
    127127            self.centroid_coordinates[i] = centroid
    128128
    … …  
    133133
    134134                #Midpoints
    135                 m0 = num.array([(x1 + x2)/2, (y1 + y2)/2], num.Float)
    136                 m1 = num.array([(x0 + x2)/2, (y0 + y2)/2], num.Float)
    137                 m2 = num.array([(x1 + x0)/2, (y1 + y0)/2], num.Float)
     135                m0 = num.array([(x1 + x2)/2, (y1 + y2)/2], num.float)
     136                m1 = num.array([(x0 + x2)/2, (y0 + y2)/2], num.float)
     137                m2 = num.array([(x1 + x0)/2, (y1 + y0)/2], num.float)
    138138
    139139                #The radius is the distance from the centroid of
    … …  
    394394            #Check that all keys in given boundary exist
    395395            for tag in tagged_elements.keys():
    396                 tagged_elements[tag] = num.array(tagged_elements[tag], num.Int)
     396                tagged_elements[tag] = num.array(tagged_elements[tag], num.int)
    397397
    398398                msg = 'Not all elements exist. '
    … …  
    506506            # Sanity check
    507507            if p0 is None:
    508                 raise Exception('Impossible')
    509 
     508                msg = 'Impossible: p0 is None!?'
     509                raise Exception, msg
    510510
    511511            # Register potential paths from A to B
    … …  
    613613                point_registry[tuple(p1)] = len(point_registry)
    614614           
    615             polygon.append(list(p1)) # De-Numeric each point :-)
     615            polygon.append(list(p1)) # De-numeric each point :-)
    616616            p0 = p1
    617617
    … …  
    11391139
    11401140            # Distances from line origin to the two intersections
    1141             z0 = num.array([x0 - xi0, y0 - eta0], num.Float)
    1142             z1 = num.array([x1 - xi0, y1 - eta0], num.Float)
     1141            z0 = num.array([x0 - xi0, y0 - eta0], num.float)
     1142            z1 = num.array([x1 - xi0, y1 - eta0], num.float)
    11431143            d0 = num.sqrt(num.sum(z0**2))
    11441144            d1 = num.sqrt(num.sum(z1**2))
    … …  
    11551155            # Normal direction:
    11561156            # Right hand side relative to line direction
    1157             vector = num.array([x1 - x0, y1 - y0], num.Float) # Segment vector
     1157            vector = num.array([x1 - x0, y1 - y0], num.float) # Segment vector
    11581158            length = num.sqrt(num.sum(vector**2))      # Segment length
    1159             normal = num.array([vector[1], -vector[0]], num.Float)/length
     1159            normal = num.array([vector[1], -vector[0]], num.float)/length
    11601160
    11611161
    … …  
    12371237        assert isinstance(segment, Triangle_intersection), msg
    12381238       
    1239         midpoint = num.sum(num.array(segment.segment, num.Float))/2
     1239        midpoint = num.sum(num.array(segment.segment, num.float))/2
    12401240        midpoints.append(midpoint)
    12411241

  • branches/numpy/anuga/abstract_2d_finite_volumes/pmesh2domain.py

    r6145 r6304  
    1 """Class pmesh2domain - Converting .tsh files to doamains
     1"""Class pmesh2domain - Converting .tsh files to domains
    22
    33
    … …  
    88
    99import sys
    10 
    11 import Numeric as num
    12 
    13 
    14 def pmesh_instance_to_domain_instance(mesh,
    15                                       DomainClass):
    16     """
    17     Convert a pmesh instance/object into a domain instance.
    18 
    19     Use pmesh_to_domain_instance to convert a mesh file to a domain instance.
    20     """
    21 
    22     vertex_coordinates, vertices, tag_dict, vertex_quantity_dict \
    23                         ,tagged_elements_dict, geo_reference = \
    24                         pmesh_to_domain(mesh_instance=mesh)
     10import numpy as num
     11
     12
     13##
     14# @brief Convert a pmesh instance to a domain instance.
     15# @param mesh The pmesh instance to convert.
     16# @param DomainClass The class to instantiate and return.
     17# @return The converted pmesh instance (as a 'DomainClass' instance).
     18def pmesh_instance_to_domain_instance(mesh, DomainClass):
     19    """Convert a pmesh instance/object into a domain instance.
     20
     21    Uses pmesh_to_domain_instance to convert a mesh file to a domain instance.
     22    """
     23
     24    (vertex_coordinates, vertices, tag_dict, vertex_quantity_dict,
     25     tagged_elements_dict, geo_reference) = pmesh_to_domain(mesh_instance=mesh)
     26
     27    # NOTE(Ole): This import cannot be at the module level
     28    #            due to mutual dependency with domain.py
     29    from anuga.abstract_2d_finite_volumes.domain import Domain
     30
     31    # ensure that the required 'DomainClass' actually is an instance of Domain
     32    msg = ('The class %s is not a subclass of the generic domain class %s'
     33           % (DomainClass, Domain))
     34    assert issubclass(DomainClass, Domain), msg
     35
     36    # instantiate the result class
     37    result = DomainClass(coordinates=vertex_coordinates,
     38                         vertices=vertices,
     39                         boundary=tag_dict,
     40                         tagged_elements=tagged_elements_dict,
     41                         geo_reference=geo_reference)
     42
     43    # set the water stage to be the elevation
     44    if (vertex_quantity_dict.has_key('elevation') and
     45        not vertex_quantity_dict.has_key('stage')):
     46        vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation']
     47    result.set_quantity_vertices_dict(vertex_quantity_dict)
     48
     49    return result
     50
     51
     52##
     53# @brief Convert a mesh file to a Domain instance.
     54# @param file_name Name of the file to convert (TSH or MSH).
     55# @param DomainClass Class of return instance.
     56# @param use_cache True if caching is to be used.
     57# @param verbose True if this function is to be verbose.
     58# @return An instance of 'DomainClass' containing the file data.
     59def pmesh_to_domain_instance(file_name, DomainClass, use_cache=False,
     60                             verbose=False):
     61    """Converts a mesh file(.tsh or .msh), to a Domain instance.
     62
     63    file_name is the name of the mesh file to convert, including the extension
     64
     65    DomainClass is the Class that will be returned.
     66    It must be a subclass of Domain, with the same interface as domain.
     67
     68    use_cache: True means that caching is attempted for the computed domain.   
     69    """
     70
     71    if use_cache is True:
     72        from caching import cache
     73        result = cache(_pmesh_to_domain_instance, (file_name, DomainClass),
     74                       dependencies=[file_name], verbose=verbose)
     75    else:
     76        result = apply(_pmesh_to_domain_instance, (file_name, DomainClass))       
     77       
     78    return result
     79
     80
     81##
     82# @brief Convert a mesh file to a Domain instance.
     83# @param file_name Name of the file to convert (TSH or MSH).
     84# @param DomainClass Class of return instance.
     85# @return The DomainClass instance containing the file data.
     86def _pmesh_to_domain_instance(file_name, DomainClass):
     87    """Converts a mesh file(.tsh or .msh), to a Domain instance.
     88
     89    Internal function. See public interface pmesh_to_domain_instance for details
     90    """
     91   
     92    (vertex_coordinates, vertices, tag_dict, vertex_quantity_dict,
     93     tagged_elements_dict, geo_reference) = pmesh_to_domain(file_name=file_name)
    2594
    2695    # NOTE(Ole): This import cannot be at the module level due to mutual
    … …  
    2897    from anuga.abstract_2d_finite_volumes.domain import Domain
    2998
    30  
    31 
    32 
    33     msg = 'The class %s is not a subclass of the generic domain class %s'\
    34           %(DomainClass, Domain)
     99    # ensure the required class is a subclass of Domain
     100    msg = ('The class %s is not a subclass of the generic domain class %s'
     101           % (DomainClass, Domain))
    35102    assert issubclass(DomainClass, Domain), msg
    36 
    37103
    38104    domain = DomainClass(coordinates = vertex_coordinates,
    … …  
    42108                         geo_reference = geo_reference )
    43109
    44     # set the water stage to be the elevation
    45     if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'):
    46         vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation']
    47 
    48     domain.set_quantity_vertices_dict(vertex_quantity_dict)
    49     #print "vertex_quantity_dict",vertex_quantity_dict
    50     return domain
    51 
    52 
    53 
    54 def pmesh_to_domain_instance(file_name, DomainClass, use_cache = False, verbose = False):
    55     """
    56     Converts a mesh file(.tsh or .msh), to a Domain instance.
    57 
    58     file_name is the name of the mesh file to convert, including the extension
    59 
    60     DomainClass is the Class that will be returned.
    61     It must be a subclass of Domain, with the same interface as domain.
    62 
    63     use_cache: True means that caching is attempted for the computed domain.   
    64     """
    65 
    66     if use_cache is True:
    67         from caching import cache
    68         result = cache(_pmesh_to_domain_instance, (file_name, DomainClass),
    69                        dependencies = [file_name],                     
    70                        verbose = verbose)
    71 
    72     else:
    73         result = apply(_pmesh_to_domain_instance, (file_name, DomainClass))       
    74        
    75     return result
    76 
    77 
    78 
    79 
    80 def _pmesh_to_domain_instance(file_name, DomainClass):
    81     """
    82     Converts a mesh file(.tsh or .msh), to a Domain instance.
    83 
    84     Internal function. See public interface pmesh_to_domain_instance for details
    85     """
    86    
    87     vertex_coordinates, vertices, tag_dict, vertex_quantity_dict, \
    88                         tagged_elements_dict, geo_reference = \
    89                         pmesh_to_domain(file_name=file_name)
    90 
    91 
    92     # NOTE(Ole): This import cannot be at the module level due to mutual
    93     # dependency with domain.py
    94     from anuga.abstract_2d_finite_volumes.domain import Domain
    95 
    96 
    97     msg = 'The class %s is not a subclass of the generic domain class %s'\
    98           %(DomainClass, Domain)
    99     assert issubclass(DomainClass, Domain), msg
    100 
    101 
    102 
    103     domain = DomainClass(coordinates = vertex_coordinates,
    104                          vertices = vertices,
    105                          boundary = tag_dict,
    106                          tagged_elements = tagged_elements_dict,
    107                          geo_reference = geo_reference )
    108 
    109 
    110 
    111     #FIXME (Ole): Is this really the right place to apply the a default
    112     #value specific to the shallow water wave equation?
    113     #The 'assert' above indicates that any subclass of Domain is acceptable.
    114     #Suggestion - module shallow_water.py will eventually take care of this
    115     #(when I get around to it) so it should be removed from here.
     110    # FIXME (Ole): Is this really the right place to apply a default
     111    # value specific to the shallow water wave equation?
     112    # The 'assert' above indicates that any subclass of Domain is acceptable.
     113    # Suggestion - module shallow_water.py will eventually take care of this
     114    # (when I get around to it) so it should be removed from here.
    116115
    117116    # This doesn't work on the domain instance.
    118117    # This is still needed so -ve elevations don't cuase 'lakes'
    119118    # The fixme we discussed was to only create a quantity when its values
    120     #are set.
     119    # are set.
    121120    # I think that's the way to go still
    122121
    123122    # set the water stage to be the elevation
    124     if vertex_quantity_dict.has_key('elevation') and not vertex_quantity_dict.has_key('stage'):
     123    if (vertex_quantity_dict.has_key('elevation') and
     124        not vertex_quantity_dict.has_key('stage')):
    125125        vertex_quantity_dict['stage'] = vertex_quantity_dict['elevation']
    126 
    127126    domain.set_quantity_vertices_dict(vertex_quantity_dict)
    128     #print "vertex_quantity_dict",vertex_quantity_dict
     127
    129128    return domain
    130129
    131130
    132 def pmesh_to_domain(file_name=None,
    133                     mesh_instance=None,
    134                     use_cache=False,
     131##
     132# @brief Convert pmesh file/instance to list(s) that can instantiate a Domain.
     133# @param file_name Path to file to convert.
     134# @param mesh_instance Instance to convert.
     135# @param use_cache True if we are to cache.
     136# @param verbose True if this function is to be verbose.
     137# @return ??
     138def pmesh_to_domain(file_name=None, mesh_instance=None, use_cache=False,
    135139                    verbose=False):
    136     """
    137     Convert a pmesh file or a pmesh mesh instance to a bunch of lists
     140    """Convert a pmesh file or a pmesh mesh instance to a bunch of lists
    138141    that can be used to instanciate a domain object.
    139142
    … …  
    144147        from caching import cache
    145148        result = cache(_pmesh_to_domain, (file_name, mesh_instance),
    146                        dependencies = [file_name],                     
    147                        verbose = verbose)
     149                       dependencies=[file_name], verbose=verbose)
    148150
    149151    else:
    … …  
    153155
    154156
    155 def _pmesh_to_domain(file_name=None,
    156                      mesh_instance=None,
    157                      use_cache=False,
     157##
     158# @brief Convert pmesh file/instance to list(s) that can instantiate a Domain.
     159# @param file_name Path to file to convert.
     160# @param mesh_instance Instance to convert.
     161# @param use_cache True if we are to cache.
     162# @param verbose True if this function is to be verbose.
     163# @return ??
     164def _pmesh_to_domain(file_name=None, mesh_instance=None, use_cache=False,
    158165                     verbose=False):
    159     """
    160     Convert a pmesh file or a pmesh mesh instance to a bunch of lists
     166    """Convert a pmesh file or a pmesh mesh instance to a bunch of lists
    161167    that can be used to instanciate a domain object.
    162168    """
    … …  
    164170    from load_mesh.loadASCII import import_mesh_file
    165171
     172    # get data from mesh instance or file
    166173    if file_name is None:
    167174        mesh_dict = mesh_instance.Mesh2IODict()
    168175    else:
    169176        mesh_dict = import_mesh_file(file_name)
    170     #print "mesh_dict",mesh_dict
     177
     178    # extract required data from the mesh dictionary
    171179    vertex_coordinates = mesh_dict['vertices']
    172180    volumes = mesh_dict['triangles']
    173181    vertex_quantity_dict = {}
    174     point_atts = num.transpose(mesh_dict['vertex_attributes'])
     182    point_atts = mesh_dict['vertex_attributes']
    175183    point_titles  = mesh_dict['vertex_attribute_titles']
    176184    geo_reference  = mesh_dict['geo_reference']
    177     if point_atts != None:
    178         for quantity, value_vector in map (None, point_titles, point_atts):
     185    if point_atts is not None:
     186        point_atts = num.transpose(mesh_dict['vertex_attributes'])
     187        for quantity, value_vector in map(None, point_titles, point_atts):
    179188            vertex_quantity_dict[quantity] = value_vector
    180189    tag_dict = pmesh_dict_to_tag_dict(mesh_dict)
    181190    tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict)
    182     return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict, tagged_elements_dict, geo_reference
    183 
     191
     192    return (vertex_coordinates, volumes, tag_dict, vertex_quantity_dict,
     193            tagged_elements_dict, geo_reference)
    184194
    185195
    186196def build_tagged_elements_dictionary(mesh_dict):
    187197    """Build the dictionary of element tags.
     198
    188199    tagged_elements is a dictionary of element arrays,
    189     keyed by tag:
    190     { (tag): [e1, e2, e3..] }
    191     """
     200    keyed by tag: { (tag): [e1, e2, e3..] }
     201    """
     202
    192203    tri_atts = mesh_dict['triangle_tags']
    193204    tagged_elements = {}
    … …  
    198209            tagged_elements.setdefault(tri_atts[tri_att_index],
    199210                                       []).append(tri_att_index)
     211
    200212    return tagged_elements
     213
    201214
    202215def pmesh_dict_to_tag_dict(mesh_dict):
    … …  
    204217    to a dictionary of tags, indexed with volume id and face number.
    205218    """
     219
    206220    triangles = mesh_dict['triangles']
    207221    sides = calc_sides(triangles)
    208222    tag_dict = {}
    209     for seg, tag in map(None, mesh_dict['segments'],
    210                         mesh_dict['segment_tags']):
     223    for seg, tag in map(None, mesh_dict['segments'], mesh_dict['segment_tags']):
    211224        v1 = int(seg[0])
    212225        v2 = int(seg[1])
    … …  
    222235
    223236def calc_sides(triangles):
    224     #Build dictionary mapping from sides (2-tuple of points)
    225     #to left hand side neighbouring triangle
     237    '''Build dictionary mapping from sides (2-tuple of points)
     238    to left hand side neighbouring triangle
     239    '''
     240
    226241    sides = {}
     242
    227243    for id, triangle in enumerate(triangles):
    228244        a = int(triangle[0])
    229245        b = int(triangle[1])
    230246        c = int(triangle[2])
     247
    231248        sides[a,b] = (id, 2) #(id, face)
    232249        sides[b,c] = (id, 0) #(id, face)
    233250        sides[c,a] = (id, 1) #(id, face)
     251
    234252    return sides
     253

  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity.py

    r6244 r6304  
    2323from anuga.caching import cache
    2424
    25 import Numeric as num
     25import anuga.utilities.numerical_tools as aunt
     26
     27import numpy as num
    2628
    2729
    … …  
    4345        if vertex_values is None:
    4446            N = len(domain)             # number_of_elements
    45             self.vertex_values = num.zeros((N, 3), num.Float)
     47            self.vertex_values = num.zeros((N, 3), num.float)
    4648        else:
    47             self.vertex_values = num.array(vertex_values, num.Float)
     49            self.vertex_values = num.array(vertex_values, num.float)
    4850
    4951            N, V = self.vertex_values.shape
    … …  
    5759
    5860        # Allocate space for other quantities
    59         self.centroid_values = num.zeros(N, num.Float)
    60         self.edge_values = num.zeros((N, 3), num.Float)
     61        self.centroid_values = num.zeros(N, num.float)
     62        self.edge_values = num.zeros((N, 3), num.float)
    6163
    6264        # Allocate space for Gradient
    63         self.x_gradient = num.zeros(N, num.Float)
    64         self.y_gradient = num.zeros(N, num.Float)
     65        self.x_gradient = num.zeros(N, num.float)
     66        self.y_gradient = num.zeros(N, num.float)
    6567
    6668        # Allocate space for Limiter Phi
    67         self.phi = num.zeros(N, num.Float)
     69        self.phi = num.zeros(N, num.float)
    6870
    6971        # Intialise centroid and edge_values
    … …  
    7274        # Allocate space for boundary values
    7375        L = len(domain.boundary)
    74         self.boundary_values = num.zeros(L, num.Float)
     76        self.boundary_values = num.zeros(L, num.float)
    7577
    7678        # Allocate space for updates of conserved quantities by
    … …  
    7880
    7981        # Allocate space for update fields
    80         self.explicit_update = num.zeros(N, num.Float )
    81         self.semi_implicit_update = num.zeros(N, num.Float )
    82         self.centroid_backup_values = num.zeros(N, num.Float)
     82        self.explicit_update = num.zeros(N, num.float )
     83        self.semi_implicit_update = num.zeros(N, num.float )
     84        self.centroid_backup_values = num.zeros(N, num.float)
    8385
    8486        self.set_beta(1.0)
    … …  
    316318
    317319        numeric:
    318           Compatible list, Numeric array (see below) or constant.
     320          Compatible list, numeric array (see below) or constant.
    319321          If callable it will treated as a function (see below)
    320322          If instance of another Quantity it will be treated as such.
    … …  
    357359
    358360                  In case of location == 'centroids' the dimension values must
    359                   be a list of a Numerical array of length N,
     361                  be a list of a numerical array of length N,
    360362                  N being the number of elements.
    361363                  Otherwise it must be of dimension Nx3
    … …  
    458460
    459461        msg = 'Indices must be a list or None'
    460         assert type(indices) in [ListType, NoneType, num.ArrayType], msg
     462        assert (type(indices) in [ListType, NoneType]
     463                or isinstance(indices, num.ndarray)), msg
    461464
    462465        # Determine which 'set_values_from_...' to use
    463466        if numeric is not None:
    464             if type(numeric) in [FloatType, IntType, LongType]:
    465                 self.set_values_from_constant(numeric, location,
    466                                               indices, verbose)
    467             elif type(numeric) in [num.ArrayType, ListType]:
     467            if isinstance(numeric, num.ndarray) or isinstance(numeric, list):
    468468                self.set_values_from_array(numeric, location, indices,
    469469                                           use_cache=use_cache, verbose=verbose)
    … …  
    479479                                                     indices, verbose=verbose,
    480480                                                     use_cache=use_cache)
     481#            if (isinstance(numeric, float) or isinstance(numeric, int)
     482#                or isinstance(numeric, long)):
    481483            else:
    482                 msg = 'Illegal type for argument numeric: %s' % str(numeric)
    483                 raise msg
     484                try:
     485                    numeric = float(numeric)
     486                except ValueError:
     487                    msg = ("Illegal type for variable 'numeric': "
     488                           "%s, shape=%s\n%s"
     489                           % (type(numeric), numeric.shape, str(numeric)))
     490                    msg += ('\ntype(numeric)==FloatType is %s'
     491                            % str(type(numeric)==FloatType))
     492                    msg += ('\nisinstance(numeric, float)=%s'
     493                            % str(isinstance(numeric, float)))
     494                    msg += ('\nisinstance(numeric, num.ndarray)=%s'
     495                            % str(isinstance(numeric, num.ndarray)))
     496                    raise Exception, msg
     497                self.set_values_from_constant(numeric, location,
     498                                              indices, verbose)
     499        elif quantity is not None:
     500            if type(numeric) in [FloatType, IntType, LongType]:
     501                self.set_values_from_constant(numeric, location,
     502                                              indices, verbose)
     503            else:
     504                msg = ("Illegal type for variable 'numeric': %s, shape=%s\n%s"
     505                       % (type(numeric), numeric.shape, str(numeric)))
     506                msg += ('\ntype(numeric)==FloatType is %s'
     507                        % str(type(numeric)==FloatType))
     508                msg += ('\nisinstance(numeric, float)=%s'
     509                        % str(isinstance(numeric, float)))
     510                msg += ('\nisinstance(numeric, num.ndarray)=%s'
     511                        % str(isinstance(numeric, num.ndarray)))
     512                raise Exception, msg
    484513        elif quantity is not None:
    485514            self.set_values_from_quantity(quantity, location, indices, verbose)
    … …  
    583612        """Set values for quantity
    584613
    585         values: Numeric array
     614        values: numeric array
    586615        location: Where values are to be stored.
    587616        Permissible options are: vertices, centroid, unique vertices
    … …  
    593622
    594623        In case of location == 'centroid' the dimension values must
    595         be a list of a Numerical array of length N, N being the number
     624        be a list of a numerical array of length N, N being the number
    596625        of elements.
    597626
    … …  
    607636        """
    608637
    609         values = num.array(values, num.Float)
     638        values = num.array(values, num.float)
    610639
    611640        if indices is not None:
    612             indices = num.array(indices, num.Int)
     641            indices = num.array(indices, num.int)
    613642            msg = ('Number of values must match number of indices: You '
    614643                   'specified %d values and %d indices'
    … …  
    637666                    'Values array must be 1d')
    638667
    639             self.set_vertex_values(values.flat, indices=indices,
     668            self.set_vertex_values(values.flatten(), indices=indices,
    640669                                   use_cache=use_cache, verbose=verbose)
    641670        else:
    … …  
    788817        from anuga.coordinate_transforms.geo_reference import Geo_reference
    789818
    790         points = ensure_numeric(points, num.Float)
    791         values = ensure_numeric(values, num.Float)
     819        points = ensure_numeric(points, num.float)
     820        values = ensure_numeric(values, num.float)
    792821
    793822        if location != 'vertices':
    … …  
    10801109
    10811110        # Ensure that interpolation points is either a list of
    1082         # points, Nx2 array, or geospatial and convert to Numeric array
     1111        # points, Nx2 array, or geospatial and convert to numeric array
    10831112        if isinstance(interpolation_points, Geospatial_data):
    10841113            # Ensure interpolation points are in absolute UTM coordinates
    … …  
    11191148        """Get values for quantity
    11201149
    1121         Extract values for quantity as a Numeric array.
     1150        Extract values for quantity as a numeric array.
    11221151
    11231152        Inputs:
    … …  
    11371166
    11381167        In case of location == 'centroids' the dimension of returned
    1139         values will be a list or a Numerical array of length N, N being
     1168        values will be a list or a numerical array of length N, N being
    11401169        the number of elements.
    11411170
    … …  
    11731202                            'edges', 'unique vertices']:
    11741203            msg = 'Invalid location: %s' % location
    1175             raise msg
     1204            raise Exception, msg
    11761205
    11771206        import types
    11781207
    11791208        assert type(indices) in [types.ListType, types.NoneType,
    1180                                  num.ArrayType], \
    1181                    'Indices must be a list or None'
     1209                                 num.ndarray], \
     1210                   'Indices must be a list or None'     #??#
    11821211
    11831212        if location == 'centroids':
    … …  
    12101239                    sum += self.vertex_values[triangle_id, vertex_id]
    12111240                vert_values.append(sum / len(triangles))
    1212             return num.array(vert_values, num.Float)
     1241            return num.array(vert_values, num.float)
    12131242        else:
    12141243            if (indices is None):
    … …  
    12361265
    12371266        # Check that A can be converted to array and is of appropriate dim
    1238         A = ensure_numeric(A, num.Float)
     1267        A = ensure_numeric(A, num.float)
    12391268        assert len(A.shape) == 1
    12401269
    … …  
    13361365
    13371366        if precision is None:
    1338             precision = num.Float
     1367            precision = num.float
    13391368
    13401369        if smooth is True:
    … …  
    13421371            V = self.domain.get_triangles()
    13431372            N = self.domain.number_of_full_nodes # Ignore ghost nodes if any
    1344             A = num.zeros(N, num.Float)
     1373            A = num.zeros(N, num.float)
    13451374            points = self.domain.get_nodes()
    13461375
    … …  
    13601389                    if current_node == N:
    13611390                        msg = 'Current node exceeding number of nodes (%d) ' % N
    1362                         raise msg
     1391                        raise Exception, msg
    13631392
    13641393                    k += 1
    … …  
    13811410            V = self.domain.get_disconnected_triangles()
    13821411            points = self.domain.get_vertex_coordinates()
    1383             A = self.vertex_values.flat.astype(precision)
     1412            A = self.vertex_values.flatten().astype(precision)
    13841413
    13851414        # Return

  • branches/numpy/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r5897 r6304  
    1111
    1212#include "Python.h"
    13 #include "Numeric/arrayobject.h"
     13#include "numpy/arrayobject.h"
    1414#include "math.h"
    1515

  • branches/numpy/anuga/abstract_2d_finite_volumes/region.py

    r6145 r6304  
    88# FIXME (DSG-DSG) add better comments
    99
    10 import Numeric as num
     10import numpy as num
    1111
    1212

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_domain.py

    r6195 r6304  
    77from anuga.config import epsilon
    88
    9 import Numeric as num
     9import numpy as num
    1010
    1111
    … …  
    6565
    6666
    67         assert domain.get_conserved_quantities(0, edge=1) == 0.
     67        assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
    6868
    6969
    … …  
    345345
    346346
    347         A = num.array([[1,2,3], [5,5,-5], [0,0,9], [-6,3,3]], num.Float)
    348         B = num.array([[2,4,4], [3,2,1], [6,-3,4], [4,5,-1]], num.Float)
     347        A = num.array([[1,2,3], [5,5,-5], [0,0,9], [-6,3,3]], num.float)
     348        B = num.array([[2,4,4], [3,2,1], [6,-3,4], [4,5,-1]], num.float)
    349349       
    350350        #print A
    … …  
    613613
    614614        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
    615         denom = num.ones(4, num.Float)-domain.timestep*sem
     615        denom = num.ones(4, num.float)-domain.timestep*sem
    616616
    617617#        x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] )

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_ermapper.py

    r6145 r6304  
    77from os import remove
    88
    9 import Numeric as num
     9import numpy as num
    1010
    1111
    … …  
    4545
    4646        # Write test data
    47         ermapper_grids.write_ermapper_data(original_grid, filename, num.Float64)
     47        ermapper_grids.write_ermapper_data(original_grid, filename, num.float64)
    4848
    4949        # Read in the test data
    50         new_grid = ermapper_grids.read_ermapper_data(filename, num.Float64)
     50        new_grid = ermapper_grids.read_ermapper_data(filename, num.float64)
    5151
    5252        # Check that the test data that has been read in matches the original data

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_general_mesh.py

    r6174 r6304  
    1010from anuga.coordinate_transforms.geo_reference import Geo_reference
    1111
    12 import Numeric as num
     12import numpy as num
    1313
    1414
    … …  
    5959       
    6060        #bac, bce, ecf, dbe, daf, dae
    61         triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.Int)
     61        triangles = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.int)
    6262
    6363        domain = General_mesh(nodes, triangles,
    … …  
    277277                       geo_reference = geo)
    278278        node = domain.get_node(2)       
    279         self.assertEqual(c, node)
     279        #self.assertEqual(c, node)
     280        self.failUnless(num.alltrue(c == node))
    280281       
    281282        node = domain.get_node(2, absolute=True)     
    282         self.assertEqual(nodes_absolute[2], node)
     283        #self.assertEqual(nodes_absolute[2], node)
     284        self.failUnless(num.alltrue(nodes_absolute[2] == node))
    283285       
    284286        node = domain.get_node(2, absolute=True)     
    285         self.assertEqual(nodes_absolute[2], node)
     287        #self.assertEqual(nodes_absolute[2], node)
     288        self.failUnless(num.alltrue(nodes_absolute[2] == node))
    286289       
    287290

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_generic_boundary_conditions.py

    r6195 r6304  
    88from anuga.config import epsilon
    99
    10 import Numeric as num
     10import numpy as num
    1111
    1212

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_ghost.py

    r6195 r6304  
    66from anuga.abstract_2d_finite_volumes.domain import *
    77from anuga.config import epsilon
     8
     9import numpy as num
    810
    911
    … …  
    4042
    4143
    42         assert domain.get_conserved_quantities(0, edge=1) == 0.
     44        assert num.alltrue(domain.get_conserved_quantities(0, edge=1) == 0.)
    4345
    4446

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_neighbour_mesh.py

    r6174 r6304  
    1818from anuga.utilities.numerical_tools import ensure_numeric
    1919
    20 import Numeric as num
     20import numpy as num
    2121
    2222
    … …  
    988988                  [  75735.4765625 ,  23762.00585938],
    989989                  [  52341.70703125,  38563.39453125]]
    990 
    991         ##points = ensure_numeric(points, Int)/1000  # Simplify for ease of interpretation       
    992990
    993991        triangles = [[19, 0,15],
    … …  
    10921090                  [  35406.3359375 ,  79332.9140625 ]]
    10931091
    1094         scaled_points = ensure_numeric(points, num.Int)/1000  # Simplify for ease of interpretation
     1092        scaled_points = ensure_numeric(points, num.int)/1000  # Simplify for ease of interpretation
    10951093
    10961094        triangles = [[ 0, 1, 2],

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_pmesh2domain.py

    r6191 r6304  
    1414from anuga.pmesh.mesh import importMeshFromFile
    1515
    16 import Numeric as num
     16import numpy as num
    1717
    1818

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_quantity.py

    r6195 r6304  
    1515from anuga.utilities.polygon import *
    1616
    17 import Numeric as num
     17import numpy as num
    1818
    1919
    … …  
    17651765        quantity = Quantity(self.mesh4)
    17661766
    1767         quantity.vertex_values = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.Float)
     1767        quantity.vertex_values = num.array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]], num.float)
    17681768
    17691769        quantity.interpolate_from_vertices_to_edges()
    … …  
    17811781                                          [3., 2.5, 1.5],
    17821782                                          [3.5, 4.5, 3.],
    1783                                           [2.5, 3.5, 2]], num.Float)
     1783                                          [2.5, 3.5, 2]], num.float)
    17841784
    17851785        quantity.interpolate_from_edges_to_vertices()
    … …  
    18351835
    18361836        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
    1837         denom = num.ones(4, num.Float)-timestep*sem
     1837        denom = num.ones(4, num.float)-timestep*sem
    18381838
    18391839        x = num.array([1, 2, 3, 4])/denom
    … …  
    18591859
    18601860        sem = num.array([1.,1.,1.,1.])/num.array([1, 2, 3, 4])
    1861         denom = num.ones(4, num.Float)-timestep*sem
     1861        denom = num.ones(4, num.float)-timestep*sem
    18621862
    18631863        x = num.array([1., 2., 3., 4.])
    … …  
    19011901
    19021902        bed = domain.quantities['elevation'].vertex_values
    1903         stage = num.zeros(bed.shape, num.Float)
     1903        stage = num.zeros(bed.shape, num.float)
    19041904
    19051905        h = 0.03
    … …  
    19851985
    19861986        bed = domain.quantities['elevation'].vertex_values
    1987         stage = num.zeros(bed.shape, num.Float)
     1987        stage = num.zeros(bed.shape, num.float)
    19881988
    19891989        h = 0.03
    … …  
    19991999        stage = domain.quantities['stage']
    20002000        A, V = stage.get_vertex_values(xy=False, smooth=False)
    2001         Q = stage.vertex_values.flat
     2001        Q = stage.vertex_values.flatten()
    20022002
    20032003        for k in range(8):
    … …  
    25112511    suite = unittest.makeSuite(Test_Quantity, 'test')   
    25122512    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')
    2513 
    2514     #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_subset_and_geo')
    2515     #print "restricted test"
    2516     #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
    25172513    runner = unittest.TextTestRunner()
    25182514    runner.run(suite)

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_region.py

    r6145 r6304  
    88#from anuga.config import epsilon
    99
    10 import Numeric as num
     10import numpy as num
    1111
    1212
    … …  
    4242
    4343    def test_region_tags(self):
    44         """
    45         get values based on triangle lists.
    46         """
    47         from mesh_factory import rectangular
    48         from shallow_water import Domain
    49 
    50         #Create basic mesh
    51         points, vertices, boundary = rectangular(1, 3)
    52 
    53         #Create shallow water domain
    54         domain = Domain(points, vertices, boundary)
    55         domain.build_tagged_elements_dictionary({'bottom':[0,1],
    56                                                  'top':[4,5],
    57                                                  'all':[0,1,2,3,4,5]})
    58 
     44        """get values based on triangle lists."""
     45
     46        from mesh_factory import rectangular
     47        from shallow_water import Domain
     48
     49        #Create basic mesh
     50        points, vertices, boundary = rectangular(1, 3)
     51
     52        #Create shallow water domain
     53        domain = Domain(points, vertices, boundary)
     54        domain.build_tagged_elements_dictionary({'bottom': [0,1],
     55                                                 'top': [4,5],
     56                                                 'all': [0,1,2,3,4,5]})
    5957
    6058        #Set friction
    … …  
    6563        b = Set_region('top', 'friction', 1.0)
    6664        domain.set_region([a, b])
    67         #print domain.quantities['friction'].get_values()
    68         assert num.allclose(domain.quantities['friction'].get_values(),\
    69                             [[ 0.09,  0.09,  0.09],
    70                              [ 0.09,  0.09,  0.09],
    71                              [ 0.07,  0.07,  0.07],
    72                              [ 0.07,  0.07,  0.07],
    73                              [ 1.0,  1.0,  1.0],
    74                              [ 1.0,  1.0,  1.0]])
     65
     66        expected = [[ 0.09,  0.09,  0.09],
     67                    [ 0.09,  0.09,  0.09],
     68                    [ 0.07,  0.07,  0.07],
     69                    [ 0.07,  0.07,  0.07],
     70                    [ 1.0,  1.0,  1.0],
     71                    [ 1.0,  1.0,  1.0]]
     72        msg = ("\ndomain.quantities['friction']=%s\nexpected value=%s"
     73               % (str(domain.quantities['friction'].get_values()),
     74                  str(expected)))
     75        assert num.allclose(domain.quantities['friction'].get_values(),
     76                            expected), msg
    7577
    7678        #c = Add_Value_To_region('all', 'friction', 10.0)
    … …  
    126128       
    127129    def test_unique_vertices(self):
    128         """
    129         get values based on triangle lists.
    130         """
     130        """get values based on triangle lists."""
     131
    131132        from mesh_factory import rectangular
    132133        from shallow_water import Domain
    … …  
    147148        a = Set_region('bottom', 'friction', 0.09, location = 'unique vertices')
    148149        domain.set_region(a)
    149         #print domain.quantities['friction'].get_values()
    150         assert num.allclose(domain.quantities['friction'].get_values(),\
     150        assert num.allclose(domain.quantities['friction'].get_values(),
    151151                            [[ 0.09,  0.09,  0.09],
    152152                             [ 0.09,  0.09,  0.09],
    … …  
    217217        #print domain.quantities['friction'].get_values()
    218218        frict_points = domain.quantities['friction'].get_values()
    219         assert num.allclose(frict_points[0],\
    220                             [ calc_frict, calc_frict, calc_frict])
    221         assert num.allclose(frict_points[1],\
    222                             [ calc_frict, calc_frict, calc_frict])
     219        expected = [calc_frict, calc_frict, calc_frict]
     220        msg = ('frict_points[0]=%s\nexpected=%s' % (str(frict_points[0]),
     221                                                    str(expected)))
     222        assert num.allclose(frict_points[0], expected), msg
     223        msg = ('frict_points[1]=%s\nexpected=%s' % (str(frict_points[1]),
     224                                                    str(expected)))
     225        assert num.allclose(frict_points[1], expected), msg
    223226 
    224227    def test_unique_vertices_average_loc_unique_vert(self):
    … …  
    262265#-------------------------------------------------------------
    263266if __name__ == "__main__":
    264     suite = unittest.makeSuite(Test_Region,'test')
     267    suite = unittest.makeSuite(Test_Region, 'test')   
    265268    runner = unittest.TextTestRunner()
    266269    runner.run(suite)

  • branches/numpy/anuga/abstract_2d_finite_volumes/test_util.py

    r6175 r6304  
    2222import string
    2323
    24 import Numeric as num
     24import numpy as num
    2525
    2626
    … …  
    630630                    q1 = F(t+60, point_id=id)
    631631
    632                 if q0 == NAN:
     632                if num.alltrue(q0 == NAN):
    633633                    actual = q0
    634634                else:
    … …  
    641641                #print "actual", actual
    642642                #print
    643                 if q0 == NAN:
    644                      self.failUnless( q == actual, 'Fail!')
     643                if num.alltrue(q0 == NAN):
     644                     self.failUnless(num.alltrue(q == actual), 'Fail!')
    645645                else:
    646646                    assert num.allclose(q, actual)
    … …  
    11811181        #FIXME: Division is not expected to work for integers.
    11821182        #This must be caught.
    1183         foo = num.array([[1,2,3], [4,5,6]], num.Float)
    1184 
    1185         bar = num.array([[-1,0,5], [6,1,1]], num.Float)                 
     1183        foo = num.array([[1,2,3], [4,5,6]], num.float)
     1184
     1185        bar = num.array([[-1,0,5], [6,1,1]], num.float)                 
    11861186
    11871187        D = {'X': foo, 'Y': bar}
    … …  
    12011201
    12021202        # make an error for zero on zero
    1203         # this is really an error in Numeric, SciPy core can handle it
     1203        # this is really an error in numeric, SciPy core can handle it
    12041204        # Z = apply_expression_to_dictionary('0/Y', D)
    12051205

  • branches/numpy/anuga/abstract_2d_finite_volumes/util.py

    r6194 r6304  
    2727from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    2828
    29 import Numeric as num
     29import numpy as num
    3030
    3131
    … …  
    554554                expression e.g. by overloading
    555555
    556     Due to a limitation with Numeric, this can not evaluate 0/0
     556    Due to a limitation with numeric, this can not evaluate 0/0
    557557    In general, the user can fix by adding 1e-30 to the numerator.
    558558    SciPy core can handle this situation.
    … …  
    789789                        - easting, northing, name , elevation?
    790790                    - OR (this is not yet done)
    791                         - structure which can be converted to a Numeric array,
     791                        - structure which can be converted to a numeric array,
    792792                          such as a geospatial data object
    793793                     
    … …  
    12551255    n0 = int(n0)
    12561256    m = len(locations)
    1257     model_time = num.zeros((n0, m, p), num.Float)
    1258     stages = num.zeros((n0, m, p), num.Float)
    1259     elevations = num.zeros((n0, m, p), num.Float)
    1260     momenta = num.zeros((n0, m, p), num.Float)
    1261     xmom = num.zeros((n0, m, p), num.Float)
    1262     ymom = num.zeros((n0, m, p), num.Float)
    1263     speed = num.zeros((n0, m, p), num.Float)
    1264     bearings = num.zeros((n0, m, p), num.Float)
    1265     due_east = 90.0*num.ones((n0, 1), num.Float)
    1266     due_west = 270.0*num.ones((n0, 1), num.Float)
    1267     depths = num.zeros((n0, m, p), num.Float)
    1268     eastings = num.zeros((n0, m, p), num.Float)
     1257    model_time = num.zeros((n0, m, p), num.float)
     1258    stages = num.zeros((n0, m, p), num.float)
     1259    elevations = num.zeros((n0, m, p), num.float)
     1260    momenta = num.zeros((n0, m, p), num.float)
     1261    xmom = num.zeros((n0, m, p), num.float)
     1262    ymom = num.zeros((n0, m, p), num.float)
     1263    speed = num.zeros((n0, m, p), num.float)
     1264    bearings = num.zeros((n0, m, p), num.float)
     1265    due_east = 90.0*num.ones((n0, 1), num.float)
     1266    due_west = 270.0*num.ones((n0, 1), num.float)
     1267    depths = num.zeros((n0, m, p), num.float)
     1268    eastings = num.zeros((n0, m, p), num.float)
    12691269    min_stages = []
    12701270    max_stages = []
    … …  
    12781278    min_speeds = []   
    12791279    max_depths = []
    1280     model_time_plot3d = num.zeros((n0, m), num.Float)
    1281     stages_plot3d = num.zeros((n0, m), num.Float)
    1282     eastings_plot3d = num.zeros((n0, m),num.Float)
     1280    model_time_plot3d = num.zeros((n0, m), num.float)
     1281    stages_plot3d = num.zeros((n0, m), num.float)
     1282    eastings_plot3d = num.zeros((n0, m),num.float)
    12831283    if time_unit is 'mins': scale = 60.0
    12841284    if time_unit is 'hours': scale = 3600.0
    … …  
    18321832    """
    18331833       
    1834     xc = num.zeros(triangles.shape[0], num.Float) # Space for centroid info
     1834    xc = num.zeros(triangles.shape[0], num.float) # Space for centroid info
    18351835   
    18361836    for k in range(triangles.shape[0]):
    … …  
    21192119                #add tide to stage if provided
    21202120                if quantity == 'stage':
    2121                      quantity_value[quantity] = num.array(quantity_value[quantity], num.Float) \
     2121                     quantity_value[quantity] = num.array(quantity_value[quantity], num.float) \
    21222122                                                          + directory_add_tide
    21232123
    … …  
    24692469       
    24702470    #convert to array for file_function
    2471     points_array = num.array(points,num.Float)
     2471    points_array = num.array(points,num.float)
    24722472       
    24732473    points_array = ensure_absolute(points_array)

  • branches/numpy/anuga/advection/advection.py

    r6146 r6304  
    3333from anuga.abstract_2d_finite_volumes.domain import *
    3434
    35 import Numeric as num
     35import numpy as num
    3636
    3737
    … …  
    252252        stage_bdry = Stage.boundary_values
    253253
    254         flux = num.zeros(1, num.Float) #Work array for summing up fluxes
     254        flux = num.zeros(1, num.float) #Work array for summing up fluxes
    255255
    256256        #Loop

  • branches/numpy/anuga/advection/advection_ext.c

    r5162 r6304  
    1010
    1111#include "Python.h"
    12 #include "Numeric/arrayobject.h"
     12#include "numpy/arrayobject.h"
    1313#include "math.h"
    1414#include "stdio.h"

  • branches/numpy/anuga/advection/test_advection.py

    r6146 r6304  
    88from anuga.advection.advection import Domain, Transmissive_boundary, Dirichlet_boundary
    99
    10 import Numeric as num
     10import numpy as num
    1111
    1212
    … …  
    142142
    143143        X = domain.quantities['stage'].explicit_update
     144        #        print 'X=%s' % str(X)
    144145        assert X[0] == -X[1]
    145146

  • branches/numpy/anuga/alpha_shape/alpha_shape.py

    r6174 r6304  
    2525from anuga.geospatial_data.geospatial_data import Geospatial_data
    2626
    27 import Numeric as num
     27import numpy as num
    2828
    2929
    … …  
    8888                raise PointError, "Three points on a straight line"
    8989       
    90         #Convert input to Numeric arrays
    91         self.points = num.array(points, num.Float)
     90        #Convert input to numeric arrays
     91        self.points = num.array(points, num.float)
    9292
    9393   
    … …  
    288288            zeroind = [k for k in range(len(denom)) if \
    289289                       (denom[k]< EPSILON and  denom[k] > -EPSILON)]
    290         try:
    291             dx = num.divide_safe(y31*dist21 - y21*dist31,denom)
    292             dy = num.divide_safe(x21*dist31 - x31*dist21,denom)
    293         except ZeroDivisionError:
    294             raise  AlphaError
    295            
     290
     291        if num.any(denom == 0.0):
     292            raise AlphaError
     293
     294        dx = num.divide(y31*dist21 - y21*dist31, denom)
     295        dy = num.divide(x21*dist31 - x31*dist21, denom)
     296
    296297        self.triradius = 0.5*num.sqrt(dx*dx + dy*dy)
    297298        #print "triangle radii", self.triradius

  • branches/numpy/anuga/alpha_shape/test_alpha_shape.py

    r6147 r6304  
    55import unittest
    66
    7 import Numeric as num
     7import numpy as num
    88
    99try:

  • branches/numpy/anuga/caching/caching.py

    r6232 r6304  
    5151  unix = True
    5252
    53 import Numeric as num
     53import numpy as num
    5454
    5555
    … …  
    13901390      I.sort()   
    13911391      val = myhash(I, ids)
    1392   elif type(T) == num.ArrayType:
     1392  elif isinstance(T, num.ndarray):
    13931393      T = num.array(T) # Ensure array is contiguous
    13941394
    … …  
    14611461            identical = compare(a, b, ids)
    14621462           
    1463     elif type(A) == num.ArrayType:
     1463    elif isinstance(A, num.ndarray):
    14641464        # Use element by element comparison
    14651465        identical = num.alltrue(A==B)
    … …  
    24352435      argstr = argstr + "'"+str(args)+"'"
    24362436    else:
    2437       # Truncate large Numeric arrays before using str()
    2438       if type(args) == num.ArrayType:
     2437      # Truncate large numeric arrays before using str()
     2438      if isinstance(args, num.ndarray):
    24392439#        if len(args.flat) > textwidth: 
    24402440#        Changed by Duncan and Nick 21/2/07 .flat has problems with

  • branches/numpy/anuga/caching/test_caching.py

    r6232 r6304  
    77from anuga.caching.dummy_classes_for_testing import Dummy, Dummy_memorytest
    88
    9 import Numeric as num
     9import numpy as num
    1010
    1111
    … …  
    2626 
    2727def f_numeric(A, B):
    28   """Operation on Numeric arrays
     28  """Operation on numeric arrays
    2929  """
    3030 
    … …  
    123123        """test_caching_of_numeric_arrays
    124124       
    125         Test that Numeric arrays can be recognised by caching even if their id's are different
     125        Test that numeric arrays can be recognised by caching even if their id's are different
    126126        """
    127127       
    … …  
    160160
    161161
    162             assert T1 == T2, 'Cached result does not match computed result'
    163             assert T2 == T3, 'Cached result does not match computed result'
     162            assert num.alltrue(T1 == T2), 'Cached result does not match computed result'
     163            assert num.alltrue(T2 == T3), 'Cached result does not match computed result'
    164164           
    165165
    166166    def test_hash_collision(self):
    167         """test_hash_collision(self):
    168        
    169         Test that hash collisons are dealt with correctly
    170         """
     167        """Test that hash collisons are dealt with correctly"""
    171168       
    172169        verbose = False
    … …  
    202199            T2 = cache(f_numeric, (A1, A1),
    203200                       compression=comp, verbose=verbose)
    204            
    205 
    206             #print T1
    207             #print T2
    208             assert T2 != T1
    209 
    210            
    211 
    212            
    213            
     201           
     202            T1_ref = f_numeric(A0, A0)
     203            T2_ref = f_numeric(A1, A1)
     204
     205            assert num.alltrue(T1 == T1_ref)
     206            assert num.alltrue(T2 == T2_ref)
     207
     208
    214209    def test_caching_of_dictionaries(self):
    215210        """test_caching_of_dictionaries
    … …  
    412407
    413408       
    414         x = num.arange(10).astype(num.Float)
     409        x = num.arange(10).astype(num.float)
    415410       
    416411        ref1 = f1(x)

  • branches/numpy/anuga/config.py

    r6108 r6304  
    44import os
    55import sys
    6 
    7 ################################################################################
    8 # Numerical constants
     6import numpy as num
     7
     8
     9################################################################################
     10# numerical constants
    911################################################################################
    1012
    … …  
    160162                             # Too large (100) creates 'flopping' water
    161163                             # Too small (0) creates 'creep'
    162                            
     164
    163165maximum_froude_number = 100.0 # To be used in limiters.
    164166
    … …  
    181183
    182184################################################################################
     185# NetCDF-specific type constants.  Used when defining NetCDF file variables.
     186################################################################################
     187
     188netcdf_char = 'c'
     189netcdf_byte = 'b'
     190netcdf_int = 'i'
     191netcdf_float = 'd'
     192netcdf_float64 = 'd'
     193netcdf_float32 = 'f'
     194
     195################################################################################
    183196# Dynamically-defined constants.
    184197################################################################################
    … …  
    191204# Code to set the write mode depending on
    192205# whether Scientific.IO supports large NetCDF files
    193 s = """from Scientific.IO.NetCDF import NetCDFFile; fid = NetCDFFile('tmpfilenamexx', 'wl')"""
     206s = """
     207from Scientific.IO.NetCDF import NetCDFFile
     208fid = NetCDFFile('tmpfilenamexx', 'wl')
     209"""
    194210
    195211# Need to run in a separate process due an

  • branches/numpy/anuga/coordinate_transforms/geo_reference.py

    r6149 r6304  
    1111from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, ParsingError, \
    1212     ShapeError
    13 
    14 import Numeric as num
     13from anuga.config import netcdf_float, netcdf_int, netcdf_float32
     14
     15import numpy as num
    1516
    1617
    … …  
    6465        self.xllcorner = xllcorner
    6566        self.yllcorner = yllcorner       
     67        #self.is_absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
    6668           
    6769        if NetCDFObject is not None:
    … …  
    99101       
    100102        # Fix some assertion failures
    101         if type(self.zone) == num.ArrayType and self.zone.shape == ():
     103        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
    102104            self.zone = self.zone[0]
    103         if type(self.xllcorner) == num.ArrayType and self.xllcorner.shape == ():
     105        if isinstance(self.xllcorner, num.ndarray) and self.xllcorner.shape == ():
    104106            self.xllcorner = self.xllcorner[0]
    105         if type(self.yllcorner) == num.ArrayType and self.yllcorner.shape == ():
     107        if isinstance(self.yllcorner, num.ndarray) and self.yllcorner.shape == ():
    106108            self.yllcorner = self.yllcorner[0]
    107109
    108         assert (type(self.xllcorner) == types.FloatType or\
    109                 type(self.xllcorner) == types.IntType)
    110         assert (type(self.yllcorner) == types.FloatType or\
    111                 type(self.yllcorner) == types.IntType)
    112         assert (type(self.zone) == types.IntType)
     110        assert (self.xllcorner.dtype.kind in num.typecodes['Float'] or
     111                self.xllcorner.dtype.kind in num.typecodes['Integer'])
     112        assert (self.yllcorner.dtype.kind in num.typecodes['Float'] or
     113                self.yllcorner.dtype.kind in num.typecodes['Integer'])
     114        assert (self.zone.dtype.kind in num.typecodes['Integer'])
    113115       
    114116        try:
    … …  
    173175           
    174176        # Fix some assertion failures
    175         if(type(self.zone) == num.ArrayType and self.zone.shape == ()):
     177        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
    176178            self.zone = self.zone[0]
    177         if type(self.xllcorner) == num.ArrayType and self.xllcorner.shape == ():
     179        if isinstance(self.xllcorner, num.ndarray) and self.xllcorner.shape == ():
    178180            self.xllcorner = self.xllcorner[0]
    179         if type(self.yllcorner) == num.ArrayType and self.yllcorner.shape == ():
     181        if isinstance(self.yllcorner, num.ndarray) and self.yllcorner.shape == ():
    180182            self.yllcorner = self.yllcorner[0]
    181    
    182         assert (type(self.xllcorner) == types.FloatType)
    183         assert (type(self.yllcorner) == types.FloatType)
    184         assert (type(self.zone) == types.IntType)
    185        
     183
     184# useless asserts - see try/except code above
     185#        assert (type(self.xllcorner) == types.FloatType)
     186#        assert (type(self.yllcorner) == types.FloatType)
     187#        assert (type(self.zone) == types.IntType)
     188
    186189       
    187190    def change_points_geo_ref(self, points,
    188191                              points_geo_ref=None):
    189192        """
    190         Change the geo reference of a list or Numeric array of points to
     193        Change the geo reference of a list or numeric array of points to
    191194        be this reference.(The reference used for this object)
    192195        If the points do not have a geo ref, assume 'absolute' values
    … …  
    197200            is_list = True
    198201
    199         points = ensure_numeric(points, num.Float)
     202        points = ensure_numeric(points, num.float)
    200203       
    201204        if len(points.shape) == 1:
    … …  
    243246
    244247       
    245    
     248    ##
     249    # @brief
     250    # @param points
     251    # @return
     252    # @note
    246253    def get_absolute(self, points):
    247         """
    248         Given a set of points geo referenced to this instance,
     254        """Given a set of points geo referenced to this instance,
    249255        return the points as absolute values.
    250256        """
    251257
    252         #if self.is_absolute():
     258        #if self.is_absolute:
    253259        #    return points
     260
    254261        is_list = False
    255262        if type(points) == types.ListType:
    256263            is_list = True
    257264
    258         points = ensure_numeric(points, num.Float)
     265        points = ensure_numeric(points, num.float)
    259266        if len(points.shape) == 1:
    260267            #One point has been passed
    … …  
    264271                #points = reshape(points, (1,2))
    265272
    266 
    267273        msg = 'Input must be an N x 2 array or list of (x,y) values. '
    268274        msg += 'I got an %d x %d array' %points.shape   
    269275        if not points.shape[1] == 2:
    270276            raise ShapeError, msg   
    271            
    272277       
    273278        # Add geo ref to points
     279        #if not self.is_absolute:
    274280        if not self.is_absolute():
    275281            points[:,0] += self.xllcorner
    276282            points[:,1] += self.yllcorner
    277 
     283            #self.is_absolute = True
    278284       
    279285        if is_list:
    … …  
    296302            is_list = True
    297303
    298         points = ensure_numeric(points, num.Float)
     304        points = ensure_numeric(points, num.float)
    299305        if len(points.shape) == 1:
    300306            #One point has been passed

  • branches/numpy/anuga/coordinate_transforms/test_geo_reference.py

    r6149 r6304  
    99from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
    1010
    11 import Numeric as num
     11import numpy as num
    1212
    1313
    … …  
    171171        #print "4 new_lofl",new_lofl
    172172
    173         self.failUnless(type(new_lofl) == num.ArrayType, ' failed')
     173        self.failUnless(isinstance(new_lofl, num.ndarray), ' failed')
    174174        self.failUnless(type(new_lofl) == type(lofl), ' failed')
    175175        lofl[:,0] -= x
    … …  
    189189        #print "5 new_lofl",new_lofl
    190190
    191         self.failUnless(type(new_lofl) == num.ArrayType, ' failed')
     191        self.failUnless(isinstance(new_lofl, num.ndarray), ' failed')
    192192        self.failUnless(type(new_lofl) == type(lofl), ' failed')
    193193
    … …  
    206206        #print "new_lofl",new_lofl
    207207
    208         self.failUnless(type(new_lofl) == num.ArrayType, ' failed')
     208        self.failUnless(isinstance(new_lofl, num.ndarray), ' failed')
    209209        self.failUnless(type(new_lofl) == type(lofl), ' failed')
    210210        for point,new_point in map(None,[lofl],new_lofl):
    … …  
    229229            self.failUnless(point[0]+point_x-x==new_point[0], ' failed')
    230230            self.failUnless(point[1]+point_y-y==new_point[1], ' failed')
    231      
    232 
    233     def test_get_absolute(self):
     231
     232    def test_get_absolute_list(self):
     233        # test with supplied offsets
    234234        x = 7.0
    235235        y = 3.0
    236236       
    237         g = Geo_reference(56,x,y)
    238         lofl = [[3.0,34.0], [64.0,6.0]]
    239         new_lofl = g.get_absolute(lofl)
    240         #print "lofl",lofl
    241         #print "new_lofl",new_lofl
    242 
    243         self.failUnless(type(new_lofl) == types.ListType, ' failed')
    244         self.failUnless(type(new_lofl) == type(lofl), ' failed')
    245         for point,new_point in map(None,lofl,new_lofl):
    246             self.failUnless(point[0]+x==new_point[0], ' failed')
    247             self.failUnless(point[1]+y==new_point[1], ' failed')
    248 
     237        g = Geo_reference(56, x, y)
     238        points = [[3.0,34.0], [64.0,6.0]]
     239        new_points = g.get_absolute(points)
     240
     241        self.failUnless(type(new_points) == types.ListType, 'failed')
     242        self.failUnless(type(new_points) == type(points), 'failed')
     243        for point, new_point in map(None, points, new_points):
     244            self.failUnless(point[0]+x == new_point[0], 'failed')
     245            self.failUnless(point[1]+y == new_point[1], 'failed')
     246
     247        # test with no supplied offsets
     248        g = Geo_reference()
     249        points = [[3.0,34.0], [64.0,6.0]]
     250        new_points = g.get_absolute(points)
     251
     252        self.failUnless(type(new_points) == types.ListType, 'failed')
     253        self.failUnless(type(new_points) == type(points), 'failed')
     254        for point, new_point in map(None, points, new_points):
     255            self.failUnless(point[0] == new_point[0], 'failed')
     256            self.failUnless(point[1] == new_point[1], 'failed')
    249257           
     258        # test that calling get_absolute twice does the right thing
     259        # first call
     260        dx = 10.0
     261        dy = 10.0
     262        g = Geo_reference(56, dx, dy)
     263        points = [[3.0,34.0], [64.0,6.0]]
     264        expected_new_points = [[3.0+dx,34.0+dy], [64.0+dx,6.0+dy]]
     265        new_points = g.get_absolute(points)
     266
     267        self.failUnless(type(new_points) == types.ListType, 'failed')
     268        self.failUnless(type(new_points) == type(points), 'failed')
     269        for point, new_point in map(None, expected_new_points, new_points):
     270            self.failUnless(point[0] == new_point[0], 'failed')
     271            self.failUnless(point[1] == new_point[1], 'failed')
     272
     273        # and repeat from 'new_points = g.get_absolute(points)' above
     274        # to see if second call with same input gives same results.
     275        new_points = g.get_absolute(points)
     276
     277        self.failUnless(type(new_points) == types.ListType, 'failed')
     278        self.failUnless(type(new_points) == type(points), 'failed')
     279        for point, new_point in map(None, expected_new_points, new_points):
     280            self.failUnless(point[0] == new_point[0], 'failed')
     281            self.failUnless(point[1] == new_point[1], 'failed')
     282
     283    def test_get_absolute_array(self):
     284        '''Same test as test_get_absolute_list(), but with numeric arrays.'''
     285
     286        # test with supplied offsets
     287        x = 7.0
     288        y = 3.0
     289       
     290        g = Geo_reference(56, x, y)
     291        points = num.array([[3.0,34.0], [64.0,6.0]])
     292        new_points = g.get_absolute(points)
     293
     294        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     295        self.failUnless(type(new_points) == type(points), 'failed')
     296        self.failUnless(num.alltrue(points == new_points), 'failed')
     297
     298        # test with no supplied offsets
    250299        g = Geo_reference()
    251         lofl = [[3.0,34.0], [64.0,6.0]]
    252         new_lofl = g.get_absolute(lofl)
    253         #print "lofl",lofl
    254         #print "new_lofl",new_lofl
    255 
    256         self.failUnless(type(new_lofl) == types.ListType, ' failed')
    257         self.failUnless(type(new_lofl) == type(lofl), ' failed')
    258         for point,new_point in map(None,lofl,new_lofl):
    259             self.failUnless(point[0]==new_point[0], ' failed')
    260             self.failUnless(point[1]==new_point[1], ' failed')
    261            
     300        points = num.array([[3.0,34.0], [64.0,6.0]])
     301        new_points = g.get_absolute(points)
     302
     303        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     304        self.failUnless(type(new_points) == type(points), 'failed')
     305        self.failUnless(num.alltrue(points == new_points), 'failed')
     306
     307        # test that calling get_absolute twice does the right thing
     308        # first call
     309        dx = 10.0
     310        dy = 10.0
     311        g = Geo_reference(56, dx, dy)
     312        points = num.array([[3.0,34.0], [64.0,6.0]])
     313        expected_new_points = num.array([[3.0+dx,34.0+dy], [64.0+dx,6.0+dy]])
     314        new_points = g.get_absolute(points)
     315
     316        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     317        self.failUnless(type(new_points) == type(points), 'failed')
     318        msg = ('First call of .get_absolute() returned %s\nexpected %s'
     319               % (str(new_points), str(expected_new_points)))
     320        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
     321
     322        # and repeat from 'new_points = g.get_absolute(points)' above
     323        # to see if second call with same input gives same results.
     324        new_points = g.get_absolute(points)
     325
     326        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     327        self.failUnless(type(new_points) == type(points), 'failed')
     328        msg = ('Second call of .get_absolute() returned %s\nexpected %s'
     329               % (str(new_points), str(expected_new_points)))
     330        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
     331
     332        # and repeat again to see if *third* call with same input
     333        # gives same results.
     334        new_points = g.get_absolute(points)
     335
     336        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     337        self.failUnless(type(new_points) == type(points), 'failed')
     338        msg = ('Second call of .get_absolute() returned %s\nexpected %s'
     339               % (str(new_points), str(expected_new_points)))
     340        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
     341
    262342    def test_is_absolute(self):
    263343       
    … …  
    427507if __name__ == "__main__":
    428508
    429     suite = unittest.makeSuite(geo_referenceTestCase,'test')
     509    #suite = unittest.makeSuite(geo_referenceTestCase,'test')
     510    suite = unittest.makeSuite(geo_referenceTestCase,'test_get_absolute')
    430511    runner = unittest.TextTestRunner() #verbosity=2)
    431512    runner.run(suite)

  • branches/numpy/anuga/coordinate_transforms/test_lat_long_UTM_conversion.py

    r6149 r6304  
    1212from anuga.utilities.anuga_exceptions import ANUGAError
    1313
    14 import Numeric as num
     14import numpy as num
    1515
    1616

  • branches/numpy/anuga/coordinate_transforms/test_redfearn.py

    r6149 r6304  
    1111from anuga.utilities.anuga_exceptions import ANUGAError
    1212
    13 import Numeric as num
     13import numpy as num
    1414
    1515

  • branches/numpy/anuga/culvert_flows/Test_Culvert_Flat_Water_Lev.py

    r6150 r6304  
    2626from math import pi,pow,sqrt
    2727
    28 import Numeric as num
     28import numpy as num
    2929
    3030

  • branches/numpy/anuga/culvert_flows/culvert_class.py

    r6150 r6304  
    1414from anuga.config import minimum_allowed_height, velocity_protection       
    1515
    16 import Numeric as num
     16import numpy as num
    1717
    1818

  • branches/numpy/anuga/culvert_flows/culvert_polygons.py

    r6150 r6304  
    66from anuga.utilities.polygon import inside_polygon, polygon_area
    77
    8 import Numeric as num
     8import numpy as num
    99
    1010

  • branches/numpy/anuga/culvert_flows/test_culvert_class.py

    r6150 r6304  
    2121from math import pi,pow,sqrt
    2222
    23 import Numeric as num
     23import numpy as num
    2424
    2525
    … …  
    470470        ref_volume = domain.get_quantity('stage').get_integral()
    471471        for t in domain.evolve(yieldstep = 1, finaltime = 25):
     472           
    472473            #print domain.timestepping_statistics()
    473474            new_volume = domain.get_quantity('stage').get_integral()
    474            
     475
    475476            msg = 'Total volume has changed'
    476477            assert num.allclose(new_volume, ref_volume), msg
    … …  
    564565#-------------------------------------------------------------
    565566if __name__ == "__main__":
    566     #suite = unittest.makeSuite(Test_Culvert, 'test_that_culvert_rating_limits_flow_in_shallow_inlet_condition')
    567567    suite = unittest.makeSuite(Test_Culvert, 'test')
    568     runner = unittest.TextTestRunner()
     568    runner = unittest.TextTestRunner() #verbosity=2)
    569569    runner.run(suite)
    570570       

  • branches/numpy/anuga/damage_modelling/inundation_damage.py

    r6151 r6304  
    1010from types import StringType
    1111
    12 import Numeric as num
     12import numpy as num
    1313
    1414
    … …  
    1717except ImportError: 
    1818    # Hand-built mockup of the things we need from the kinds package, since it
    19     # was recently removed from the standard Numeric distro.  Some users may 
     19    # was recently removed from the standard numeric distro.  Some users may 
    2020    # not have it by default. 
    2121    class _bunch: 
    … …  
    314314
    315315        # the data being created
    316         struct_damage = num.zeros(self.structure_count, num.Float)
    317         contents_damage = num.zeros(self.structure_count, num.Float)
     316        struct_damage = num.zeros(self.structure_count, num.float)
     317        contents_damage = num.zeros(self.structure_count, num.float)
    318318        self.struct_inundated = ['']* self.structure_count
    319319

  • branches/numpy/anuga/damage_modelling/test_inundation_damage.py

    r6151 r6304  
    1717from anuga.shallow_water.data_manager import get_dataobject
    1818
    19 import Numeric as num
     19import numpy as num
    2020
    2121
    … …  
    8383        #Initial condition - with jumps
    8484        bed = domain.quantities['elevation'].vertex_values
    85         stage = num.zeros(bed.shape, num.Float)
     85        stage = num.zeros(bed.shape, num.float)
    8686
    8787        h = 0.3
    … …  
    153153        #Initial condition - with jumps
    154154        bed = domain.quantities['elevation'].vertex_values
    155         stage = num.zeros(bed.shape, num.Float)
     155        stage = num.zeros(bed.shape, num.float)
    156156
    157157        h = 30.
    … …  
    476476        edm = EventDamageModel([0.0]*17, [0.0]*17, [0.0]*17,
    477477                               [0.0]*17, [0.0]*17)
    478         edm.struct_damage = num.zeros(17,num.Float)
    479         edm.contents_damage = num.zeros(17,num.Float)
     478        edm.struct_damage = num.zeros(17,num.float)
     479        edm.contents_damage = num.zeros(17,num.float)
    480480        collapse_probability = {0.4:[0], #0
    481481                                0.6:[1], #1
    … …  
    597597        pass
    598598    suite = unittest.makeSuite(Test_inundation_damage,'test')
    599     #suite = unittest.makeSuite(Test_inundation_damage,'test_inundation_damage_list_as_input')
    600599    runner = unittest.TextTestRunner()
    601600    runner.run(suite)

  • branches/numpy/anuga/fit_interpolate/fit.py

    r6244 r6304  
    4545class VertsWithNoTrianglesError(exceptions.Exception): pass
    4646
    47 import Numeric as num
     47import numpy as num
    4848
    4949
    … …  
    6666
    6767          vertex_coordinates: List of coordinate pairs [xi, eta] of
    68               points constituting a mesh (or an m x 2 Numeric array or
     68              points constituting a mesh (or an m x 2 numeric array or
    6969              a geospatial object)
    7070              Points may appear multiple times
    7171              (e.g. if vertices have discontinuities)
    7272
    73           triangles: List of 3-tuples (or a Numeric array) of
     73          triangles: List of 3-tuples (or a numeric array) of
    7474              integers representing indices of all vertices in the mesh.
    7575
    … …  
    251251            if len(z.shape) > 1:
    252252                att_num = z.shape[1]
    253                 self.Atz = num.zeros((m,att_num), num.Float)
     253                self.Atz = num.zeros((m,att_num), num.float)
    254254            else:
    255255                att_num = 1
    256                 self.Atz = num.zeros((m,), num.Float)
     256                self.Atz = num.zeros((m,), num.float)
    257257            assert z.shape[0] == point_coordinates.shape[0]
    258258
    … …  
    336336        point_coordinates: The co-ordinates of the data points.
    337337              List of coordinate pairs [x, y] of
    338               data points or an nx2 Numeric array or a Geospatial_data object
     338              data points or an nx2 numeric array or a Geospatial_data object
    339339              or points file filename
    340340          z: Single 1d vector or array of data at the point_coordinates.
    … …  
    382382        if point_coordinates is None:
    383383            if verbose: print 'Warning: no data points in fit'
    384             assert self.AtA <> None, 'no interpolation matrix'
    385             assert self.Atz <> None
     384            #assert self.AtA != None, 'no interpolation matrix'
     385            #assert self.Atz != None
     386            assert not self.AtA is None, 'no interpolation matrix'
     387            assert not self.Atz is None
    386388           
    387389            # FIXME (DSG) - do  a message
    … …  
    433435        point_coordinates: The co-ordinates of the data points.
    434436              List of coordinate pairs [x, y] of
    435               data points or an nx2 Numeric array or a Geospatial_data object
     437              data points or an nx2 numeric array or a Geospatial_data object
    436438        z: Single 1d vector or array of data at the point_coordinates.
    437439        attribute_name: Used to get the z values from the
    … …  
    448450                absolute = True)
    449451       
    450         # Convert input to Numeric arrays
     452        # Convert input to numeric arrays
    451453        if z is not None:
    452             z = ensure_numeric(z, num.Float)
     454            z = ensure_numeric(z, num.float)
    453455        else:
    454456            msg = 'z not specified'
    … …  
    456458            z = point_coordinates.get_attributes(attribute_name)
    457459
    458         point_coordinates = ensure_numeric(point_coordinates, num.Float)
     460        point_coordinates = ensure_numeric(point_coordinates, num.float)
    459461        self._build_matrix_AtA_Atz(point_coordinates, z, verbose)
    460462
    … …  
    556558        Inputs:
    557559        vertex_coordinates: List of coordinate pairs [xi, eta] of
    558               points constituting a mesh (or an m x 2 Numeric array or
     560              points constituting a mesh (or an m x 2 numeric array or
    559561              a geospatial object)
    560562              Points may appear multiple times
    561563              (e.g. if vertices have discontinuities)
    562564
    563           triangles: List of 3-tuples (or a Numeric array) of
     565          triangles: List of 3-tuples (or a numeric array) of
    564566          integers representing indices of all vertices in the mesh.
    565567
    566568          point_coordinates: List of coordinate pairs [x, y] of data points
    567           (or an nx2 Numeric array). This can also be a .csv/.txt/.pts
     569          (or an nx2 numeric array). This can also be a .csv/.txt/.pts
    568570          file name.
    569571
    … …  
    593595        # are None
    594596           
    595         #Convert input to Numeric arrays
    596         triangles = ensure_numeric(triangles, num.Int)
     597        #Convert input to numeric arrays
     598        triangles = ensure_numeric(triangles, num.int)
    597599        vertex_coordinates = ensure_absolute(vertex_coordinates,
    598600                                             geo_reference = mesh_origin)
    … …  
    662664    vertex_coordinates = mesh_dict['vertices']
    663665    triangles = mesh_dict['triangles']
    664     if type(mesh_dict['vertex_attributes']) == num.ArrayType:
     666    if isinstance(mesh_dict['vertex_attributes'], num.ndarray):
    665667        old_point_attributes = mesh_dict['vertex_attributes'].tolist()
    666668    else:
    667669        old_point_attributes = mesh_dict['vertex_attributes']
    668670
    669     if type(mesh_dict['vertex_attribute_titles']) == num.ArrayType:
     671    if isinstance(mesh_dict['vertex_attribute_titles'], num.ndarray):
    670672        old_title_list = mesh_dict['vertex_attribute_titles'].tolist()
    671673    else:

  • branches/numpy/anuga/fit_interpolate/general_fit_interpolate.py

    r6152 r6304  
    3535from anuga.fit_interpolate.search_functions import set_last_triangle
    3636
    37 import Numeric as num
     37import numpy as num
    3838
    3939
    … …  
    6767
    6868          vertex_coordinates: List of coordinate pairs [xi, eta] of
    69               points constituting a mesh (or an m x 2 Numeric array or
     69              points constituting a mesh (or an m x 2 numeric array or
    7070              a geospatial object)
    7171              Points may appear multiple times
    7272              (e.g. if vertices have discontinuities)
    7373
    74           triangles: List of 3-tuples (or a Numeric array) of
     74          triangles: List of 3-tuples (or a numeric array) of
    7575              integers representing indices of all vertices in the mesh.
    7676
    … …  
    9696                # are None
    9797           
    98                 #Convert input to Numeric arrays
    99                 triangles = ensure_numeric(triangles, num.Int)
     98                #Convert input to numeric arrays
     99                triangles = ensure_numeric(triangles, num.int)
    100100                vertex_coordinates = ensure_absolute(vertex_coordinates,
    101101                                                 geo_reference = mesh_origin)

  • branches/numpy/anuga/fit_interpolate/interpolate.py

    r6223 r6304  
    4141
    4242
    43 import Numeric as num
     43import numpy as num
    4444
    4545
    … …  
    7777    vertex_coordinates: List of coordinate pairs [xi, eta] of
    7878                        points constituting a mesh
    79                         (or an m x 2 Numeric array or
     79                        (or an m x 2 numeric array or
    8080                        a geospatial object)
    8181                        Points may appear multiple times
    8282                        (e.g. if vertices have discontinuities)
    8383
    84     triangles: List of 3-tuples (or a Numeric array) of
     84    triangles: List of 3-tuples (or a numeric array) of
    8585               integers representing indices of all vertices
    8686               in the mesh.
    … …  
    9292    interpolation_points: Interpolate mesh data to these positions.
    9393                          List of coordinate pairs [x, y] of
    94                           data points or an nx2 Numeric array or a
     94                          data points or an nx2 numeric array or a
    9595                          Geospatial_data object
    9696               
    … …  
    132132
    133133    # Create interpolation object with matrix
    134     args = (ensure_numeric(vertex_coordinates, num.Float),
     134    args = (ensure_numeric(vertex_coordinates, num.float),
    135135            ensure_numeric(triangles))
    136136    kwargs = {'mesh_origin': mesh_origin,
    … …  
    181181        Inputs:
    182182          vertex_coordinates: List of coordinate pairs [xi, eta] of
    183               points constituting a mesh (or an m x 2 Numeric array or
     183              points constituting a mesh (or an m x 2 numeric array or
    184184              a geospatial object)
    185185              Points may appear multiple times
    186186              (e.g. if vertices have discontinuities)
    187187
    188           triangles: List of 3-tuples (or a Numeric array) of
     188          triangles: List of 3-tuples (or a numeric array) of
    189189              integers representing indices of all vertices in the mesh.
    190190
    … …  
    243243          point_coordinates: Interpolate mesh data to these positions.
    244244              List of coordinate pairs [x, y] of
    245               data points or an nx2 Numeric array or a Geospatial_data object
     245              data points or an nx2 numeric array or a Geospatial_data object
    246246             
    247247              If point_coordinates is absent, the points inputted last time
    … …  
    294294                # creating a dummy array to concatenate to.
    295295               
    296                 f = ensure_numeric(f, num.Float)
     296                f = ensure_numeric(f, num.float)
    297297                if len(f.shape) > 1:
    298                     z = num.zeros((0, f.shape[1]))
     298                    z = num.zeros((0, f.shape[1]), num.int)     #array default#
    299299                else:
    300                     z = num.zeros((0,))
     300                    z = num.zeros((0,), num.int)        #array default#
    301301                   
    302302                for end in range(start_blocking_len,
    … …  
    340340            point_coordinates = point_coordinates.get_data_points(absolute=True)
    341341
    342         # Convert lists to Numeric arrays if necessary
    343         point_coordinates = ensure_numeric(point_coordinates, num.Float)
    344         f = ensure_numeric(f, num.Float)               
     342        # Convert lists to numeric arrays if necessary
     343        point_coordinates = ensure_numeric(point_coordinates, num.float)
     344        f = ensure_numeric(f, num.float)               
    345345
    346346        from anuga.caching import myhash
    … …  
    447447        if verbose: print 'Building interpolation matrix'
    448448
    449         # Convert point_coordinates to Numeric arrays, in case it was a list.
    450         point_coordinates = ensure_numeric(point_coordinates, num.Float)
     449        # Convert point_coordinates to numeric arrays, in case it was a list.
     450        point_coordinates = ensure_numeric(point_coordinates, num.float)
    451451
    452452        if verbose: print 'Getting indices inside mesh boundary'
    … …  
    526526    points: Interpolate mesh data to these positions.
    527527            List of coordinate pairs [x, y] of
    528             data points or an nx2 Numeric array or a Geospatial_data object
     528            data points or an nx2 numeric array or a Geospatial_data object
    529529             
    530530    No test for this yet.
    … …  
    694694
    695695    Mandatory input
    696         time:                 px1 array of monotonously increasing times (Float)
    697         quantities:           Dictionary of arrays or 1 array (Float)
     696        time:                 px1 array of monotonously increasing times (float)
     697        quantities:           Dictionary of arrays or 1 array (float)
    698698                              The arrays must either have dimensions pxm or mx1.
    699699                              The resulting function will be time dependent in
    … …  
    704704        quantity_names:       List of keys into the quantities dictionary for
    705705                              imposing a particular order on the output vector.
    706         vertex_coordinates:   mx2 array of coordinates (Float)
     706        vertex_coordinates:   mx2 array of coordinates (float)
    707707        triangles:            nx3 array of indices into vertex_coordinates (Int)
    708708        interpolation_points: Nx2 array of coordinates to be interpolated to
    … …  
    810810        self.quantities_range = {}
    811811        for name in quantity_names:
    812             q = quantities[name][:].flat
     812            q = quantities[name][:].flatten()
    813813            self.quantities_range[name] = [min(q), max(q)]
    814814       
    … …  
    830830                    interpolation_points = ensure_numeric(interpolation_points)
    831831            except:
    832                 msg = 'Interpolation points must be an N x 2 Numeric array ' \
     832                msg = 'Interpolation points must be an N x 2 numeric array ' \
    833833                      'or a list of points\n'
    834834                msg += 'Got: %s.' %(str(self.interpolation_points)[:60] + '...')
    … …  
    913913           
    914914            for name in quantity_names:
    915                 self.precomputed_values[name] = num.zeros((p, m), num.Float)
     915                self.precomputed_values[name] = num.zeros((p, m), num.float)
    916916
    917917            if verbose is True:
    … …  
    10561056
    10571057        # Compute interpolated values
    1058         q = num.zeros(len(self.quantity_names), num.Float)
     1058        q = num.zeros(len(self.quantity_names), num.float)
    10591059        for i, name in enumerate(self.quantity_names):
    10601060            Q = self.precomputed_values[name]
    … …  
    11051105                    res = []
    11061106                    for col in q:
    1107                         res.append(col*num.ones(N, num.Float))
     1107                        res.append(col*num.ones(N, num.float))
    11081108                       
    11091109                return res
    … …  
    11451145            minq, maxq = self.quantities_range[name]
    11461146            str += '    %s in [%f, %f]\n' %(name, minq, maxq)           
    1147             #q = quantities[name][:].flat
     1147            #q = quantities[name][:].flatten()
    11481148            #str += '    %s in [%f, %f]\n' %(name, min(q), max(q))
    11491149
    … …  
    11581158       
    11591159            for name in quantity_names:
    1160                 q = precomputed_values[name][:].flat
     1160                q = precomputed_values[name][:].flatten()
    11611161                str += '    %s at interpolation points in [%f, %f]\n'\
    11621162                       %(name, min(q), max(q))
    … …  
    11901190
    11911191    #Add the x and y together
    1192     vertex_coordinates = num.concatenate((x[:,num.NewAxis], y[:,num.NewAxis]),axis=1)
     1192    vertex_coordinates = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]),axis=1)
    11931193
    11941194    #Will return the quantity values at the specified times and locations
    … …  
    12621262    keys.remove("volumes")
    12631263    keys.remove("time")
    1264      #Turn NetCDF objects into Numeric arrays
     1264     #Turn NetCDF objects into numeric arrays
    12651265    quantities = {}
    12661266    for name in keys:

  • branches/numpy/anuga/fit_interpolate/search_functions.py

    r6174 r6304  
    1111from anuga.config import max_float
    1212
    13 import Numeric as num
     13import numpy as num
    1414
    1515
    … …  
    2222                        num.array([max_float,max_float]),
    2323                        num.array([max_float,max_float])),
    24                        (num.array([1,1], num.Int),      #array default#
    25                         num.array([0,0], num.Int),      #array default#
     24                       (num.array([1,1]),
     25                        num.array([0,0]),
    2626                        num.array([-1.1,-1.1]))]]]
    2727

  • branches/numpy/anuga/fit_interpolate/test_fit.py

    r6174 r6304  
    2121from anuga.shallow_water import Domain
    2222
    23 import Numeric as num
     23import numpy as num
    2424
    2525
    … …  
    582582
    583583        #Define f(x,y) = x
    584         f = num.array([0,0,2], num.Int) #Value at global vertex 2      #array default#
     584        f = num.array([0,0,2]) #Value at global vertex 2
    585585
    586586        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
    … …  
    589589
    590590        #Define f(x,y) = y
    591         f = num.array([0,2,0], num.Int)  #Value at global vertex 1      #array default#
     591        f = num.array([0,2,0])  #Value at global vertex 1
    592592
    593593        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
    … …  
    596596
    597597        #Define f(x,y) = x+y
    598         f = num.array([0,2,2], num.Int)  #Values at global vertex 1 and 2      #array default#
     598        f = num.array([0,2,2])  #Values at global vertex 1 and 2
    599599
    600600        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
    … …  
    623623
    624624        #Define f(x,y) = x
    625         f = num.array([0,0,2,0,2,4], num.Int) #f evaluated at points a-f      #array default#
     625        f = num.array([0,0,2,0,2,4]) #f evaluated at points a-f
    626626
    627627        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
    … …  
    630630
    631631        #Define f(x,y) = y
    632         f = num.array([0,2,0,4,2,0], num.Int) #f evaluated at points a-f      #array default#
     632        f = num.array([0,2,0,4,2,0]) #f evaluated at points a-f
    633633
    634634        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
    … …  
    637637
    638638        #Define f(x,y) = x+y
    639         f = num.array([0,2,2,4,4,4], num.Int)  #f evaluated at points a-f     #array default#
     639        f = num.array([0,2,2,4,4,4])  #f evaluated at points a-f
    640640
    641641        #Check that int (df/dx)**2 + (df/dy)**2 dx dy =
    … …  
    10991099if __name__ == "__main__":
    11001100    suite = unittest.makeSuite(Test_Fit,'test')
    1101     #suite = unittest.makeSuite(Test_Fit,'test_smooth_attributes_to_mesh_function')
    1102     #suite = unittest.makeSuite(Test_Fit,'')
    1103     runner = unittest.TextTestRunner(verbosity=1)
     1101    runner = unittest.TextTestRunner() #verbosity=1)
    11041102    runner.run(suite)
    11051103

  • branches/numpy/anuga/fit_interpolate/test_interpolate.py

    r6189 r6304  
    1515from Scientific.IO.NetCDF import NetCDFFile
    1616
    17 import Numeric as num
     17import numpy as num
    1818
    1919
    … …  
    6767
    6868        bed = domain.quantities['elevation'].vertex_values
    69         stage = num.zeros(bed.shape, num.Float)
     69        stage = num.zeros(bed.shape, num.float)
    7070
    7171        h = 0.3
    … …  
    129129
    130130        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
    131         vertex_coordinates = num.concatenate( (x[:, num.NewAxis], y[:, num.NewAxis]), axis=1 )
     131        vertex_coordinates = num.concatenate( (x[:, num.newaxis], y[:, num.newaxis]), axis=1 )
    132132        # FIXME: This concat should roll into get_vertex_values
    133133
    … …  
    149149       
    150150        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
    151         vertex_coordinates = num.concatenate( (x[:, num.NewAxis], y[:, num.NewAxis]), axis=1 )
     151        vertex_coordinates = num.concatenate( (x[:, num.newaxis], y[:, num.newaxis]), axis=1 )
    152152        # FIXME: This concat should roll into get_vertex_values
    153153
    … …  
    182182
    183183        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
    184         vertex_coordinates = num.concatenate( (x[:, num.NewAxis], y[:, num.NewAxis]), axis=1 )
     184        vertex_coordinates = num.concatenate( (x[:, num.newaxis], y[:, num.newaxis]), axis=1 )
    185185        # FIXME: This concat should roll into get_vertex_values
    186186
    … …  
    201201       
    202202        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
    203         vertex_coordinates = num.concatenate( (x[:, num.NewAxis], y[:, num.NewAxis]), axis=1 )
     203        vertex_coordinates = num.concatenate( (x[:, num.newaxis], y[:, num.newaxis]), axis=1 )
    204204        # FIXME: This concat should roll into get_vertex_values
    205205
    … …  
    233233       
    234234        x, y, vertex_values, triangles = quantity.get_vertex_values(xy=True, smooth=False)
    235         vertex_coordinates = num.concatenate( (x[:, num.NewAxis], y[:, num.NewAxis]), axis=1 )
     235        vertex_coordinates = num.concatenate( (x[:, num.newaxis], y[:, num.newaxis]), axis=1 )
    236236        # FIXME: This concat should roll into get_vertex_values
    237237
    … …  
    10781078
    10791079        #One quantity
    1080         Q = num.zeros( (3,6), num.Float )
     1080        Q = num.zeros( (3,6), num.float )
    10811081
    10821082        #Linear in time and space
    … …  
    12251225
    12261226        #One quantity
    1227         Q = num.zeros( (3,6), num.Float )
     1227        Q = num.zeros( (3,6), num.float )
    12281228
    12291229        #Linear in time and space
    … …  
    12851285
    12861286        # One quantity
    1287         Q = num.zeros((8,6), num.Float)
     1287        Q = num.zeros((8,6), num.float)
    12881288
    12891289        # Linear in time and space
    … …  
    13531353
    13541354        #One quantity
    1355         Q = num.zeros( (2,6), num.Float )
     1355        Q = num.zeros( (2,6), num.float )
    13561356
    13571357        #Linear in time and space
    … …  
    14191419
    14201420        # One quantity
    1421         Q = num.zeros( (3,6), num.Float )
     1421        Q = num.zeros( (3,6), num.float )
    14221422
    14231423        # Linear in time and space
    … …  
    16081608
    16091609        #One quantity
    1610         Q = num.zeros( (len(time),6), num.Float )
     1610        Q = num.zeros( (len(time),6), num.float )
    16111611
    16121612        #Linear in time and space
    … …  
    18501850if __name__ == "__main__":
    18511851    suite = unittest.makeSuite(Test_Interpolate,'test')
    1852     #suite = unittest.makeSuite(Test_Interpolate,'test_interpolate_one_point_many_triangles')
    1853     runner = unittest.TextTestRunner(verbosity=1)
     1852    runner = unittest.TextTestRunner() #verbosity=1)
    18541853    runner.run(suite)
    18551854

  • branches/numpy/anuga/fit_interpolate/test_search_functions.py

    r6152 r6304  
    211211if __name__ == "__main__":
    212212    suite = unittest.makeSuite(Test_search_functions,'test')
    213     #suite = unittest.makeSuite(Test_search_functions,'expanding_search')
    214     runner = unittest.TextTestRunner(verbosity=1)
     213    runner = unittest.TextTestRunner() #verbosity=1)
    215214    runner.run(suite)
    216215   

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

    r6166 r6304  
    1313from Scientific.IO.NetCDF import NetCDFFile
    1414
    15 import Numeric as num
     15import numpy as num
    1616
    1717from anuga.coordinate_transforms.lat_long_UTM_conversion import UTMtoLL
    … …  
    2323from anuga.config import points_file_block_line_size as MAX_READ_LINES
    2424from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
     25from anuga.config import netcdf_float
    2526
    2627DEFAULT_ATTRIBUTE = 'elevation'
    … …  
    5657                 load_file_now=True,
    5758                 verbose=False):
    58         """
    59         Create instance from data points and associated attributes
     59        """Create instance from data points and associated attributes
    6060
    6161        data_points: x,y coordinates in meters. Type must be either a
    62         sequence of 2-tuples or an Mx2 Numeric array of floats.  A file name
     62        sequence of 2-tuples or an Mx2 numeric array of floats.  A file name
    6363        with extension .txt, .cvs or .pts can also be passed in here.
    6464
    … …  
    131131
    132132        verbose:
    133 
    134133        """
    135134
    136135        if isinstance(data_points, basestring):
    137             # assume data point is really a file name
     136            # assume data_points is really a file name
    138137            file_name = data_points
    139138
    140139        self.set_verbose(verbose)
    141         self.geo_reference = None #create the attribute
     140        self.geo_reference = None # create the attribute
    142141        self.file_name = file_name
    143142
    … …  
    157156                                        data_points=data_points,
    158157                                        points_are_lats_longs=\
    159                                             points_are_lats_longs)
     158                                        points_are_lats_longs)
    160159            self.check_data_points(data_points)
    161160            self.set_attributes(attributes)
    … …  
    188187                        # This message was misleading.
    189188                        # FIXME (Ole): Are we blocking here or not?
    190                         #print 'ASCII formats are not that great for '
    191                         #print 'blockwise reading. Consider storing this'
    192                         #print 'data as a pts NetCDF format'
     189                        # print 'ASCII formats are not that great for '
     190                        # print 'blockwise reading. Consider storing this'
     191                        # print 'data as a pts NetCDF format'
    193192
    194193    ##
    … …  
    203202
    204203    ##
    205     # @brief Save points data in instance.
    206     # @param data_points Points data to store in instance and check.
     204    # @brief Check data points.
     205    # @param data_points Points data to check and store in instance.
    207206    # @note Throws ValueError exception if no data.
    208207    def check_data_points(self, data_points):
    209         """Checks data points
    210         """
     208        """Checks data points"""
    211209
    212210        if data_points is None:
    … …  
    217215        else:
    218216            self.data_points = ensure_numeric(data_points)
     217            return
     218
     219            print 'self.data_points=%s' % str(self.data_points)
     220            print 'self.data_points.shape=%s' % str(self.data_points.shape)
    219221            if not (0,) == self.data_points.shape:
    220222                assert len(self.data_points.shape) == 2
    … …  
    226228    # @note Throws exception if unable to convert dict keys to numeric.
    227229    def set_attributes(self, attributes):
    228         """Check and assign attributes dictionary
    229         """
     230        """Check and assign attributes dictionary"""
    230231
    231232        if attributes is None:
    … …  
    234235
    235236        if not isinstance(attributes, DictType):
    236             #Convert single attribute into dictionary
     237            # Convert single attribute into dictionary
    237238            attributes = {DEFAULT_ATTRIBUTE: attributes}
    238239
    239         #Check input attributes
     240        # Check input attributes
    240241        for key in attributes.keys():
    241242            try:
    242243                attributes[key] = ensure_numeric(attributes[key])
    243244            except:
    244                 msg = 'Attribute %s could not be converted' %key
    245                 msg += 'to a numeric vector'
    246                 raise msg
     245                msg = ("Attribute '%s' (%s) could not be converted to a"
     246                       "numeric vector" % (str(key), str(attributes[key])))
     247                raise Exception, msg
    247248
    248249        self.attributes = attributes
    249250
    250     #def set_geo_reference(self, geo_reference):
    251     #    # FIXME (Ole): Backwards compatibility - deprecate
    252     #    self.setgeo_reference(geo_reference)
    253 
    254251    ##
    255252    # @brief Set the georeference of geospatial data.
    256     # @param geo_reference The georeference data    to set.
     253    # @param geo_reference The georeference data to set.
    257254    # @note Will raise exception if param not instance of Geo_reference.
    258255    def set_geo_reference(self, geo_reference):
    … …  
    277274            msg = 'Argument geo_reference must be a valid Geo_reference '
    278275            msg += 'object or None.'
    279             raise msg
     276            raise Expection, msg
    280277
    281278        # If a geo_reference already exists, change the point data according to
    … …  
    511508                        raise Exception, msg
    512509        else:
    513             #other is None:
     510            # other is None:
    514511            new_points = self.get_data_points(absolute=True)
    515512            new_attributes = self.attributes
    … …  
    525522    # @return The new geospatial object.
    526523    def __radd__(self, other):
    527         """Handle cases like None + Geospatial_data(...)
    528         """
     524        """Handle cases like None + Geospatial_data(...)"""
    529525
    530526        return self + other
    … …  
    542538    def import_points_file(self, file_name, delimiter=None, verbose=False):
    543539        """ load an .txt, .csv or .pts file
     540
    544541        Note: will throw an IOError/SyntaxError if it can't load the file.
    545542        Catch these!
    … …  
    571568            except SyntaxError, e:
    572569                # This should only be if there is a format error
    573                 msg = 'Could not open file %s. \n' %file_name
     570                msg = 'Problem with format of file %s. \n' %file_name
    574571                msg += Error_message['IOError']
    575572                raise SyntaxError, msg
    … …  
    591588                           as_lat_long=False, isSouthHemisphere=True):
    592589
    593         """
    594         write a points file, file_name, as a text (.csv) or binary (.pts) file
     590        """write a points file as a text (.csv) or binary (.pts) file
     591
    595592        file_name is the file name, including the extension
    596593        The point_dict is defined at the top of this file.
    … …  
    659656        """
    660657
    661         #FIXME: add the geo_reference to this
     658        # FIXME: add the geo_reference to this
    662659        points = self.get_data_points()
    663660        sampled_points = num.take(points, indices)
    … …  
    679676    # @note Points in each result object are selected randomly.
    680677    def split(self, factor=0.5, seed_num=None, verbose=False):
    681         """Returns two
    682         geospatial_data object, first is the size of the 'factor'
     678        """Returns two geospatial_data object, first is the size of the 'factor'
    683679        smaller the original and the second is the remainder. The two
    684         new object are disjoin set of each other.
     680        new objects are disjoint sets of each other.
    685681
    686682        Points of the two new object have selected RANDOMLY.
    … …  
    707703        if verbose: print "make unique random number list and get indices"
    708704
    709         total=num.array(range(self_size), num.Int)     #array default#
     705        total=num.array(range(self_size))
    710706        total_list = total.tolist()
    711707
    … …  
    731727        random_num = random_num.tolist()
    732728
    733         #need to sort and reverse so the pop() works correctly
     729        # need to sort and reverse so the pop() works correctly
    734730        random_num.sort()
    735731        random_num.reverse()
    … …  
    748744            random_list.append(remainder_list.pop(i))
    749745            j += 1
    750             #prints progress
     746            # prints progress
    751747            if verbose and round(random_num_len/10*k) == j:
    752748                print '(%s/%s)' % (j, random_num_len)
    … …  
    779775        """
    780776        from Scientific.IO.NetCDF import NetCDFFile
    781         #FIXME - what to do if the file isn't there
     777        # FIXME - what to do if the file isn't there
    782778
    783779        # FIXME (Ole): Shouldn't this go into the constructor?
    … …  
    941937                        data_points,
    942938                        points_are_lats_longs):
    943     """
    944     if the points has lat long info, assume it is in (lat, long) order.
    945     """
     939    """If the points has lat long info, assume it is in (lat, long) order."""
    946940
    947941    if geo_reference is not None:
    … …  
    10551049                                                 header,
    10561050                                                 max_read_lines=1e30)
    1057                                     #If the file is bigger that this, block..
     1051                                    # If the file is bigger that this, block..
    10581052                                    # FIXME (Ole) What's up here?
    10591053    except ANUGAError:
    … …  
    10821076    """
    10831077
    1084     line = file_pointer.readline()
     1078    line = file_pointer.readline().strip()
    10851079    header = clean_line(line, delimiter)
    10861080
    … …  
    10951089# @param verbose True if this function is to be verbose.
    10961090# @note Will throw IndexError, SyntaxError exceptions.
    1097 def _read_csv_file_blocking(file_pointer, header,
     1091def _read_csv_file_blocking(file_pointer,
     1092                            header,
    10981093                            delimiter=CSV_DELIMITER,
    10991094                            max_read_lines=MAX_READ_LINES,
    … …  
    11501145        raise StopIteration
    11511146
    1152     pointlist = num.array(points, num.Float)
     1147    pointlist = num.array(points, num.float)
    11531148    for key in att_dict.keys():
    1154         att_dict[key] = num.array(att_dict[key], num.Float)
     1149        att_dict[key] = num.array(att_dict[key], num.float)
    11551150
    11561151    # Do stuff here so the info is in lat's and longs
    … …  
    11831178# @note Will throw IOError and AttributeError exceptions.
    11841179def _read_pts_file_header(fid, verbose=False):
    1185     """Read the geo_reference and number_of_points from a .pts file
    1186     """
     1180    """Read the geo_reference and number_of_points from a .pts file"""
    11871181
    11881182    keys = fid.variables.keys()
    … …  
    12121206# @return Tuple of (pointlist, attributes).
    12131207def _read_pts_file_blocking(fid, start_row, fin_row, keys):
    1214     """Read the body of a .csv file.
    1215     """
     1208    """Read the body of a .csv file."""
    12161209
    12171210    pointlist = num.array(fid.variables['points'][start_row:fin_row])
    … …  
    12391232
    12401233    WARNING: This function mangles the point_atts data structure
    1241     #F??ME: (DSG)This format has issues.
     1234    # F??ME: (DSG)This format has issues.
    12421235    # There can't be an attribute called points
    12431236    # consider format change
    … …  
    12641257    shape = write_data_points.shape[0]
    12651258    outfile.createDimension('number_of_points', shape)
    1266     outfile.createDimension('number_of_dimensions', 2) #This is 2d data
     1259    outfile.createDimension('number_of_dimensions', 2) # This is 2d data
    12671260
    12681261    # Variable definition
    1269     outfile.createVariable('points', num.Float, ('number_of_points',
    1270                                                  'number_of_dimensions'))
    1271 
    1272     #create variables
    1273     outfile.variables['points'][:] = write_data_points #.astype(Float32)
     1262    outfile.createVariable('points', netcdf_float, ('number_of_points',
     1263                                                    'number_of_dimensions'))
     1264
     1265    # create variables
     1266    outfile.variables['points'][:] = write_data_points
    12741267
    12751268    if write_attributes is not None:
    12761269        for key in write_attributes.keys():
    1277             outfile.createVariable(key, num.Float, ('number_of_points',))
    1278             outfile.variables[key][:] = write_attributes[key] #.astype(Float32)
     1270            outfile.createVariable(key, netcdf_float, ('number_of_points',))
     1271            outfile.variables[key][:] = write_attributes[key]
    12791272
    12801273    if write_geo_reference is not None:
    … …  
    12961289                    as_lat_long=False,
    12971290                    delimiter=','):
    1298     """Write a .csv file.
    1299     """
     1291    """Write a .csv file."""
    13001292
    13011293    points = write_data_points
    … …  
    13611353# @return ??
    13621354def _point_atts2array(point_atts):
    1363     point_atts['pointlist'] = num.array(point_atts['pointlist'], num.Float)
     1355    point_atts['pointlist'] = num.array(point_atts['pointlist'], num.float)
    13641356
    13651357    for key in point_atts['attributelist'].keys():
    13661358        point_atts['attributelist'][key] = \
    1367                 num.array(point_atts['attributelist'][key], num.Float)
     1359                num.array(point_atts['attributelist'][key], num.float)
    13681360
    13691361    return point_atts
    … …  
    13751367# @return A points dictionary.
    13761368def geospatial_data2points_dictionary(geospatial_data):
    1377     """Convert geospatial data to points_dictionary
    1378     """
     1369    """Convert geospatial data to points_dictionary"""
    13791370
    13801371    points_dictionary = {}
    … …  
    13961387# @param points_dictionary A points dictionary to convert.
    13971388def points_dictionary2geospatial_data(points_dictionary):
    1398     """Convert points_dictionary to geospatial data object
    1399     """
     1389    """Convert points_dictionary to geospatial data object"""
    14001390
    14011391    msg = 'Points dictionary must have key pointlist'
    … …  
    14171407##
    14181408# @brief Split a string into 'clean' fields.
    1419 # @param line The string to process.
     1409# @param str The string to process.
    14201410# @param delimiter The delimiter string to split 'line' with.
    14211411# @return A list of 'cleaned' field strings.
    14221412# @note Any fields that were initially zero length will be removed.
    1423 def clean_line(line,delimiter):
    1424     """Remove whitespace
    1425     """
    1426 
    1427     line = line.strip()                 # probably unnecessary RW
    1428     numbers = line.split(delimiter)
    1429 
    1430     i = len(numbers) - 1
    1431     while i >= 0:
    1432         if numbers[i] == '':
    1433             numbers.pop(i)
    1434         else:
    1435             numbers[i] = numbers[i].strip()
    1436         i += -1
    1437 
    1438     return numbers
     1413# @note If a field contains '\n' it isn't zero length.
     1414def clean_line(str, delimiter):
     1415    """Split string on given delimiter, remove whitespace from each field."""
     1416
     1417    return [x.strip() for x in str.split(delimiter) if x != '']
    14391418
    14401419
    … …  
    14621441    # Input check
    14631442    if isinstance(points, basestring):
    1464         #It's a string - assume it is a point file
     1443        # It's a string - assume it is a point file
    14651444        points = Geospatial_data(file_name=points)
    14661445
    … …  
    14701449        assert geo_reference == None, msg
    14711450    else:
    1472         points = ensure_numeric(points, num.Float)
     1451        points = ensure_numeric(points, num.float)
    14731452
    14741453    # Sort of geo_reference and convert points
    14751454    if geo_reference is None:
    1476         geo = None #Geo_reference()
     1455        geo = None # Geo_reference()
    14771456    else:
    14781457        if isinstance(geo_reference, Geo_reference):
    … …  
    14931472# @return A geospatial object.
    14941473def ensure_geospatial(points, geo_reference=None):
    1495     """
    1496     This function inputs several formats and
    1497     outputs one format. - a geospatial_data instance.
     1474    """Convert various data formats to a geospatial_data instance.
    14981475
    14991476    Inputed formats are;
    … …  
    15141491    else:
    15151492        # List or numeric array of absolute points
    1516         points = ensure_numeric(points, num.Float)
     1493        points = ensure_numeric(points, num.float)
    15171494
    15181495    # Sort out geo reference
    … …  
    15631540                                     cache=False,
    15641541                                     verbose=False):
    1565     """
    1566     Removes a small random sample of points from 'data_file'. Then creates
    1567     models with different alpha values from 'alpha_list' and cross validates
    1568     the predicted value to the previously removed point data. Returns the
    1569     alpha value which has the smallest covariance.
     1542    """Removes a small random sample of points from 'data_file'.
     1543    Then creates models with different alpha values from 'alpha_list' and
     1544    cross validates the predicted value to the previously removed point data.
     1545    Returns the alpha value which has the smallest covariance.
    15701546
    15711547    data_file: must not contain points outside the boundaries defined
    … …  
    16281604        if no_boundary is True:
    16291605            msg = 'All boundaries must be defined'
    1630             raise msg
     1606            raise Expection, msg
    16311607
    16321608        poly_topo = [[east_boundary,south_boundary],
    … …  
    16451621
    16461622    else: # if mesh file provided
    1647         #test mesh file exists?
     1623        # test mesh file exists?
    16481624        if verbose: "reading from file: %s" % mesh_file
    16491625        if access(mesh_file,F_OK) == 0:
    … …  
    16511627            raise IOError, msg
    16521628
    1653     #split topo data
     1629    # split topo data
    16541630    if verbose: print 'Reading elevation file: %s' % data_file
    16551631    G = Geospatial_data(file_name = data_file)
    … …  
    16681644        alphas = alpha_list
    16691645
    1670     #creates array with columns 1 and 2 are x, y. column 3 is elevation
    1671     #4 onwards is the elevation_predicted using the alpha, which will
    1672     #be compared later against the real removed data
    1673     data = num.array([], typecode=num.Float)
     1646    # creates array with columns 1 and 2 are x, y. column 3 is elevation
     1647    # 4 onwards is the elevation_predicted using the alpha, which will
     1648    # be compared later against the real removed data
     1649    data = num.array([], dtype=num.float)
    16741650
    16751651    data=num.resize(data, (len(points), 3+len(alphas)))
    16761652
    1677     #gets relative point from sample
     1653    # gets relative point from sample
    16781654    data[:,0] = points[:,0]
    16791655    data[:,1] = points[:,1]
    … …  
    16811657    data[:,2] = elevation_sample
    16821658
    1683     normal_cov=num.array(num.zeros([len(alphas), 2]), typecode=num.Float)
     1659    normal_cov=num.array(num.zeros([len(alphas), 2]), dtype=num.float)
    16841660
    16851661    if verbose: print 'Setup computational domains with different alphas'
    16861662
    16871663    for i, alpha in enumerate(alphas):
    1688         #add G_other data to domains with different alphas
     1664        # add G_other data to domains with different alphas
    16891665        if verbose:
    16901666            print '\n Calculating domain and mesh for Alpha = ', alpha, '\n'
    … …  
    17001676        points_geo = Geospatial_data(points, domain.geo_reference)
    17011677
    1702         #returns the predicted elevation of the points that were "split" out
    1703         #of the original data set for one particular alpha
     1678        # returns the predicted elevation of the points that were "split" out
     1679        # of the original data set for one particular alpha
    17041680        if verbose: print 'Get predicted elevation for location to be compared'
    17051681        elevation_predicted = \
    … …  
    17071683                    get_values(interpolation_points=points_geo)
    17081684
    1709         #add predicted elevation to array that starts with x, y, z...
     1685        # add predicted elevation to array that starts with x, y, z...
    17101686        data[:,i+3] = elevation_predicted
    17111687
    … …  
    18341810        if no_boundary is True:
    18351811            msg = 'All boundaries must be defined'
    1836             raise msg
     1812            raise Expection, msg
    18371813
    18381814        poly_topo = [[east_boundary,south_boundary],
    … …  
    18511827
    18521828    else: # if mesh file provided
    1853         #test mesh file exists?
     1829        # test mesh file exists?
    18541830        if access(mesh_file,F_OK) == 0:
    18551831            msg = "file %s doesn't exist!" % mesh_file
    18561832            raise IOError, msg
    18571833
    1858     #split topo data
     1834    # split topo data
    18591835    G = Geospatial_data(file_name=data_file)
    18601836    if verbose: print 'start split'
    … …  
    18631839    points = G_small.get_data_points()
    18641840
    1865     #FIXME: Remove points outside boundary polygon
     1841    # FIXME: Remove points outside boundary polygon
    18661842#    print 'new point',len(points)
    18671843#
    18681844#    new_points=[]
    1869 #    new_points=array([],typecode=Float)
     1845#    new_points=array([],dtype=float)
    18701846#    new_points=resize(new_points,(len(points),2))
    18711847#    print "BOUNDARY", boundary_poly
    … …  
    18901866
    18911867    for alpha in alphas:
    1892         #add G_other data to domains with different alphas
     1868        # add G_other data to domains with different alphas
    18931869        if verbose:
    18941870            print '\n Calculating domain and mesh for Alpha = ', alpha, '\n'
    … …  
    19021878        domains[alpha] = domain
    19031879
    1904     #creates array with columns 1 and 2 are x, y. column 3 is elevation
    1905     #4 onwards is the elevation_predicted using the alpha, which will
    1906     #be compared later against the real removed data
    1907     data = num.array([], typecode=num.Float)
     1880    # creates array with columns 1 and 2 are x, y. column 3 is elevation
     1881    # 4 onwards is the elevation_predicted using the alpha, which will
     1882    # be compared later against the real removed data
     1883    data = num.array([], dtype=num.float)
    19081884
    19091885    data = num.resize(data, (len(points), 3+len(alphas)))
    19101886
    1911     #gets relative point from sample
     1887    # gets relative point from sample
    19121888    data[:,0] = points[:,0]
    19131889    data[:,1] = points[:,1]
    … …  
    19151891    data[:,2] = elevation_sample
    19161892
    1917     normal_cov = num.array(num.zeros([len(alphas), 2]), typecode=num.Float)
     1893    normal_cov = num.array(num.zeros([len(alphas), 2]), dtype=num.float)
    19181894
    19191895    if verbose:
    … …  
    19231899
    19241900        points_geo = domains[alpha].geo_reference.change_points_geo_ref(points)
    1925         #returns the predicted elevation of the points that were "split" out
    1926         #of the original data set for one particular alpha
     1901        # returns the predicted elevation of the points that were "split" out
     1902        # of the original data set for one particular alpha
    19271903        elevation_predicted = \
    19281904                domains[alpha].quantities[attribute_smoothed].\
    19291905                        get_values(interpolation_points=points_geo)
    19301906
    1931         #add predicted elevation to array that starts with x, y, z...
     1907        # add predicted elevation to array that starts with x, y, z...
    19321908        data[:,i+3] = elevation_predicted
    19331909

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

    r6153 r6304  
    11#!/usr/bin/env python
    2 
    32
    43import unittest
    54import os
     5import tempfile
    66from math import sqrt, pi
    7 import tempfile
    87from sets import ImmutableSet
    98
    10 import Numeric as num
     9import numpy as num
    1110
    1211from anuga.geospatial_data.geospatial_data import *
    … …  
    1615from anuga.utilities.system_tools import get_host_name
    1716from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
     17from anuga.config import netcdf_float
     18
    1819
    1920class Test_Geospatial_data(unittest.TestCase):
    … …  
    2425        pass
    2526
    26 
    2727    def test_0(self):
    2828        #Basic points
    2929        from anuga.coordinate_transforms.geo_reference import Geo_reference
    30        
     30
    3131        points = [[1.0, 2.1], [3.0, 5.3]]
    3232        G = Geospatial_data(points)
    33 
    3433        assert num.allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]])
    3534
    … …  
    3837        rep = `G`
    3938        ref = '[[ 1.   2.1]\n [ 3.   5.3]]'
    40 
    41         msg = 'Representation %s is not equal to %s' %(rep, ref)
     39        msg = 'Representation %s is not equal to %s' % (rep, ref)
    4240        assert rep == ref, msg
    4341
    4442        #Check getter
    4543        assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3]])
    46        
     44
    4745        #Check defaults
    4846        assert G.attributes is None
    49        
    5047        assert G.geo_reference.zone == Geo_reference().zone
    5148        assert G.geo_reference.xllcorner == Geo_reference().xllcorner
    5249        assert G.geo_reference.yllcorner == Geo_reference().yllcorner
    53        
    5450
    5551    def test_1(self):
    5652        points = [[1.0, 2.1], [3.0, 5.3]]
    5753        attributes = [2, 4]
    58         G = Geospatial_data(points, attributes)       
     54        G = Geospatial_data(points, attributes)
    5955        assert G.attributes.keys()[0] == DEFAULT_ATTRIBUTE
    6056        assert num.allclose(G.attributes.values()[0], [2, 4])
    61        
    6257
    6358    def test_2(self):
    6459        from anuga.coordinate_transforms.geo_reference import Geo_reference
     60
    6561        points = [[1.0, 2.1], [3.0, 5.3]]
    6662        attributes = [2, 4]
    … …  
    7268        assert G.geo_reference.yllcorner == 200
    7369
    74 
    7570    def test_get_attributes_1(self):
    7671        from anuga.coordinate_transforms.geo_reference import Geo_reference
     72
    7773        points = [[1.0, 2.1], [3.0, 5.3]]
    7874        attributes = [2, 4]
    … …  
    8076                            geo_reference=Geo_reference(56, 100, 200))
    8177
    82 
    8378        P = G.get_data_points(absolute=False)
    84         assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
     79        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])
    8580
    8681        P = G.get_data_points(absolute=True)
    87         assert num.allclose(P, [[101.0, 202.1], [103.0, 205.3]])       
     82        assert num.allclose(P, [[101.0, 202.1], [103.0, 205.3]])
    8883
    8984        V = G.get_attributes() #Simply get them
    … …  
    9590    def test_get_attributes_2(self):
    9691        #Multiple attributes
    97        
    98        
    9992        from anuga.coordinate_transforms.geo_reference import Geo_reference
     93
    10094        points = [[1.0, 2.1], [3.0, 5.3]]
    10195        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
    … …  
    10498                            default_attribute_name='a1')
    10599
    106 
    107100        P = G.get_data_points(absolute=False)
    108         assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
    109        
     101        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])
     102
    110103        V = G.get_attributes() #Get default attribute
    111104        assert num.allclose(V, [2, 4])
    … …  
    125118            pass
    126119        else:
    127             raise 'Should have raised exception'
     120            raise Exception, 'Should have raised exception'
    128121
    129122    def test_get_data_points(self):
    130         points_ab = [[12.5,34.7],[-4.5,-60.0]]
     123        points_ab = [[12.5, 34.7], [-4.5, -60.0]]
    131124        x_p = -10
    132125        y_p = -40
    133126        geo_ref = Geo_reference(56, x_p, y_p)
    134127        points_rel = geo_ref.change_points_geo_ref(points_ab)
    135        
     128
    136129        spatial = Geospatial_data(points_rel, geo_reference=geo_ref)
    137 
    138130        results = spatial.get_data_points(absolute=False)
    139        
    140131        assert num.allclose(results, points_rel)
    141        
     132
    142133        x_p = -1770
    143134        y_p = 4.01
    144135        geo_ref = Geo_reference(56, x_p, y_p)
    145136        points_rel = geo_ref.change_points_geo_ref(points_ab)
    146         results = spatial.get_data_points \
    147                   ( geo_reference=geo_ref)
    148        
     137        results = spatial.get_data_points(geo_reference=geo_ref)
     138
    149139        assert num.allclose(results, points_rel)
    150140
    151  
    152141    def test_get_data_points_lat_long(self):
    153         # lat long [-30.],[130]
    154         #Zone:   52   
    155         #Easting:  596450.153  Northing: 6680793.777
    156         # lat long [-32.],[131]
    157         #Zone:   52   
    158         #Easting:  688927.638  Northing: 6457816.509
    159        
    160         points_Lat_long = [[-30.,130], [-32,131]]
    161        
    162         spatial = Geospatial_data(latitudes=[-30, -32.],
    163                                   longitudes=[130, 131])
    164 
     142        # lat long [-30.], [130]
     143        # Zone:   52
     144        # Easting:  596450.153  Northing: 6680793.777
     145        # lat long [-32.], [131]
     146        # Zone:   52
     147        # Easting:  688927.638  Northing: 6457816.509
     148
     149        points_Lat_long = [[-30., 130], [-32, 131]]
     150
     151        spatial = Geospatial_data(latitudes=[-30, -32.], longitudes=[130, 131])
    165152        results = spatial.get_data_points(as_lat_long=True)
    166         #print "test_get_data_points_lat_long - results", results
    167         #print "points_Lat_long",points_Lat_long
    168153        assert num.allclose(results, points_Lat_long)
    169      
     154
    170155    def test_get_data_points_lat_longII(self):
    171156        # x,y  North,east long,lat
    172157        boundary_polygon = [[ 250000, 7630000]]
    173158        zone = 50
    174        
     159
    175160        geo_reference = Geo_reference(zone=zone)
    176         geo = Geospatial_data(boundary_polygon,geo_reference=geo_reference)
     161        geo = Geospatial_data(boundary_polygon ,geo_reference=geo_reference)
    177162        seg_lat_long = geo.get_data_points(as_lat_long=True)
    178         lat_result = degminsec2decimal_degrees(-21,24,54)
    179         long_result = degminsec2decimal_degrees(114,35,17.89)
    180         #print "seg_lat_long", seg_lat_long [0][0]
    181         #print "lat_result",lat_result
     163        lat_result = degminsec2decimal_degrees(-21, 24, 54)
     164        long_result = degminsec2decimal_degrees(114, 35, 17.89)
     165        assert num.allclose(seg_lat_long[0][0], lat_result)     #lat
     166        assert num.allclose(seg_lat_long[0][1], long_result)    #long
     167
     168    def test_get_data_points_lat_longIII(self):
     169        # x,y  North,east long,lat
     170        # for northern hemisphere
     171        boundary_polygon = [[419944.8, 918642.4]]
     172        zone = 47
     173
     174        geo_reference = Geo_reference(zone=zone)
     175        geo = Geospatial_data(boundary_polygon, geo_reference=geo_reference)
     176        seg_lat_long = geo.get_data_points(as_lat_long=True,
     177                                           isSouthHemisphere=False)
     178        lat_result = degminsec2decimal_degrees(8.31, 0, 0)
     179        long_result = degminsec2decimal_degrees(98.273, 0, 0)
    182180        assert num.allclose(seg_lat_long[0][0], lat_result)#lat
    183181        assert num.allclose(seg_lat_long[0][1], long_result)#long
    184182
    185 
    186     def test_get_data_points_lat_longIII(self):
    187         # x,y  North,east long,lat
    188         #for northern hemisphere
    189         boundary_polygon = [[419944.8, 918642.4]]
    190         zone = 47
    191        
    192         geo_reference = Geo_reference(zone=zone)
    193         geo = Geospatial_data(boundary_polygon,
    194                               geo_reference=geo_reference)
    195                              
    196         seg_lat_long = geo.get_data_points(as_lat_long=True,
    197                                            isSouthHemisphere=False)
    198                                            
    199         lat_result = degminsec2decimal_degrees(8.31,0,0)
    200         long_result = degminsec2decimal_degrees(98.273,0,0)
    201         #print "seg_lat_long", seg_lat_long [0]
    202         #print "lat_result",lat_result
    203         assert num.allclose(seg_lat_long[0][0], lat_result)#lat
    204         assert num.allclose(seg_lat_long[0][1], long_result)#long
    205 
    206 
    207              
    208183    def test_set_geo_reference(self):
    209         """test_set_georeference
    210        
    211         Test that georeference can be changed without changing the 
     184        '''test_set_georeference
     185
     186        Test that georeference can be changed without changing the
    212187        absolute values.
    213         """
    214            
    215         points_ab = [[12.5,34.7],[-4.5,-60.0]]
     188        '''
     189
     190        points_ab = [[12.5, 34.7], [-4.5, -60.0]]
    216191        x_p = -10
    217192        y_p = -40
    218193        geo_ref = Geo_reference(56, x_p, y_p)
    219194        points_rel = geo_ref.change_points_geo_ref(points_ab)
    220        
     195
    221196        # Create without geo_ref properly set
    222         G = Geospatial_data(points_rel)       
     197        G = Geospatial_data(points_rel)
    223198        assert not num.allclose(points_ab, G.get_data_points(absolute=True))
    224        
     199
    225200        # Create the way it should be
    226201        G = Geospatial_data(points_rel, geo_reference=geo_ref)
    227202        assert num.allclose(points_ab, G.get_data_points(absolute=True))
    228        
     203
    229204        # Change georeference and check that absolute values are unchanged.
    230205        x_p = 10
    … …  
    232207        new_geo_ref = Geo_reference(56, x_p, y_p)
    233208        G.set_geo_reference(new_geo_ref)
    234         assert num.allclose(points_ab, G.get_data_points(absolute=True))
    235        
    236 
    237                
    238        
     209        msg = ('points_ab=%s, but\nG.get_data_points(absolute=True)=%s'
     210               % (str(points_ab), str(G.get_data_points(absolute=True))))
     211        assert num.allclose(points_ab, G.get_data_points(absolute=True)), msg
     212
    239213    def test_conversions_to_points_dict(self):
    240214        #test conversions to points_dict
    241        
    242        
    243215        from anuga.coordinate_transforms.geo_reference import Geo_reference
     216
    244217        points = [[1.0, 2.1], [3.0, 5.3]]
    245218        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
    … …  
    248221                            default_attribute_name='a1')
    249222
    250 
    251223        points_dict = geospatial_data2points_dictionary(G)
    252224
    253225        assert points_dict.has_key('pointlist')
    254         assert points_dict.has_key('attributelist')       
     226        assert points_dict.has_key('attributelist')
    255227        assert points_dict.has_key('geo_reference')
    256228
    … …  
    260232        assert A.has_key('a0')
    261233        assert A.has_key('a1')
    262         assert A.has_key('a2')       
     234        assert A.has_key('a2')
    263235
    264236        assert num.allclose( A['a0'], [0, 0] )
    265         assert num.allclose( A['a1'], [2, 4] )       
     237        assert num.allclose( A['a1'], [2, 4] )
    266238        assert num.allclose( A['a2'], [79.4, -7] )
    267 
    268239
    269240        geo = points_dict['geo_reference']
    270241        assert geo is G.geo_reference
    271242
    272 
    273243    def test_conversions_from_points_dict(self):
    274         """test conversions from points_dict
    275         """
     244        '''test conversions from points_dict'''
    276245
    277246        from anuga.coordinate_transforms.geo_reference import Geo_reference
    278        
     247
    279248        points = [[1.0, 2.1], [3.0, 5.3]]
    280249        attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]}
    … …  
    284253        points_dict['attributelist'] = attributes
    285254        points_dict['geo_reference'] = Geo_reference(56, 100, 200)
    286        
    287255
    288256        G = points_dictionary2geospatial_data(points_dict)
    289 
    290257        P = G.get_data_points(absolute=False)
    291         assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])       
    292        
    293         #V = G.get_attribute_values() #Get default attribute
    294         #assert allclose(V, [2, 4])
     258        assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]])
    295259
    296260        V = G.get_attributes('a0') #Get by name
    … …  
    304268
    305269    def test_add(self):
    306         """ test the addition of two geospatical objects
    307             no geo_reference see next test
    308         """
     270        '''test the addition of two geospatical objects
     271        no geo_reference see next test
     272        '''
     273
    309274        points = [[1.0, 2.1], [3.0, 5.3]]
    310         attributes = {'depth':[2, 4], 'elevation':[6.1, 5]}
    311         attributes1 = {'depth':[2, 4], 'elevation':[2.5, 1]}
    312         G1 = Geospatial_data(points, attributes)       
    313         G2 = Geospatial_data(points, attributes1)
    314        
    315 #        g3 = geospatial_data2points_dictionary(G1)
    316 #        print 'g3=', g3
    317        
     275        attributes = {'depth': [2, 4], 'elevation': [6.1, 5]}
     276        attributes1 = {'depth': [2, 4], 'elevation': [2.5, 1]}
     277        G1 = Geospatial_data(points, attributes)
     278        G2 = Geospatial_data(points, attributes1)
     279
    318280        G = G1 + G2
    319281
    … …  
    326288
    327289    def test_addII(self):
    328         """ test the addition of two geospatical objects
    329             no geo_reference see next test
    330         """
     290        '''test the addition of two geospatical objects
     291        no geo_reference see next test
     292        '''
     293
    331294        points = [[1.0, 2.1], [3.0, 5.3]]
    332         attributes = {'depth':[2, 4]}
    333         G1 = Geospatial_data(points, attributes) 
    334        
     295        attributes = {'depth': [2, 4]}
     296        G1 = Geospatial_data(points, attributes)
     297
    335298        points = [[5.0, 2.1], [3.0, 50.3]]
    336         attributes = {'depth':[200, 400]}
     299        attributes = {'depth': [200, 400]}
    337300        G2 = Geospatial_data(points, attributes)
    338        
    339 #        g3 = geospatial_data2points_dictionary(G1)
    340 #        print 'g3=', g3
    341        
     301
    342302        G = G1 + G2
    343303
    344         assert G.attributes.has_key('depth') 
     304        assert G.attributes.has_key('depth')
    345305        assert G.attributes.keys(), ['depth']
    346306        assert num.allclose(G.attributes['depth'], [2, 4, 200, 400])
    347307        assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3],
    348308                                                  [5.0, 2.1], [3.0, 50.3]])
     309
    349310    def test_add_with_geo (self):
    350         """
    351         Difference in Geo_reference resolved
    352         """
     311        '''Difference in Geo_reference resolved'''
     312
    353313        points1 = [[1.0, 2.1], [3.0, 5.3]]
    354314        points2 = [[5.0, 6.1], [6.0, 3.3]]
    … …  
    362322                                 projection='UTM',
    363323                                 units='m')
    364                                
     324
    365325        G1 = Geospatial_data(points1, attributes1, geo_ref1)
    366326        G2 = Geospatial_data(points2, attributes2, geo_ref2)
    … …  
    371331
    372332        P2 = G2.get_data_points(absolute=True)
    373         assert num.allclose(P2, [[5.1, 9.1], [6.1, 6.3]])       
    374        
     333        assert num.allclose(P2, [[5.1, 9.1], [6.1, 6.3]])
     334
    375335        G = G1 + G2
    376336
    … …  
    381341        P = G.get_data_points(absolute=True)
    382342
    383         #P_relative = G.get_data_points(absolute=False)
    384         #
    385         #assert allclose(P_relative, P - [0.1, 2.0])
    386 
    387343        assert num.allclose(P, num.concatenate( (P1,P2) ))
     344        msg = 'P=%s' % str(P)
    388345        assert num.allclose(P, [[2.0, 4.1], [4.0, 7.3],
    389                                 [5.1, 9.1], [6.1, 6.3]])
    390        
    391 
    392 
    393        
     346                                [5.1, 9.1], [6.1, 6.3]]), msg
    394347
    395348    def test_add_with_geo_absolute (self):
    396         """
    397         Difference in Geo_reference resolved
    398         """
     349        '''Difference in Geo_reference resolved'''
     350
    399351        points1 = num.array([[2.0, 4.1], [4.0, 7.3]])
    400         points2 = num.array([[5.1, 9.1], [6.1, 6.3]])       
     352        points2 = num.array([[5.1, 9.1], [6.1, 6.3]])
    401353        attributes1 = [2, 4]
    402354        attributes2 = [5, 76]
    … …  
    404356        geo_ref2 = Geo_reference(55, 2.0, 3.0)
    405357
    406        
    407                                
    408         G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()],
     358        G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(),
     359                                        geo_ref1.get_yllcorner()],
    409360                             attributes1, geo_ref1)
    410        
    411         G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()],
     361
     362        G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(),
     363                                        geo_ref2.get_yllcorner()],
    412364                             attributes2, geo_ref2)
    413365
    … …  
    417369
    418370        P1 = G1.get_data_points(absolute=False)
    419         assert num.allclose(P1, points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()])       
     371        msg = ('P1=%s, but\npoints1 - <...>=%s'
     372               % (str(P1),
     373                  str(points1 - [geo_ref1.get_xllcorner(),
     374                                 geo_ref1.get_yllcorner()])))
     375        assert num.allclose(P1, points1 - [geo_ref1.get_xllcorner(),
     376                                           geo_ref1.get_yllcorner()]), msg
    420377
    421378        P2 = G2.get_data_points(absolute=True)
    … …  
    423380
    424381        P2 = G2.get_data_points(absolute=False)
    425         assert num.allclose(P2, points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()])               
    426        
     382        assert num.allclose(P2, points2 - [geo_ref2.get_xllcorner(),
     383                                           geo_ref2.get_yllcorner()])
     384
    427385        G = G1 + G2
    428        
    429         #assert allclose(G.get_geo_reference().get_xllcorner(), 1.0)
    430         #assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
    431 
    432386        P = G.get_data_points(absolute=True)
    433387
    434         #P_relative = G.get_data_points(absolute=False)
    435         #
    436         #assert allclose(P_relative, [[1.0, 2.1], [3.0, 5.3], [4.1, 7.1], [5.1, 4.3]])
    437 
    438         assert num.allclose(P, num.concatenate( (points1,points2) ))
    439 
     388        assert num.allclose(P, num.concatenate((points1, points2)))
    440389
    441390    def test_add_with_None(self):
    442         """ test that None can be added to a geospatical objects
    443         """
    444        
     391        '''test that None can be added to a geospatical objects'''
     392
    445393        points1 = num.array([[2.0, 4.1], [4.0, 7.3]])
    446         points2 = num.array([[5.1, 9.1], [6.1, 6.3]])       
     394        points2 = num.array([[5.1, 9.1], [6.1, 6.3]])
    447395
    448396        geo_ref1= Geo_reference(55, 1.0, 2.0)
    … …  
    453401                                 projection='UTM',
    454402                                 units='m')
    455        
    456 
    457         attributes1 = {'depth':[2, 4.7], 'elevation':[6.1, 5]}
    458         attributes2 = {'depth':[-2.3, 4], 'elevation':[2.5, 1]}
    459 
     403
     404        attributes1 = {'depth': [2, 4.7], 'elevation': [6.1, 5]}
     405        attributes2 = {'depth': [-2.3, 4], 'elevation': [2.5, 1]}
    460406
    461407        G1 = Geospatial_data(points1, attributes1, geo_ref1)
    … …  
    465411        assert G1.attributes.has_key('elevation')
    466412        assert num.allclose(G1.attributes['depth'], [2, 4.7])
    467         assert num.allclose(G1.attributes['elevation'], [6.1, 5])       
    468        
     413        assert num.allclose(G1.attributes['elevation'], [6.1, 5])
     414
    469415        G2 = Geospatial_data(points2, attributes2, geo_ref2)
    470416        assert num.allclose(G2.get_geo_reference().get_xllcorner(), 0.1)
    … …  
    473419        assert G2.attributes.has_key('elevation')
    474420        assert num.allclose(G2.attributes['depth'], [-2.3, 4])
    475         assert num.allclose(G2.attributes['elevation'], [2.5, 1])       
     421        assert num.allclose(G2.attributes['elevation'], [2.5, 1])
    476422
    477423        #Check that absolute values are as expected
    … …  
    480426
    481427        P2 = G2.get_data_points(absolute=True)
    482         assert num.allclose(P2, [[5.2, 12.1], [6.2, 9.3]])       
     428        assert num.allclose(P2, [[5.2, 12.1], [6.2, 9.3]])
    483429
    484430        # Normal add
    485431        G = G1 + None
    486 
    487432        assert G.attributes.has_key('depth')
    488433        assert G.attributes.has_key('elevation')
    489434        assert num.allclose(G.attributes['depth'], [2, 4.7])
    490         assert num.allclose(G.attributes['elevation'], [6.1, 5])       
     435        assert num.allclose(G.attributes['elevation'], [6.1, 5])
    491436
    492437        # Points are now absolute.
    493438        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
    494439        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
    495         P = G.get_data_points(absolute=True)       
    496         assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]])
    497 
     440        P = G.get_data_points(absolute=True)
     441        msg = 'P=%s' % str(P)
     442        assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]]), msg
    498443
    499444        G = G2 + None
    … …  
    501446        assert G.attributes.has_key('elevation')
    502447        assert num.allclose(G.attributes['depth'], [-2.3, 4])
    503         assert num.allclose(G.attributes['elevation'], [2.5, 1])       
    504 
     448        assert num.allclose(G.attributes['elevation'], [2.5, 1])
    505449        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
    506450        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
    507         P = G.get_data_points(absolute=True)       
     451
     452        P = G.get_data_points(absolute=True)
    508453        assert num.allclose(P, [[5.2, 12.1], [6.2, 9.3]])
    509        
    510 
    511454
    512455        # Reverse add
    513456        G = None + G1
    514 
    515457        assert G.attributes.has_key('depth')
    516458        assert G.attributes.has_key('elevation')
    517459        assert num.allclose(G.attributes['depth'], [2, 4.7])
    518         assert num.allclose(G.attributes['elevation'], [6.1, 5])       
     460        assert num.allclose(G.attributes['elevation'], [6.1, 5])
    519461
    520462        # Points are now absolute.
    521463        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
    522464        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
    523         P = G.get_data_points(absolute=True)       
    524         assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]])       
    525 
     465
     466        P = G.get_data_points(absolute=True)
     467        assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]])
    526468
    527469        G = None + G2
    … …  
    529471        assert G.attributes.has_key('elevation')
    530472        assert num.allclose(G.attributes['depth'], [-2.3, 4])
    531         assert num.allclose(G.attributes['elevation'], [2.5, 1])       
     473        assert num.allclose(G.attributes['elevation'], [2.5, 1])
    532474
    533475        assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0)
    534476        assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0)
    535         P = G.get_data_points(absolute=True)       
     477
     478        P = G.get_data_points(absolute=True)
    536479        assert num.allclose(P, [[5.2, 12.1], [6.2, 9.3]])
    537480
    538        
    539 
    540        
    541                            
    542        
    543481    def test_clip0(self):
    544         """test_clip0(self):
    545        
     482        '''test_clip0(self):
     483
    546484        Test that point sets can be clipped by a polygon
    547         """
    548        
     485        '''
     486
    549487        from anuga.coordinate_transforms.geo_reference import Geo_reference
    550        
     488
    551489        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
    552490                  [0, 0], [2.4, 3.3]]
    553491        G = Geospatial_data(points)
    554492
    555         # First try the unit square   
    556         U = [[0,0], [1,0], [1,1], [0,1]]
    557         assert num.allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]])
     493        # First try the unit square
     494        U = [[0,0], [1,0], [1,1], [0,1]]
     495        assert num.allclose(G.clip(U).get_data_points(),
     496                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
    558497
    559498        # Then a more complex polygon
    560499        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
    561         points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
     500        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
     501                  [0.5, 1.5], [0.5, -0.5]]
    562502        G = Geospatial_data(points)
    563503
    564504        assert num.allclose(G.clip(polygon).get_data_points(),
    565                         [[0.5, 0.5], [1, -0.5], [1.5, 0]])
     505                            [[0.5, 0.5], [1, -0.5], [1.5, 0]])
    566506
    567507    def test_clip0_with_attributes(self):
    568         """test_clip0_with_attributes(self):
    569        
     508        '''test_clip0_with_attributes(self):
     509
    570510        Test that point sets with attributes can be clipped by a polygon
    571         """
    572        
     511        '''
     512
    573513        from anuga.coordinate_transforms.geo_reference import Geo_reference
    574        
     514
    575515        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
    576516                  [0, 0], [2.4, 3.3]]
    … …  
    579519        att_dict = {'att1': attributes,
    580520                    'att2': num.array(attributes)+1}
    581        
     521
    582522        G = Geospatial_data(points, att_dict)
    583523
    584         # First try the unit square   
    585         U = [[0,0], [1,0], [1,1], [0,1]]
    586         assert num.allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]])
     524        # First try the unit square
     525        U = [[0,0], [1,0], [1,1], [0,1]]
     526        assert num.allclose(G.clip(U).get_data_points(),
     527                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
    587528        assert num.allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
    588         assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])               
     529        assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])
    589530
    590531        # Then a more complex polygon
    591532        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
    592         points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
     533        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
     534                  [0.5, 1.5], [0.5, -0.5]]
    593535
    594536        # This time just one attribute
    595537        attributes = [2, -4, 5, 76, -2, 0.1]
    596538        G = Geospatial_data(points, attributes)
    597 
    598539        assert num.allclose(G.clip(polygon).get_data_points(),
    599                         [[0.5, 0.5], [1, -0.5], [1.5, 0]])
     540                            [[0.5, 0.5], [1, -0.5], [1.5, 0]])
    600541        assert num.allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
    601        
    602542
    603543    def test_clip1(self):
    604         """test_clip1(self):
    605        
     544        '''test_clip1(self):
     545
    606546        Test that point sets can be clipped by a polygon given as
    607547        another Geospatial dataset
    608         """
    609        
     548        '''
     549
    610550        from anuga.coordinate_transforms.geo_reference import Geo_reference
    611        
     551
    612552        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
    613553                  [0, 0], [2.4, 3.3]]
    … …  
    616556                    'att2': num.array(attributes)+1}
    617557        G = Geospatial_data(points, att_dict)
    618        
    619         # First try the unit square   
    620         U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
     558
     559        # First try the unit square
     560        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
    621561        assert num.allclose(G.clip(U).get_data_points(),
    622                         [[0.2, 0.5], [0.4, 0.3], [0, 0]])
    623 
     562                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
    624563        assert num.allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
    625         assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])                       
    626        
     564        assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])
     565
    627566        # Then a more complex polygon
    628         points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
    629         attributes = [2, -4, 5, 76, -2, 0.1]       
     567        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
     568                  [0.5, 1.5], [0.5, -0.5]]
     569        attributes = [2, -4, 5, 76, -2, 0.1]
    630570        G = Geospatial_data(points, attributes)
    631         polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
    632        
     571        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1],
     572                                   [2,1], [0,1]])
    633573
    634574        assert num.allclose(G.clip(polygon).get_data_points(),
    635575                            [[0.5, 0.5], [1, -0.5], [1.5, 0]])
    636576        assert num.allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
    637        
    638577
    639578    def test_clip0_outside(self):
    640         """test_clip0_outside(self):
    641        
     579        '''test_clip0_outside(self):
     580
    642581        Test that point sets can be clipped outside of a polygon
    643         """
    644        
     582        '''
     583
    645584        from anuga.coordinate_transforms.geo_reference import Geo_reference
    646        
     585
    647586        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
    648587                  [0, 0], [2.4, 3.3]]
    649         attributes = [2, -4, 5, 76, -2, 0.1, 3]       
     588        attributes = [2, -4, 5, 76, -2, 0.1, 3]
    650589        G = Geospatial_data(points, attributes)
    651590
    652         # First try the unit square   
     591        # First try the unit square
    653592        U = [[0,0], [1,0], [1,1], [0,1]]
    654593        assert num.allclose(G.clip_outside(U).get_data_points(),
    655594                            [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
    656         #print G.clip_outside(U).get_attributes()
    657         assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])       
    658        
     595        assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])
    659596
    660597        # Then a more complex polygon
    661598        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
    662         points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
    663         attributes = [2, -4, 5, 76, -2, 0.1]       
     599        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
     600                  [0.5, 1.5], [0.5, -0.5]]
     601        attributes = [2, -4, 5, 76, -2, 0.1]
    664602        G = Geospatial_data(points, attributes)
    665 
    666603        assert num.allclose(G.clip_outside(polygon).get_data_points(),
    667604                            [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
    668         assert num.allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])               
    669 
     605        assert num.allclose(G.clip_outside(polygon).get_attributes(),
     606                            [2, -2, 0.1])
    670607
    671608    def test_clip1_outside(self):
    672         """test_clip1_outside(self):
    673        
     609        '''test_clip1_outside(self):
     610
    674611        Test that point sets can be clipped outside of a polygon given as
    675612        another Geospatial dataset
    676         """
    677        
     613        '''
     614
    678615        from anuga.coordinate_transforms.geo_reference import Geo_reference
    679        
     616
    680617        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
    681618                  [0, 0], [2.4, 3.3]]
    682         attributes = [2, -4, 5, 76, -2, 0.1, 3]       
    683         G = Geospatial_data(points, attributes)       
    684 
    685         # First try the unit square   
    686         U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
     619        attributes = [2, -4, 5, 76, -2, 0.1, 3]
     620        G = Geospatial_data(points, attributes)
     621
     622        # First try the unit square
     623        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
    687624        assert num.allclose(G.clip_outside(U).get_data_points(),
    688625                            [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
    689         assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])       
     626        assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])
    690627
    691628        # Then a more complex polygon
    692         points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
    693         attributes = [2, -4, 5, 76, -2, 0.1]       
     629        points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0],
     630                  [0.5, 1.5], [0.5, -0.5]]
     631        attributes = [2, -4, 5, 76, -2, 0.1]
    694632        G = Geospatial_data(points, attributes)
    695 
    696         polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
    697        
    698 
     633        polygon = Geospatial_data([[0, 0], [1, 0], [0.5, -1], [2, -1],
     634                                   [2, 1], [0, 1]])
    699635        assert num.allclose(G.clip_outside(polygon).get_data_points(),
    700636                            [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
    701         assert num.allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])
    702        
    703 
     637        assert num.allclose(G.clip_outside(polygon).get_attributes(),
     638                            [2, -2, 0.1])
    704639
    705640    def test_clip1_inside_outside(self):
    706         """test_clip1_inside_outside(self):
    707        
     641        '''test_clip1_inside_outside(self):
     642
    708643        Test that point sets can be clipped outside of a polygon given as
    709644        another Geospatial dataset
    710         """
    711        
     645        '''
     646
    712647        from anuga.coordinate_transforms.geo_reference import Geo_reference
    713        
     648
    714649        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
    715650                  [0, 0], [2.4, 3.3]]
    716         attributes = [2, -4, 5, 76, -2, 0.1, 3]       
     651        attributes = [2, -4, 5, 76, -2, 0.1, 3]
    717652        G = Geospatial_data(points, attributes)
    718653
    719         # First try the unit square   
    720         U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
     654        # First try the unit square
     655        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
    721656        G1 = G.clip(U)
    722         assert num.allclose(G1.get_data_points(),[[0.2, 0.5], [0.4, 0.3], [0, 0]])
     657        assert num.allclose(G1.get_data_points(),
     658                            [[0.2, 0.5], [0.4, 0.3], [0, 0]])
    723659        assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])
    724        
    725660        G2 = G.clip_outside(U)
    726661        assert num.allclose(G2.get_data_points(),[[-1, 4], [1.0, 2.1],
    727662                                                  [3.0, 5.3], [2.4, 3.3]])
    728         assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])               
    729 
    730        
     663        assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])
     664
    731665        # New ordering
    732         new_points = [[0.2, 0.5], [0.4, 0.3], [0, 0]] +\
    733                      [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]
    734 
    735         new_attributes = [-4, 76, 0.1, 2, 5, -2, 3]                 
    736        
     666        new_points = [[0.2, 0.5], [0.4, 0.3], [0, 0], [-1, 4],
     667                      [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]
     668        new_attributes = [-4, 76, 0.1, 2, 5, -2, 3]
     669
    737670        assert num.allclose((G1+G2).get_data_points(), new_points)
    738671        assert num.allclose((G1+G2).get_attributes(), new_attributes)
    … …  
    742675        G.export_points_file(FN)
    743676
    744 
    745677        # Read it back in
    746678        G3 = Geospatial_data(FN)
    747679
    748 
    749680        # Check result
    750         assert num.allclose(G3.get_data_points(), new_points)       
    751         assert num.allclose(G3.get_attributes(), new_attributes)       
    752        
     681        assert num.allclose(G3.get_data_points(), new_points)
     682        assert num.allclose(G3.get_attributes(), new_attributes)
     683
    753684        os.remove(FN)
    754685
    755        
    756686    def test_load_csv(self):
    757        
    758         import os
    759         import tempfile
    760        
    761         fileName = tempfile.mktemp(".csv")
    762         file = open(fileName,"w")
    763         file.write("x,y,elevation speed \n\
     687        fileName = tempfile.mktemp('.csv')
     688        file = open(fileName,'w')
     689        file.write('x,y,elevation speed \n\
    7646901.0 0.0 10.0 0.0\n\
    7656910.0 1.0 0.0 10.0\n\
    766 1.0 0.0 10.4 40.0\n")
     6921.0 0.0 10.4 40.0\n')
    767693        file.close()
    768         #print fileName
     694
    769695        results = Geospatial_data(fileName)
    770696        os.remove(fileName)
    771 #        print 'data', results.get_data_points()
    772         assert num.allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],
    773                                                         [1.0, 0.0]])
     697        assert num.allclose(results.get_data_points(),
     698                            [[1.0, 0.0],[0.0, 1.0], [1.0, 0.0]])
    774699        assert num.allclose(results.get_attributes(attribute_name='elevation'),
    775700                            [10.0, 0.0, 10.4])
    … …  
    777702                            [0.0, 10.0, 40.0])
    778703
    779 
    780   ###################### .CSV ##############################
     704################################################################################
     705# Test CSV files
     706################################################################################
    781707
    782708    def test_load_csv_lat_long_bad_blocking(self):
    783         """
    784         test_load_csv_lat_long_bad_blocking(self):
     709        '''test_load_csv_lat_long_bad_blocking(self):
    785710        Different zones in Geo references
    786         """
    787         fileName = tempfile.mktemp(".csv")
    788         file = open(fileName,"w")
    789         file.write("Lati,LONG,z \n\
     711        '''
     712
     713        fileName = tempfile.mktemp('.csv')
     714        file = open(fileName, 'w')
     715        file.write('Lati,LONG,z \n\
    790716-25.0,180.0,452.688000\n\
    791 -34,150.0,459.126000\n")
     717-34,150.0,459.126000\n')
    792718        file.close()
    793        
     719
    794720        results = Geospatial_data(fileName, max_read_lines=1,
    795721                                  load_file_now=False)
    796        
    797         #for i in results:
    798         #    pass
     722
    799723        try:
    800724            for i in results:
    … …  
    804728        else:
    805729            msg = 'Different zones in Geo references not caught.'
    806             raise msg       
    807        
    808         os.remove(fileName)
    809        
     730            raise Exception, msg
     731
     732        os.remove(fileName)
     733
    810734    def test_load_csv(self):
    811        
    812         fileName = tempfile.mktemp(".txt")
    813         file = open(fileName,"w")
    814         file.write(" x,y, elevation ,  speed \n\
    815 1.0, 0.0, 10.0, 0.0\n\
    816 0.0, 1.0, 0.0, 10.0\n\
    817 1.0, 0.0 ,10.4, 40.0\n")
     735        fileName = tempfile.mktemp('.txt')
     736        file = open(fileName, 'w')
     737        file.write(' x,y, elevation ,  speed \n\
     7381.0, 0.0, 10.0, 0.0\n\
     7390.0, 1.0, 0.0, 10.0\n\
     7401.0, 0.0 ,10.4, 40.0\n')
    818741        file.close()
    819742
    820743        results = Geospatial_data(fileName, max_read_lines=2)
    821744
    822 
    823         assert num.allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
    824         assert num.allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4])
    825         assert num.allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0])
     745        assert num.allclose(results.get_data_points(),
     746                            [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])
     747        assert num.allclose(results.get_attributes(attribute_name='elevation'),
     748                            [10.0, 0.0, 10.4])
     749        assert num.allclose(results.get_attributes(attribute_name='speed'),
     750                            [0.0, 10.0, 40.0])
    826751
    827752        # Blocking
    … …  
    829754        for i in results:
    830755            geo_list.append(i)
    831            
     756
    832757        assert num.allclose(geo_list[0].get_data_points(),
    833758                            [[1.0, 0.0],[0.0, 1.0]])
    834 
    835759        assert num.allclose(geo_list[0].get_attributes(attribute_name='elevation'),
    836760                            [10.0, 0.0])
    837761        assert num.allclose(geo_list[1].get_data_points(),
    838                             [[1.0, 0.0]])       
     762                            [[1.0, 0.0]])
    839763        assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'),
    840764                            [10.4])
    841            
    842         os.remove(fileName)         
    843        
     765
     766        os.remove(fileName)
     767
    844768    def test_load_csv_bad(self):
    845         """ test_load_csv_bad(self):
     769        '''test_load_csv_bad(self):
    846770        header column, body column missmatch
    847771        (Format error)
    848         """
    849         import os
    850        
    851         fileName = tempfile.mktemp(".txt")
    852         file = open(fileName,"w")
    853         file.write(" elevation ,  speed \n\
    854 1.0, 0.0, 10.0, 0.0\n\
    855 0.0, 1.0, 0.0, 10.0\n\
    856 1.0, 0.0 ,10.4, 40.0\n")
     772        '''
     773
     774        fileName = tempfile.mktemp('.txt')
     775        file = open(fileName, 'w')
     776        file.write(' elevation ,  speed \n\
     7771.0, 0.0, 10.0, 0.0\n\
     7780.0, 1.0, 0.0, 10.0\n\
     7791.0, 0.0 ,10.4, 40.0\n')
    857780        file.close()
    858781
    … …  
    862785        # Blocking
    863786        geo_list = []
    864         #for i in results:
    865         #    geo_list.append(i)
    866787        try:
    867788            for i in results:
    … …  
    870791            pass
    871792        else:
    872             msg = 'bad file did not raise error!'
    873             raise msg
     793            msg = 'Bad file did not raise error!'
     794            raise Exception, msg
    874795        os.remove(fileName)
    875796
    876797    def test_load_csv_badII(self):
    877         """ test_load_csv_bad(self):
     798        '''test_load_csv_bad(self):
    878799        header column, body column missmatch
    879800        (Format error)
    880         """
    881         import os
    882        
    883         fileName = tempfile.mktemp(".txt")
    884         file = open(fileName,"w")
    885         file.write(" x,y,elevation ,  speed \n\
     801        '''
     802
     803        fileName = tempfile.mktemp('.txt')
     804        file = open(fileName, 'w')
     805        file.write(' x,y,elevation ,  speed \n\
    8868061.0, 0.0, 10.0, 0.0\n\
    8878070.0, 1.0, 0.0, 10\n\
    888 1.0, 0.0 ,10.4 yeah\n")
     8081.0, 0.0 ,10.4 yeah\n')
    889809        file.close()
    890810
    … …  
    894814        # Blocking
    895815        geo_list = []
    896         #for i in results:
    897         #    geo_list.append(i)
    898816        try:
    899817            for i in results:
    … …  
    902820            pass
    903821        else:
    904             msg = 'bad file did not raise error!'
    905             raise msg
     822            msg = 'Bad file did not raise error!'
     823            raise Exception, msg
    906824        os.remove(fileName)
    907825
    908826    def test_load_csv_badIII(self):
    909         """ test_load_csv_bad(self):
     827        '''test_load_csv_bad(self):
    910828        space delimited
    911         """
    912         import os
    913        
    914         fileName = tempfile.mktemp(".txt")
    915         file = open(fileName,"w")
    916         file.write(" x y elevation   speed \n\
     829        '''
     830
     831        fileName = tempfile.mktemp('.txt')
     832        file = open(fileName, 'w')
     833        file.write(' x y elevation   speed \n\
    9178341. 0.0 10.0 0.0\n\
    9188350.0 1.0 0.0 10.0\n\
    919 1.0 0.0 10.4 40.0\n")
     8361.0 0.0 10.4 40.0\n')
    920837        file.close()
    921838
    … …  
    926843            pass
    927844        else:
    928             msg = 'bad file did not raise error!'
    929             raise msg
    930         os.remove(fileName)
    931        
     845            msg = 'Bad file did not raise error!'
     846            raise Exception, msg
     847        os.remove(fileName)
     848
    932849    def test_load_csv_badIV(self):
    933         """ test_load_csv_bad(self):
     850        ''' test_load_csv_bad(self):
    934851        header column, body column missmatch
    935852        (Format error)
    936         """
    937         import os
    938        
    939         fileName = tempfile.mktemp(".txt")
    940         file = open(fileName,"w")
    941         file.write(" x,y,elevation ,  speed \n\
     853        '''
     854
     855        fileName = tempfile.mktemp('.txt')
     856        file = open(fileName, 'w')
     857        file.write(' x,y,elevation ,  speed \n\
    9428581.0, 0.0, 10.0, wow\n\
    9438590.0, 1.0, 0.0, ha\n\
    944 1.0, 0.0 ,10.4, yeah\n")
     8601.0, 0.0 ,10.4, yeah\n')
    945861        file.close()
    946862
    … …  
    950866        # Blocking
    951867        geo_list = []
    952         #for i in results:
    953          #   geo_list.append(i)
    954868        try:
    955869            for i in results:
    … …  
    958872            pass
    959873        else:
    960             msg = 'bad file did not raise error!'
    961             raise msg
     874            msg = 'Bad file did not raise error!'
     875            raise Exception, msg
    962876        os.remove(fileName)
    963877
    964878    def test_load_pts_blocking(self):
    965879        #This is pts!
    966        
    967         import os
    968        
    969         fileName = tempfile.mktemp(".txt")
    970         file = open(fileName,"w")
    971         file.write(" x,y, elevation ,  speed \n\
    972 1.0, 0.0, 10.0, 0.0\n\
    973 0.0, 1.0, 0.0, 10.0\n\
    974 1.0, 0.0 ,10.4, 40.0\n")
     880        fileName = tempfile.mktemp('.txt')
     881        file = open(fileName, 'w')
     882        file.write(' x,y, elevation ,  speed \n\
     8831.0, 0.0, 10.0, 0.0\n\
     8840.0, 1.0, 0.0, 10.0\n\
     8851.0, 0.0 ,10.4, 40.0\n')
    975886        file.close()
    976887
    977         pts_file = tempfile.mktemp(".pts")
    978        
     888        pts_file = tempfile.mktemp('.pts')
     889
    979890        convert = Geospatial_data(fileName)
    980891        convert.export_points_file(pts_file)
    981892        results = Geospatial_data(pts_file, max_read_lines=2)
    982893
    983         assert num.allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],
    984                                                         [1.0, 0.0]])
     894        assert num.allclose(results.get_data_points(),
     895                            [[1.0, 0.0],[0.0, 1.0], [1.0, 0.0]])
    985896        assert num.allclose(results.get_attributes(attribute_name='elevation'),
    986897                            [10.0, 0.0, 10.4])
    … …  
    991902        geo_list = []
    992903        for i in results:
    993             geo_list.append(i) 
     904            geo_list.append(i)
    994905        assert num.allclose(geo_list[0].get_data_points(),
    995906                            [[1.0, 0.0],[0.0, 1.0]])
    … …  
    997908                            [10.0, 0.0])
    998909        assert num.allclose(geo_list[1].get_data_points(),
    999                             [[1.0, 0.0]])       
     910                            [[1.0, 0.0]])
    1000911        assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'),
    1001912                            [10.4])
    1002            
    1003         os.remove(fileName) 
    1004         os.remove(pts_file)               
     913
     914        os.remove(fileName)
     915        os.remove(pts_file)
    1005916
    1006917    def verbose_test_load_pts_blocking(self):
    1007        
    1008         import os
    1009        
    1010         fileName = tempfile.mktemp(".txt")
    1011         file = open(fileName,"w")
    1012         file.write(" x,y, elevation ,  speed \n\
    1013