Ignore:
Timestamp:
Jun 30, 2009, 2:07:41 PM (16 years ago)
Author:
ole
Message:

Merged numpy branch back into the trunk.

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

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

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

For all merges, I used numpy version where conflicts existed

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

Location:
anuga_core/source/anuga/coordinate_transforms
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/coordinate_transforms/geo_reference.py

    r7062 r7276  
    77#and unit test
    88
    9 
    109import types, sys
    1110import copy
     11
    1212from anuga.utilities.numerical_tools import ensure_numeric
    13 from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, ParsingError, \
    14      ShapeError
    15 
    16 import Numeric as num
     13from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, \
     14                                             ParsingError, ShapeError
     15from anuga.config import netcdf_float, netcdf_int, netcdf_float32
     16
     17import numpy as num
    1718
    1819
    1920DEFAULT_ZONE = -1
    20 TITLE = '#geo reference' + "\n" #this title is referred to in the test format
     21TITLE = '#geo reference' + "\n" # this title is referred to in the test format
    2122
    2223DEFAULT_PROJECTION = 'UTM'
    … …  
    2425DEFAULT_UNITS = 'm'
    2526DEFAULT_FALSE_EASTING = 500000
    26 DEFAULT_FALSE_NORTHING = 10000000 #Default for southern hemisphere
    27 
     27DEFAULT_FALSE_NORTHING = 10000000    # Default for southern hemisphere
     28
     29
     30##
     31# @brief A class for ...
    2832class Geo_reference:
    2933    """
     34    Attributes of the Geo_reference class:
     35        .zone           The UTM zone (default is -1)
     36        .false_easting  ??
     37        .false_northing ??
     38        .datum          The Datum used (default is wgs84)
     39        .projection     The projection used (default is 'UTM')
     40        .units          The units of measure used (default metres)
     41        .xllcorner      The X coord of origin (default is 0.0 wrt UTM grid)
     42        .yllcorner      The y coord of origin (default is 0.0 wrt UTM grid)
     43        .is_absolute    ??
     44
    3045    """
    3146
     47    ##
     48    # @brief Instantiate an instance of class Geo_reference.
     49    # @param zone The UTM zone.
     50    # @param xllcorner X coord of origin of georef.
     51    # @param yllcorner Y coord of origin of georef.
     52    # @param datum ??
     53    # @param projection The projection used (default UTM).
     54    # @param units Units used in measuring distance (default m).
     55    # @param false_easting ??
     56    # @param false_northing ??
     57    # @param NetCDFObject NetCDF file *handle* to write to.
     58    # @param ASCIIFile ASCII text file *handle* to write to.
     59    # @param read_title Title of the georeference text.
    3260    def __init__(self,
    33                  zone = DEFAULT_ZONE,
    34                  xllcorner = 0.0,
    35                  yllcorner = 0.0,
    36                  datum = DEFAULT_DATUM,
    37                  projection = DEFAULT_PROJECTION,
    38                  units = DEFAULT_UNITS,
    39                  false_easting = DEFAULT_FALSE_EASTING,
    40                  false_northing = DEFAULT_FALSE_NORTHING,
     61                 zone=DEFAULT_ZONE,
     62                 xllcorner=0.0,
     63                 yllcorner=0.0,
     64                 datum=DEFAULT_DATUM,
     65                 projection=DEFAULT_PROJECTION,
     66                 units=DEFAULT_UNITS,
     67                 false_easting=DEFAULT_FALSE_EASTING,
     68                 false_northing=DEFAULT_FALSE_NORTHING,
    4169                 NetCDFObject=None,
    4270                 ASCIIFile=None,
    … …  
    4472        """
    4573        input:
    46         NetCDFObject - a handle to the netCDF file to be written to 
     74        NetCDFObject - a handle to the netCDF file to be written to
    4775        ASCIIFile - a handle to the text file
    4876        read_title - the title of the georeference text, if it was read in.
    … …  
    5583         must be the default info, since ANUGA assumes it isn't
    5684         changing.
    57          """
     85        """
     86
    5887        if zone is None:
    5988            zone = DEFAULT_ZONE
    … …  
    6998        if NetCDFObject is not None:
    7099            self.read_NetCDF(NetCDFObject)
    71  
     100
    72101        if ASCIIFile is not None:
    73102            self.read_ASCII(ASCIIFile, read_title=read_title)
    … …  
    79108    def get_xllcorner(self):
    80109        return self.xllcorner
    81    
     110
     111    ##
     112    # @brief Get the Y coordinate of the origin of this georef.
    82113    def get_yllcorner(self):
    83114        return self.yllcorner
    84    
     115
     116    ##
     117    # @brief Get the zone of this georef.
    85118    def get_zone(self):
    86119        return self.zone
    87        
     120
     121    ##
     122    # @brief Write <something> to an open NetCDF file.
     123    # @param outfile Handle to open NetCDF file.
    88124    def write_NetCDF(self, outfile):
    89         outfile.xllcorner = self.xllcorner 
     125        outfile.xllcorner = self.xllcorner
    90126        outfile.yllcorner = self.yllcorner
    91127        outfile.zone = self.zone
    … …  
    94130        outfile.false_northing = self.false_northing
    95131
    96         outfile.datum = self.datum       
     132        outfile.datum = self.datum
    97133        outfile.projection = self.projection
    98134        outfile.units = self.units
    99135
     136    ##
     137    # @brief Read data from an open NetCDF file.
     138    # @param infile Handle to open NetCDF file.
    100139    def read_NetCDF(self, infile):
    101140        self.xllcorner = float(infile.xllcorner[0])
    … …  
    113152            pass
    114153
    115         if (self.false_easting != DEFAULT_FALSE_EASTING):
    116             print "WARNING: False easting of %f specified." %self.false_easting
    117             print "Default false easting is %f." %DEFAULT_FALSE_EASTING
     154        if self.false_easting != DEFAULT_FALSE_EASTING:
     155            print "WARNING: False easting of %f specified." % self.false_easting
     156            print "Default false easting is %f." % DEFAULT_FALSE_EASTING
    118157            print "ANUGA does not correct for differences in False Eastings."
    119            
    120         if (self.false_northing != DEFAULT_FALSE_NORTHING):
    121             print "WARNING: False northing of %f specified." \
    122                   %self.false_northing
    123             print "Default false northing is %f." %DEFAULT_FALSE_NORTHING
     158
     159        if self.false_northing != DEFAULT_FALSE_NORTHING:
     160            print ("WARNING: False northing of %f specified."
     161                   % self.false_northing)
     162            print "Default false northing is %f." % DEFAULT_FALSE_NORTHING
    124163            print "ANUGA does not correct for differences in False Northings."
    125            
    126         if (self.datum.upper() != DEFAULT_DATUM.upper()):
    127             print "WARNING: Datum of %s specified." \
    128                   %self.datum
    129             print "Default Datum is %s." %DEFAULT_DATUM
     164
     165        if self.datum.upper() != DEFAULT_DATUM.upper():
     166            print "WARNING: Datum of %s specified." % self.datum
     167            print "Default Datum is %s." % DEFAULT_DATUM
    130168            print "ANUGA does not correct for differences in datums."
    131            
    132         if (self.projection.upper() != DEFAULT_PROJECTION.upper()):
    133             print "WARNING: Projection of %s specified." \
    134                   %self.projection
    135             print "Default Projection is %s." %DEFAULT_PROJECTION
     169
     170        if self.projection.upper() != DEFAULT_PROJECTION.upper():
     171            print "WARNING: Projection of %s specified." % self.projection
     172            print "Default Projection is %s." % DEFAULT_PROJECTION
    136173            print "ANUGA does not correct for differences in Projection."
    137            
    138         if (self.units.upper() != DEFAULT_UNITS.upper()):
    139             print "WARNING: Units of %s specified." \
    140                   %self.units
    141             print "Default units is %s." %DEFAULT_UNITS
     174
     175        if self.units.upper() != DEFAULT_UNITS.upper():
     176            print "WARNING: Units of %s specified." % self.units
     177            print "Default units is %s." % DEFAULT_UNITS
    142178            print "ANUGA does not correct for differences in units."
    143        
    144     ### ASCII files with geo-refs are currently not used   
     179
     180################################################################################
     181# ASCII files with geo-refs are currently not used
     182################################################################################
     183
     184    ##
     185    # @brief Write georef data to an open text file.
     186    # @param fd Handle to open text file.
    145187    def write_ASCII(self, fd):
    146188        fd.write(TITLE)
    147         fd.write(str(self.zone) + "\n") 
    148         fd.write(str(self.xllcorner) + "\n") 
     189        fd.write(str(self.zone) + "\n")
     190        fd.write(str(self.xllcorner) + "\n")
    149191        fd.write(str(self.yllcorner) + "\n")
    150192
    151     def read_ASCII(self, fd,read_title=None):
     193    ##
     194    # @brief Read georef data from an open text file.
     195    # @param fd Handle to open text file.
     196    def read_ASCII(self, fd, read_title=None):
    152197        try:
    153198            if read_title == None:
    154                 read_title = fd.readline() # remove the title line
     199                read_title = fd.readline()     # remove the title line
    155200            if read_title[0:2].upper() != TITLE[0:2].upper():
    156                 msg = 'File error.  Expecting line: %s.  Got this line: %s' \
    157                       %(TITLE, read_title)
    158                 raise TitleError, msg 
     201                msg = ('File error.  Expecting line: %s.  Got this line: %s'
     202                       % (TITLE, read_title))
     203                raise TitleError, msg
    159204            self.zone = int(fd.readline())
    160205            self.xllcorner = float(fd.readline())
    161206            self.yllcorner = float(fd.readline())
    162207        except SyntaxError:
    163                 msg = 'File error.  Got syntax error while parsing geo reference'
    164                 raise ParsingError, msg
    165            
     208            msg = 'File error.  Got syntax error while parsing geo reference'
     209            raise ParsingError, msg
     210
    166211        # Fix some assertion failures
    167         if(type(self.zone) == num.ArrayType and self.zone.shape == ()):
     212        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
    168213            self.zone = self.zone[0]
    169         if type(self.xllcorner) == num.ArrayType and self.xllcorner.shape == ():
     214        if (isinstance(self.xllcorner, num.ndarray) and
     215                self.xllcorner.shape == ()):
    170216            self.xllcorner = self.xllcorner[0]
    171         if type(self.yllcorner) == num.ArrayType and self.yllcorner.shape == ():
     217        if (isinstance(self.yllcorner, num.ndarray) and
     218                self.yllcorner.shape == ()):
    172219            self.yllcorner = self.yllcorner[0]
    173    
     220
    174221        assert (type(self.xllcorner) == types.FloatType)
    175222        assert (type(self.yllcorner) == types.FloatType)
    176223        assert (type(self.zone) == types.IntType)
    177        
    178        
    179     def change_points_geo_ref(self, points,
    180                               points_geo_ref=None):
    181         """
    182         Change the geo reference of a list or Numeric array of points to
     224
     225################################################################################
     226
     227    ##
     228    # @brief Change points to be absolute wrt new georef 'points_geo_ref'.
     229    # @param points The points to change.
     230    # @param points_geo_ref The new georef to make points absolute wrt.
     231    # @return The changed points.
     232    # @note If 'points' is a list then a changed list is returned.
     233    def change_points_geo_ref(self, points, points_geo_ref=None):
     234        """Change the geo reference of a list or numeric array of points to
    183235        be this reference.(The reference used for this object)
    184236        If the points do not have a geo ref, assume 'absolute' values
    185237        """
    186238        import copy
    187         
    188         is_list = False
    189         if type(points) == types.ListType:
    190             is_list = True
    191 
    192         points = ensure_numeric(points, num.Float)
    193        
     239       
     240        # remember if we got a list
     241        is_list = isinstance(points, list)
     242
     243        points = ensure_numeric(points, num.float)
     244
     245        # sanity checks
    194246        if len(points.shape) == 1:
    195247            #One point has been passed
    … …  
    206258        assert points.shape[1] == 2, msg               
    207259
    208        
    209260        # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
    210261        # are identical in the two geo refs.   
    211262        if points_geo_ref is not self:
    212263            # If georeferences are different
    213        
    214264            points = copy.copy(points) # Don't destroy input                   
    215265            if not points_geo_ref is None:
    … …  
    222272            points[:,1] -= self.yllcorner
    223273
    224            
    225274        if is_list:
    226275            points = points.tolist()
    … …  
    228277        return points
    229278
    230    
    231279    def is_absolute(self):
    232280        """Return True if xllcorner==yllcorner==0 indicating that points
    … …  
    241289        # With the flag method fitting is much faster (18 Mar 2009).
    242290
    243        
    244291        # FIXME(Ole): HACK to be able to reuse data already cached (18 Mar 2009).
    245292        # Remove at some point
    246293        if not hasattr(self, 'absolute'):
    247294            self.absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
    248 
    249295           
    250296        # Return absolute flag   
    251297        return self.absolute
    252        
    253    
     298
    254299    def get_absolute(self, points):
    255         """
    256         Given a set of points geo referenced to this instance,
     300        """Given a set of points geo referenced to this instance,
    257301        return the points as absolute values.
    258302        """
    259303
    260         is_list = False
    261         if type(points) == types.ListType:
    262             is_list = True
    263 
    264         points = ensure_numeric(points, num.Float)
     304        # remember if we got a list
     305        is_list = isinstance(points, list)
     306
     307        points = ensure_numeric(points, num.float)
    265308        if len(points.shape) == 1:
    266309            # One point has been passed
    … …  
    288331        return points
    289332
    290 
     333    ##
     334    # @brief Convert points to relative measurement.
     335    # @param points Points to convert to relative measurements.
     336    # @return A set of points relative to the geo_reference instance.
    291337    def get_relative(self, points):
    292         """
    293         Given a set of points in absolute UTM coordinates,
     338        """Given a set of points in absolute UTM coordinates,
    294339        make them relative to this geo_reference instance,
    295340        return the points as relative values.
    … …  
    298343        """
    299344
    300         is_list = False
    301         if type(points) == types.ListType:
    302             is_list = True
    303 
    304         points = ensure_numeric(points, num.Float)
     345        # remember if we got a list
     346        is_list = isinstance(points, list)
     347
     348        points = ensure_numeric(points, num.float)
    305349        if len(points.shape) == 1:
    306350            #One point has been passed
    … …  
    308352            if not len(points) == 2:
    309353                raise ShapeError, msg   
    310                 #points = reshape(points, (1,2))
    311 
    312 
    313         msg = 'Input must be an N x 2 array or list of (x,y) values. '
    314         msg += 'I got an %d x %d array' %points.shape   
     354
    315355        if not points.shape[1] == 2:
     356            msg = ('Input must be an N x 2 array or list of (x,y) values. '
     357                   'I got an %d x %d array' % points.shape)
    316358            raise ShapeError, msg   
    317            
    318        
     359
    319360        # Subtract geo ref from points
    320361        if not self.is_absolute():
    … …  
    323364            points[:,1] -= self.yllcorner
    324365
    325        
    326366        if is_list:
    327367            points = points.tolist()
    … …  
    329369        return points
    330370
    331 
    332 
     371    ##
     372    # @brief ??
     373    # @param other ??
    333374    def reconcile_zones(self, other):
    334375        if other is None:
    335376            other = Geo_reference()
    336         if self.zone == other.zone or\
    337                self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE:
     377        if (self.zone == other.zone or
     378            self.zone == DEFAULT_ZONE and
     379            other.zone == DEFAULT_ZONE):
    338380            pass
    339381        elif self.zone == DEFAULT_ZONE:
    340382            self.zone = other.zone
    341383        elif other.zone == DEFAULT_ZONE:
    342             other.zone = self.zone           
    343         else:   
    344             msg = 'Geospatial data must be in the same '+\
    345                   'ZONE to allow reconciliation. I got zone %d and %d'\
    346                   %(self.zone, other.zone)
     384            other.zone = self.zone
     385        else:
     386            msg = ('Geospatial data must be in the same '
     387                   'ZONE to allow reconciliation. I got zone %d and %d'
     388                   % (self.zone, other.zone))
    347389            raise ANUGAError, msg
    348    
     390
    349391    #def easting_northing2geo_reffed_point(self, x, y):
    350392    #    return [x-self.xllcorner, y - self.xllcorner]
    … …  
    353395    #    return [x-self.xllcorner, y - self.xllcorner]
    354396
     397    ##
     398    # @brief Get origin of this geo_reference.
     399    # @return (zone, xllcorner, yllcorner).
    355400    def get_origin(self):
    356401        return (self.zone, self.xllcorner, self.yllcorner)
    357    
     402
     403    ##
     404    # @brief Get a string representation of this geo_reference instance.
    358405    def __repr__(self):
    359         return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner)
    360    
    361     def __cmp__(self,other):
     406        return ('(zone=%i easting=%f, northing=%f)'
     407                % (self.zone, self.xllcorner, self.yllcorner))
     408
     409    ##
     410    # @brief Compare two geo_reference instances.
     411    # @param self This geo_reference instance.
     412    # @param other Another geo_reference instance to compare against.
     413    # @return 0 if instances have the same attributes, else 1.
     414    # @note Attributes are: zone, xllcorner, yllcorner.
     415    def __cmp__(self, other):
    362416        # FIXME (DSG) add a tolerence
    363417        if other is None:
    … …  
    372426        return cmp
    373427
     428
     429##
     430# @brief Write a geo_reference to a NetCDF file (usually SWW).
     431# @param origin A georef instance or parameters to create a georef instance.
     432# @param outfile Path to file to write.
     433# @return A normalized geo_reference.
    374434def write_NetCDF_georeference(origin, outfile):
    375     """
    376     Write georeferrence info to a netcdf file, usually sww.
    377 
    378     The origin can be a georef instance or parrameters for a geo_ref instance
     435    """Write georeference info to a netcdf file, usually sww.
     436
     437    The origin can be a georef instance or parameters for a geo_ref instance
    379438
    380439    outfile is the name of the file to be written to.
    381440    """
     441
    382442    geo_ref = ensure_geo_reference(origin)
    383443    geo_ref.write_NetCDF(outfile)
    384444    return geo_ref
    385445
     446
     447##
     448# @brief Convert an object to a georeference instance.
     449# @param origin A georef instance or (zone, xllcorner, yllcorner)
     450# @return A georef object, or None if 'origin' was None.
    386451def ensure_geo_reference(origin):
    387452    """
    … …  
    392457    effect code logic
    393458    """
    394     if isinstance(origin, Geo_reference):
     459
     460    if isinstance(origin, Geo_reference):
    395461        geo_ref = origin
    396462    elif origin is None:
    397463        geo_ref = None
    398464    else:
    399         geo_ref = apply(Geo_reference, origin)       
     465        geo_ref = apply(Geo_reference, origin)
     466
    400467    return geo_ref
    401468
    402    
     469
    403470#-----------------------------------------------------------------------
    404471

  • anuga_core/source/anuga/coordinate_transforms/projection_test_points.lic

    r6521 r7276  
    1 <?xml version="1.0" encoding="iso-8859-1"?>
    2 
     1<?xml version='1.0' encoding='iso-8859-1'?>
    32<ga_license_file>
    43  <metadata>
    … …  
    87  <datafile>
    98    <filename>projection_test_points.csv</filename>
    10     <checksum>-1184955128</checksum>
     9    <checksum>3110012168</checksum>
    1110    <publishable>Yes</publishable>
    1211    <accountable>Ole Nielsen</accountable>

  • anuga_core/source/anuga/coordinate_transforms/projection_test_points_z53.lic

    r6521 r7276  
    1 <?xml version="1.0" encoding="iso-8859-1"?>
    2 
     1<?xml version='1.0' encoding='iso-8859-1'?>
    32<ga_license_file>
    43  <metadata>

  • anuga_core/source/anuga/coordinate_transforms/projection_test_points_z54.lic

    r6521 r7276  
    1 <?xml version="1.0" encoding="iso-8859-1"?>
    2 
     1<?xml version='1.0' encoding='iso-8859-1'?>
    32<ga_license_file>
    43  <metadata>
    … …  
    87  <datafile>
    98    <filename>projection_test_points_z54.csv</filename>
    10     <checksum>-2053861105</checksum>
     9    <checksum>2241106191</checksum>
    1110    <publishable>Yes</publishable>
    1211    <accountable>Ole Nielsen</accountable>

  • anuga_core/source/anuga/coordinate_transforms/redfearn.lic

    r5056 r7276  
    1 <?xml version="1.0" encoding="iso-8859-1"?>
    2 
     1<?xml version='1.0' encoding='iso-8859-1'?>
    32<ga_license_file>
    43  <metadata>

  • anuga_core/source/anuga/coordinate_transforms/test_geo_reference.py

    r7063 r7276  
    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(isinstance(new_points, list), '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(isinstance(new_points, list), '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 = 12.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(isinstance(new_points, list), 'failed')
     268        self.failUnless(type(new_points) == type(points), 'failed')
     269        self.failUnless(new_points == expected_new_points, 'failed')
     270
     271        # and repeat from 'new_points = g.get_absolute(points)' above
     272        # to see if second call with same input gives same results.
     273        new_points = g.get_absolute(points)
     274
     275        self.failUnless(isinstance(new_points, list), 'failed')
     276        self.failUnless(type(new_points) == type(points), 'failed')
     277        self.failUnless(new_points == expected_new_points, 'failed')
     278
     279    def test_get_absolute_array(self):
     280        '''Same test as test_get_absolute_list(), but with numeric arrays.'''
     281
     282        # test with supplied offsets
     283        x = 7.0
     284        y = 3.0
     285       
     286        g = Geo_reference(56, x, y)
     287        points = num.array([[3.0,34.0], [64.0,6.0]])
     288        new_points = g.get_absolute(points)
     289
     290        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     291        self.failUnless(type(new_points) == type(points), 'failed')
     292        msg = 'points=\n%s\nnew_points=\n%s' % (str(points), str(new_points))
     293        for point, new_point in map(None, points, new_points):
     294            self.failUnless(point[0]+x == new_point[0], msg)
     295            self.failUnless(point[1]+y == new_point[1], msg)
     296
     297        # test with no supplied offsets
    250298        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')
     299        points = num.array([[3.0,34.0], [64.0,6.0]])
     300        new_points = g.get_absolute(points)
     301
     302        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     303        self.failUnless(type(new_points) == type(points), 'failed')
     304        self.failUnless(num.alltrue(points == new_points), 'failed')
     305
     306        # test that calling get_absolute twice does the right thing
     307        # first call
     308        dx = 11.0
     309        dy = 13.0
     310        g = Geo_reference(56, dx, dy)
     311        points = num.array([[3.0,34.0], [64.0,6.0]])
     312        expected_new_points = num.array([[3.0+dx,34.0+dy], [64.0+dx,6.0+dy]])
     313        new_points = g.get_absolute(points)
     314
     315        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     316        self.failUnless(type(new_points) == type(points), 'failed')
     317        msg = ('First call of .get_absolute() returned %s\nexpected %s'
     318               % (str(new_points), str(expected_new_points)))
     319        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
     320
     321        # and repeat from 'new_points = g.get_absolute(points)' above
     322        # to see if second call with same input gives same results.
     323        new_points = g.get_absolute(points)
     324
     325        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     326        self.failUnless(type(new_points) == type(points), 'failed')
     327        msg = ('Second call of .get_absolute() returned\n%s\nexpected\n%s'
     328               % (str(new_points), str(expected_new_points)))
     329        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
     330
     331        # and repeat again to see if *third* call with same input
     332        # gives same results.
     333        new_points = g.get_absolute(points)
     334
     335        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     336        self.failUnless(type(new_points) == type(points), 'failed')
     337        msg = ('Third call of .get_absolute() returned %s\nexpected %s'
     338               % (str(new_points), str(expected_new_points)))
     339        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
     340
     341    def test_get_relative_list(self):
     342        # test with supplied offsets
     343        x = 7.0
     344        y = 3.0
     345       
     346        g = Geo_reference(56, x, y)
     347        points = [[3.0,34.0], [64.0,6.0]]
     348        new_points = g.get_relative(points)
     349
     350        self.failUnless(isinstance(new_points, list), 'failed')
     351        self.failUnless(type(new_points) == type(points), 'failed')
     352        for point, new_point in map(None, points, new_points):
     353            self.failUnless(point[0]-x == new_point[0], 'failed')
     354            self.failUnless(point[1]-y == new_point[1], 'failed')
     355
     356        # test with no supplied offsets
     357        g = Geo_reference()
     358        points = [[3.0,34.0], [64.0,6.0]]
     359        new_points = g.get_relative(points)
     360
     361        self.failUnless(isinstance(new_points, list), 'failed')
     362        self.failUnless(type(new_points) == type(points), 'failed')
     363        for point, new_point in map(None, points, new_points):
     364            self.failUnless(point[0] == new_point[0], 'failed')
     365            self.failUnless(point[1] == new_point[1], 'failed')
    261366           
     367        # test that calling get_absolute twice does the right thing
     368        # first call
     369        dx = 10.0
     370        dy = 12.0
     371        g = Geo_reference(56, dx, dy)
     372        points = [[3.0,34.0], [64.0,6.0]]
     373        expected_new_points = [[3.0-dx,34.0-dy], [64.0-dx,6.0-dy]]
     374        new_points = g.get_relative(points)
     375
     376        self.failUnless(isinstance(new_points, list), 'failed')
     377        self.failUnless(type(new_points) == type(points), 'failed')
     378        self.failUnless(new_points == expected_new_points, 'failed')
     379
     380        # and repeat from 'new_points = g.get_absolute(points)' above
     381        # to see if second call with same input gives same results.
     382        new_points = g.get_relative(points)
     383
     384        self.failUnless(isinstance(new_points, list), 'failed')
     385        self.failUnless(type(new_points) == type(points), 'failed')
     386        self.failUnless(new_points == expected_new_points, 'failed')
     387
     388    def test_get_relative_array(self):
     389        '''Same test as test_get_relative_list(), but with numeric arrays.'''
     390
     391        # test with supplied offsets
     392        x = 7.0
     393        y = 3.0
     394       
     395        g = Geo_reference(56, x, y)
     396        points = num.array([[3.0,34.0], [64.0,6.0]])
     397        new_points = g.get_relative(points)
     398
     399        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     400        self.failUnless(type(new_points) == type(points), 'failed')
     401        msg = 'points=\n%s\nnew_points=\n%s' % (str(points), str(new_points))
     402        for point, new_point in map(None, points, new_points):
     403            self.failUnless(point[0]-x == new_point[0], msg)
     404            self.failUnless(point[1]-y == new_point[1], msg)
     405
     406        # test with no supplied offsets
     407        g = Geo_reference()
     408        points = num.array([[3.0,34.0], [64.0,6.0]])
     409        new_points = g.get_relative(points)
     410
     411        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     412        self.failUnless(type(new_points) == type(points), 'failed')
     413        self.failUnless(num.alltrue(points == new_points), 'failed')
     414
     415        # test that calling get_relative twice does the right thing
     416        # first call
     417        dx = 11.0
     418        dy = 13.0
     419        g = Geo_reference(56, dx, dy)
     420        points = num.array([[3.0,34.0], [64.0,6.0]])
     421        expected_new_points = num.array([[3.0-dx,34.0-dy], [64.0-dx,6.0-dy]])
     422        new_points = g.get_relative(points)
     423
     424        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     425        self.failUnless(type(new_points) == type(points), 'failed')
     426        msg = ('First call of .get_relative() returned %s\nexpected %s'
     427               % (str(new_points), str(expected_new_points)))
     428        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
     429
     430        # and repeat from 'new_points = g.get_relative(points)' above
     431        # to see if second call with same input gives same results.
     432        new_points = g.get_relative(points)
     433
     434        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     435        self.failUnless(type(new_points) == type(points), 'failed')
     436        msg = ('Second call of .get_relative() returned\n%s\nexpected\n%s'
     437               % (str(new_points), str(expected_new_points)))
     438        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
     439
     440        # and repeat again to see if *third* call with same input
     441        # gives same results.
     442        new_points = g.get_relative(points)
     443
     444        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
     445        self.failUnless(type(new_points) == type(points), 'failed')
     446        msg = ('Third call of .get_relative() returned %s\nexpected %s'
     447               % (str(new_points), str(expected_new_points)))
     448        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
     449
    262450    def test_is_absolute(self):
    263451       
    … …  
    440628
    441629        # points in num.array()
    442         points = num.array(((2,3), (3,1), (5,2)), num.Float)
     630        points = num.array(((2,3), (3,1), (5,2)), num.float)
    443631        abs_points = geo.get_absolute(points)
    444632        # check we haven't changed 'points' itself
    … …  
    482670        file_name = tempfile.mktemp(".geo_referenceTest")
    483671       
    484         out_file = NetCDFFile(file_name, 'w')
     672        out_file = NetCDFFile(file_name, netcdf_mode_w)
    485673        g.write_NetCDF(out_file)
    486674        out_file.close()
    487675       
    488         in_file = NetCDFFile(file_name, 'r')
     676        in_file = NetCDFFile(file_name, netcdf_mode_r)
    489677        new_g = Geo_reference(NetCDFObject=in_file)
    490678        in_file.close()
    … …  
    537725       
    538726#-------------------------------------------------------------
     727
    539728if __name__ == "__main__":
    540 
    541     suite = unittest.makeSuite(geo_referenceTestCase,'test')
     729    suite = unittest.makeSuite(geo_referenceTestCase, 'test')
    542730    runner = unittest.TextTestRunner() #verbosity=2)
    543731    runner.run(suite)

  • anuga_core/source/anuga/coordinate_transforms/test_lat_long_UTM_conversion.py

    r6149 r7276  
    1212from anuga.utilities.anuga_exceptions import ANUGAError
    1313
    14 import Numeric as num
     14import numpy as num
    1515
    1616
    … …  
    121121        assert num.allclose(lat,  lat_calced)
    122122        assert num.allclose(lon, long_calced)
     123
    123124#-------------------------------------------------------------
     125
    124126if __name__ == "__main__":
    125 
    126     #mysuite = unittest.makeSuite(TestCase,'test_convert_latlon_to_UTM1')
    127127    mysuite = unittest.makeSuite(TestCase,'test')
    128128    runner = unittest.TextTestRunner()

  • anuga_core/source/anuga/coordinate_transforms/test_point.py

    r2253 r7276  
    117117                        %(self.RSISE.BearingTo(self.Kobenhavn), B))
    118118
     119#-------------------------------------------------------------
    119120
    120 
    121 #-------------------------------------------------------------
    122121if __name__ == "__main__":
    123 
    124122    mysuite = unittest.makeSuite(TestCase,'test')
    125123    runner = unittest.TextTestRunner()
    126124    runner.run(mysuite)
    127125
    128 
    129 
    130 
    131 
    132 
    133 
    134 

  • anuga_core/source/anuga/coordinate_transforms/test_redfearn.py

    r6595 r7276  
    1010from redfearn import *
    1111from anuga.utilities.anuga_exceptions import ANUGAError
     12
    1213from anuga.utilities.system_tools import get_pathname_from_package
    1314from os.path import join
    14 import Numeric as num
     15import numpy as num
    1516
    1617
    … …  
    139140        assert num.allclose(northing, 6762559.15978)
    140141
    141         # Testing using the native zone
     142        #Testing using the native zone
    142143        zone, easting, northing = redfearn(-29.233299999,114.05, zone=50)
    143144
    … …  
    146147        assert num.allclose(northing, 6762559.15978)
    147148
    148         # Then project to zone 49
    149         zone, easting, northing = redfearn(-29.233299999,114.05, zone=49)
     149        #Then project to zone 49
     150        zone, easting, northing = redfearn(-29.233299999,114.05,zone=49)
    150151
    151152        assert zone == 49
    … …  
    155156       
    156157
    157         # First test native projection (zone 49)
     158       
     159
     160        #First test native projection (zone 49)
    158161        zone, easting, northing = redfearn(-29.1333,113.9667)
    159162
    … …  
    162165        assert num.allclose(northing, 6773605.46384)
    163166
    164         # Then project to zone 50
    165         zone, easting, northing = redfearn(-29.1333,113.9667, zone=50)
     167        #Then project to zone 50
     168        zone, easting, northing = redfearn(-29.1333,113.9667,zone=50)
    166169
    167170        assert zone == 50
    … …  
    169172        assert num.allclose(northing, 6773440.04726)
    170173
    171         # Testing point on zone boundary
    172         # First test native projection (zone 50)
     174        #Testing point on zone boundary
     175        #First test native projection (zone 50)
    173176        zone, easting, northing = redfearn(-29.1667,114)
    174177
    … …  
    177180        assert num.allclose(northing, 6769820.01453)
    178181
    179         # Then project to zone 49
    180         zone, easting, northing = redfearn(-29.1667,114, zone=49)
     182        #Then project to zone 49
     183        zone, easting, northing = redfearn(-29.1667,114,zone=49)
    181184
    182185        assert zone == 49
    … …  
    184187        assert num.allclose(northing, 6769820.01453)
    185188
    186         # Testing furthest point in Geraldton scenario)
    187         # First test native projection (zone 49)
     189        #Testing furthest point in Geraldton scenario)
     190        #First test native projection (zone 49)
    188191        zone, easting, northing = redfearn(-28.2167,113.4167)
    189192
    … …  
    192195        assert num.allclose(northing, 6876426.38578)
    193196
    194         # Then project to zone 50
     197        #Then project to zone 50
    195198        zone, easting, northing = redfearn(-28.2167,113.4167,zone=50)
    196199
    … …  
    199202        assert num.allclose(northing, 6873587.50926)
    200203
    201         # Testing outside GDA zone (New Zeland)
    202         # First test native projection (zone 60)
     204        #Testing outside GDA zone (New Zeland)
     205        #First test native projection (zone 60)
    203206        zone, easting, northing = redfearn(-44,178)
    204207
    … …  
    207210        assert num.allclose(northing, 5127641.114461)
    208211
    209         # Then project to zone 59
     212        #Then project to zone 59
    210213        zone, easting, northing = redfearn(-44,178,zone=59)
    211214
    … …  
    214217        assert num.allclose(northing, 5104249.395469)
    215218
    216         # Then skip three zones 57 (native 60)
     219        #Then skip three zones 57 (native 60)
    217220        zone, easting, northing = redfearn(-44,178,zone=57)
    218221
    … …  
    229232        # Google Earth (he he)
    230233       
    231         # Testing outside GDA zone (Northern Hemisphere)
    232         # First test native projection (zone 57)
     234        #Testing outside GDA zone (Northern Hemisphere)
     235        #First test native projection (zone 57)
    233236        zone, easting, northing = redfearn(44,156)
    234237
    … …  
    243246        #assert num.allclose(northing, 14876249.1268)
    244247       
    245         # Then project to zone 56
     248        #Then project to zone 56
    246249        zone, easting, northing = redfearn(44,156,zone=56)
    247250
    … …  
    273276    #    #assert allclose(northing, 6181725.1724276)
    274277
    275 
    276278   
    277279    def test_nonstandard_meridian_coinciding_with_native(self):
    … …  
    505507           
    506508#-------------------------------------------------------------
     509
    507510if __name__ == "__main__":
    508 
    509511    mysuite = unittest.makeSuite(TestCase,'test')
    510512    runner = unittest.TextTestRunner()
Note: See TracChangeset for help on using the changeset viewer.