Ignore:
Timestamp:
Feb 18, 2009, 2:44:41 PM (16 years ago)
Author:
rwilson
Message:

Ongoing conversion changes.

Location:
branches/numpy/anuga/coordinate_transforms
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r6304 r6360  
    99import types, sys
    1010from anuga.utilities.numerical_tools import ensure_numeric
    11 from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, ParsingError, \
    12      ShapeError
     11from anuga.utilities.anuga_exceptions import ANUGAError, TitleError, \
     12                                             ParsingError, ShapeError
    1313from anuga.config import netcdf_float, netcdf_int, netcdf_float32
    1414
    … …  
    1717
    1818DEFAULT_ZONE = -1
    19 TITLE = '#geo reference' + "\n" #this title is referred to in the test format
     19TITLE = '#geo reference' + "\n" # this title is referred to in the test format
    2020
    2121DEFAULT_PROJECTION = 'UTM'
    … …  
    2323DEFAULT_UNITS = 'm'
    2424DEFAULT_FALSE_EASTING = 500000
    25 DEFAULT_FALSE_NORTHING = 10000000 #Default for southern hemisphere
    26 
     25DEFAULT_FALSE_NORTHING = 10000000    # Default for southern hemisphere
     26
     27
     28##
     29# @brief A class for ...
    2730class Geo_reference:
    2831    """
     32    Attributes of the Geo_reference class:
     33        .zone           The UTM zone (default is -1)
     34        .false_easting  ??
     35        .false_northing ??
     36        .datum          The Datum used (default is wgs84)
     37        .projection     The projection used (default is 'UTM')
     38        .units          The units of measure used (default metres)
     39        .xllcorner      The X coord of origin (default is 0.0 wrt UTM grid)
     40        .yllcorner      The y coord of origin (default is 0.0 wrt UTM grid)
     41        .is_absolute    ??
     42
    2943    """
    3044
     45    ##
     46    # @brief Instantiate an instance of class Geo_reference.
     47    # @param zone The UTM zone.
     48    # @param xllcorner X coord of origin of georef.
     49    # @param yllcorner Y coord of origin of georef.
     50    # @param datum ??
     51    # @param projection The projection used (default UTM).
     52    # @param units Units used in measuring distance (default m).
     53    # @param false_easting ??
     54    # @param false_northing ??
     55    # @param NetCDFObject NetCDF file *handle* to write to.
     56    # @param ASCIIFile ASCII text file *handle* to write to.
     57    # @param read_title Title of the georeference text.
    3158    def __init__(self,
    32                  zone = DEFAULT_ZONE,
    33                  xllcorner = 0.0,
    34                  yllcorner = 0.0,
    35                  datum = DEFAULT_DATUM,
    36                  projection = DEFAULT_PROJECTION,
    37                  units = DEFAULT_UNITS,
    38                  false_easting = DEFAULT_FALSE_EASTING,
    39                  false_northing = DEFAULT_FALSE_NORTHING,
     59                 zone=DEFAULT_ZONE,
     60                 xllcorner=0.0,
     61                 yllcorner=0.0,
     62                 datum=DEFAULT_DATUM,
     63                 projection=DEFAULT_PROJECTION,
     64                 units=DEFAULT_UNITS,
     65                 false_easting=DEFAULT_FALSE_EASTING,
     66                 false_northing=DEFAULT_FALSE_NORTHING,
    4067                 NetCDFObject=None,
    4168                 ASCIIFile=None,
    … …  
    4370        """
    4471        input:
    45         NetCDFObject - a handle to the netCDF file to be written to 
     72        NetCDFObject - a handle to the netCDF file to be written to
    4673        ASCIIFile - a handle to the text file
    4774        read_title - the title of the georeference text, if it was read in.
    … …  
    5481         must be the default info, since ANUGA assumes it isn't
    5582         changing.
    56          """
     83        """
     84
    5785        if zone is None:
    5886            zone = DEFAULT_ZONE
    5987        self.false_easting = false_easting
    60         self.false_northing = false_northing       
     88        self.false_northing = false_northing
    6189        self.datum = datum
    6290        self.projection = projection
    63         self.zone = zone       
     91        self.zone = zone
    6492        self.units = units
    6593        self.xllcorner = xllcorner
    66         self.yllcorner = yllcorner       
    67         #self.is_absolute = num.allclose([self.xllcorner, self.yllcorner], 0) 
    68            
     94        self.yllcorner = yllcorner
     95        #self.is_absolute = num.allclose([self.xllcorner, self.yllcorner], 0)
     96
    6997        if NetCDFObject is not None:
    7098            self.read_NetCDF(NetCDFObject)
    71  
     99
    72100        if ASCIIFile is not None:
    73101            self.read_ASCII(ASCIIFile, read_title=read_title)
    74102
     103#    # Might be better to have this method instead of the following 3.
     104#    def get_origin(self):
     105#        return (self.zone, self.xllcorner, self.yllcorner)
     106
     107    ##
     108    # @brief Get the X coordinate of the origin of this georef.
    75109    def get_xllcorner(self):
    76110        return self.xllcorner
    77    
     111
     112    ##
     113    # @brief Get the Y coordinate of the origin of this georef.
    78114    def get_yllcorner(self):
    79115        return self.yllcorner
    80    
     116
     117    ##
     118    # @brief Get the zone of this georef.
    81119    def get_zone(self):
    82120        return self.zone
    83        
     121
     122    ##
     123    # @brief Write <something> to an open NetCDF file.
     124    # @param outfile Handle to open NetCDF file.
    84125    def write_NetCDF(self, outfile):
    85         outfile.xllcorner = self.xllcorner 
     126        outfile.xllcorner = self.xllcorner
    86127        outfile.yllcorner = self.yllcorner
    87128        outfile.zone = self.zone
    … …  
    90131        outfile.false_northing = self.false_northing
    91132
    92         outfile.datum = self.datum       
     133        outfile.datum = self.datum
    93134        outfile.projection = self.projection
    94135        outfile.units = self.units
    95136
     137    ##
     138    # @brief Read data from an open NetCDF file.
     139    # @param infile Handle to open NetCDF file.
    96140    def read_NetCDF(self, infile):
    97141        self.xllcorner = infile.xllcorner[0]
    98         self.yllcorner = infile.yllcorner[0] 
     142        self.yllcorner = infile.yllcorner[0]
    99143        self.zone = infile.zone[0]
    100144
    101        
    102145        # Fix some assertion failures
    103146        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
    104147            self.zone = self.zone[0]
    105         if isinstance(self.xllcorner, num.ndarray) and self.xllcorner.shape == ():
     148        if (isinstance(self.xllcorner, num.ndarray) and
     149                self.xllcorner.shape == ()):
    106150            self.xllcorner = self.xllcorner[0]
    107         if isinstance(self.yllcorner, num.ndarray) and self.yllcorner.shape == ():
     151        if (isinstance(self.yllcorner, num.ndarray) and
     152                self.yllcorner.shape == ()):
    108153            self.yllcorner = self.yllcorner[0]
    109154
    … …  
    113158                self.yllcorner.dtype.kind in num.typecodes['Integer'])
    114159        assert (self.zone.dtype.kind in num.typecodes['Integer'])
    115        
     160
    116161        try:
    117162            self.false_easting = infile.false_easting[0]
    118163            self.false_northing = infile.false_northing[0]
    119        
    120             self.datum = infile.datum       
     164
     165            self.datum = infile.datum
    121166            self.projection = infile.projection
    122167            self.units = infile.units
    123168        except:
    124169            pass
    125         if (self.false_easting != DEFAULT_FALSE_EASTING):
    126             print "WARNING: False easting of %f specified." %self.false_easting
    127             print "Default false easting is %f." %DEFAULT_FALSE_EASTING
     170
     171        if self.false_easting != DEFAULT_FALSE_EASTING:
     172            print "WARNING: False easting of %f specified." % self.false_easting
     173            print "Default false easting is %f." % DEFAULT_FALSE_EASTING
    128174            print "ANUGA does not correct for differences in False Eastings."
    129            
    130         if (self.false_northing != DEFAULT_FALSE_NORTHING):
    131             print "WARNING: False northing of %f specified." \
    132                   %self.false_northing
    133             print "Default false northing is %f." %DEFAULT_FALSE_NORTHING
     175
     176        if self.false_northing != DEFAULT_FALSE_NORTHING:
     177            print ("WARNING: False northing of %f specified."
     178                   % self.false_northing)
     179            print "Default false northing is %f." % DEFAULT_FALSE_NORTHING
    134180            print "ANUGA does not correct for differences in False Northings."
    135            
    136         if (self.datum.upper() != DEFAULT_DATUM.upper()):
    137             print "WARNING: Datum of %s specified." \
    138                   %self.datum
    139             print "Default Datum is %s." %DEFAULT_DATUM
     181
     182        if self.datum.upper() != DEFAULT_DATUM.upper():
     183            print "WARNING: Datum of %s specified." % self.datum
     184            print "Default Datum is %s." % DEFAULT_DATUM
    140185            print "ANUGA does not correct for differences in datums."
    141            
    142         if (self.projection.upper() != DEFAULT_PROJECTION.upper()):
    143             print "WARNING: Projection of %s specified." \
    144                   %self.projection
    145             print "Default Projection is %s." %DEFAULT_PROJECTION
     186
     187        if self.projection.upper() != DEFAULT_PROJECTION.upper():
     188            print "WARNING: Projection of %s specified." % self.projection
     189            print "Default Projection is %s." % DEFAULT_PROJECTION
    146190            print "ANUGA does not correct for differences in Projection."
    147            
    148         if (self.units.upper() != DEFAULT_UNITS.upper()):
    149             print "WARNING: Units of %s specified." \
    150                   %self.units
    151             print "Default units is %s." %DEFAULT_UNITS
     191
     192        if self.units.upper() != DEFAULT_UNITS.upper():
     193            print "WARNING: Units of %s specified." % self.units
     194            print "Default units is %s." % DEFAULT_UNITS
    152195            print "ANUGA does not correct for differences in units."
    153        
    154     ### ASCII files with geo-refs are currently not used   
     196
     197################################################################################
     198# ASCII files with geo-refs are currently not used
     199################################################################################
     200
     201    ##
     202    # @brief Write georef data to an open text file.
     203    # @param fd Handle to open text file.
    155204    def write_ASCII(self, fd):
    156205        fd.write(TITLE)
    157         fd.write(str(self.zone) + "\n") 
    158         fd.write(str(self.xllcorner) + "\n") 
     206        fd.write(str(self.zone) + "\n")
     207        fd.write(str(self.xllcorner) + "\n")
    159208        fd.write(str(self.yllcorner) + "\n")
    160209
    161     def read_ASCII(self, fd,read_title=None):
     210    ##
     211    # @brief Read georef data from an open text file.
     212    # @param fd Handle to open text file.
     213    def read_ASCII(self, fd, read_title=None):
    162214        try:
    163215            if read_title == None:
    164                 read_title = fd.readline() # remove the title line
     216                read_title = fd.readline()     # remove the title line
    165217            if read_title[0:2].upper() != TITLE[0:2].upper():
    166                 msg = 'File error.  Expecting line: %s.  Got this line: %s' \
    167                       %(TITLE, read_title)
    168                 raise TitleError, msg 
     218                msg = ('File error.  Expecting line: %s.  Got this line: %s'
     219                       % (TITLE, read_title))
     220                raise TitleError, msg
    169221            self.zone = int(fd.readline())
    170222            self.xllcorner = float(fd.readline())
    171223            self.yllcorner = float(fd.readline())
    172224        except SyntaxError:
    173                 msg = 'File error.  Got syntax error while parsing geo reference'
    174                 raise ParsingError, msg
    175            
     225            msg = 'File error.  Got syntax error while parsing geo reference'
     226            raise ParsingError, msg
     227
    176228        # Fix some assertion failures
    177229        if isinstance(self.zone, num.ndarray) and self.zone.shape == ():
    178230            self.zone = self.zone[0]
    179         if isinstance(self.xllcorner, num.ndarray) and self.xllcorner.shape == ():
     231        if (isinstance(self.xllcorner, num.ndarray) and
     232                self.xllcorner.shape == ()):
    180233            self.xllcorner = self.xllcorner[0]
    181         if isinstance(self.yllcorner, num.ndarray) and self.yllcorner.shape == ():
     234        if (isinstance(self.yllcorner, num.ndarray) and
     235                self.yllcorner.shape == ()):
    182236            self.yllcorner = self.yllcorner[0]
    183237
    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 
    189        
    190     def change_points_geo_ref(self, points,
    191                               points_geo_ref=None):
     238################################################################################
     239
     240    def change_points_geo_ref(self, points, points_geo_ref=None):
    192241        """
    193242        Change the geo reference of a list or numeric array of points to
    … …  
    201250
    202251        points = ensure_numeric(points, num.float)
    203        
     252
    204253        if len(points.shape) == 1:
    205             #One point has been passed
     254            # One point has been passed
    206255            msg = 'Single point must have two elements'
    207256            assert len(points) == 2, msg
    208257            points = num.reshape(points, (1,2))
    209258
    210         msg = 'Points array must be two dimensional.\n'
    211         msg += 'I got %d dimensions' %len(points.shape)
     259        msg = ('Points array must be two dimensional.\n'
     260               'I got %d dimensions' % len(points.shape))
    212261        assert len(points.shape) == 2, msg
    213262
    214         msg = 'Input must be an N x 2 array or list of (x,y) values. '
    215         msg += 'I got an %d x %d array' %points.shape   
    216         assert points.shape[1] == 2, msg               
    217 
    218        
    219         # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
    220         # are identical in the two geo refs.   
     263        msg = ('Input must be an N x 2 array or list of (x,y) values. '
     264               'I got an %d x %d array' % points.shape)
     265        assert points.shape[1] == 2, msg
     266
     267        # FIXME (Ole): Could also check if zone, xllcorner, yllcorner
     268        # are identical in the two geo refs.
    221269        if points_geo_ref is not self:
    222270            # If georeferences are different
    223        
    224271            if not points_geo_ref is None:
    225272                # Convert points to absolute coordinates
    226                 points[:,0] += points_geo_ref.xllcorner 
    227                 points[:,1] += points_geo_ref.yllcorner 
    228        
     273                points[:,0] += points_geo_ref.xllcorner
     274                points[:,1] += points_geo_ref.yllcorner
     275
    229276            # Make points relative to primary geo reference
    230             points[:,0] -= self.xllcorner 
     277            points[:,0] -= self.xllcorner
    231278            points[:,1] -= self.yllcorner
    232279
    233            
     280        # if we got a list, return a list
    234281        if is_list:
    235282            points = points.tolist()
    236            
     283
    237284        return points
    238285
    239    
     286    ##
     287    # @brief Is this georef absolute?
     288    # @return True if ???
    240289    def is_absolute(self):
    241         """Return True if xllcorner==yllcorner==0 indicating that points
    242         in question are absolute.
    243         """
    244 
    245         return num.allclose([self.xllcorner, self.yllcorner], 0)
    246 
    247        
    248     ##
    249     # @brief
    250     # @param points
    251     # @return
    252     # @note
     290        """Return True if xllcorner==yllcorner==0
     291        indicating that points in question are absolute.
     292        """
     293
     294        return num.allclose([self.xllcorner, self.yllcorner], 0)
     295
     296    ##
     297    # @brief Convert points to absolute points in this georef instance.
     298    # @param points
     299    # @return
     300    # @note
    253301    def get_absolute(self, points):
    254302        """Given a set of points geo referenced to this instance,
    … …  
    259307        #    return points
    260308
     309        # remember if we got a list
    261310        is_list = False
    262311        if type(points) == types.ListType:
    … …  
    266315        if len(points.shape) == 1:
    267316            #One point has been passed
    268             msg = 'Single point must have two elements'
    269317            if not len(points) == 2:
    270                 raise ShapeError, msg   
     318                msg = 'Single point must have two elements'
     319                raise ShapeError, msg
    271320                #points = reshape(points, (1,2))
    272321
    273         msg = 'Input must be an N x 2 array or list of (x,y) values. '
    274         msg += 'I got an %d x %d array' %points.shape   
    275322        if not points.shape[1] == 2:
    276             raise ShapeError, msg   
    277        
     323            msg = ('Input must be an N x 2 array or list of (x,y) values. '
     324                   'I got an %d x %d array' % points.shape)
     325            raise ShapeError, msg
     326
    278327        # Add geo ref to points
    279328        #if not self.is_absolute:
    280329        if not self.is_absolute():
    281             points[:,0] += self.xllcorner 
     330            points[:,0] += self.xllcorner
    282331            points[:,1] += self.yllcorner
    283332            #self.is_absolute = True
    284        
     333
     334        # if we got a list, return a list
    285335        if is_list:
    286336            points = points.tolist()
    287              
     337
    288338        return points
    289339
    290 
     340    ##
     341    # @brief Convert points to relative measurement.
     342    # @param points Points to convert to relative measurements.
     343    # @return A set of points relative to the geo_reference instance.
    291344    def get_relative(self, points):
    292         """
    293         Given a set of points in absolute UTM coordinates,
     345        """Given a set of points in absolute UTM coordinates,
    294346        make them relative to this geo_reference instance,
    295347        return the points as relative values.
    … …  
    305357        if len(points.shape) == 1:
    306358            #One point has been passed
    307             msg = 'Single point must have two elements'
    308359            if not len(points) == 2:
    309                 raise ShapeError, msg   
     360                msg = 'Single point must have two elements'
     361                raise ShapeError, msg
    310362                #points = reshape(points, (1,2))
    311363
    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   
    315364        if not points.shape[1] == 2:
    316             raise ShapeError, msg   
    317            
    318        
     365            msg = ('Input must be an N x 2 array or list of (x,y) values. '
     366                   'I got an %d x %d array' % points.shape)
     367            raise ShapeError, msg
     368
    319369        # Subtract geo ref from points
    320370        if not self.is_absolute():
    321             points[:,0] -= self.xllcorner 
     371            points[:,0] -= self.xllcorner
    322372            points[:,1] -= self.yllcorner
    323373
    324        
     374        # if we got a list, return a list
    325375        if is_list:
    326376            points = points.tolist()
    327              
     377
    328378        return points
    329379
    330 
    331 
     380    ##
     381    # @brief ??
     382    # @param other ??
    332383    def reconcile_zones(self, other):
    333384        if other is None:
    334385            other = Geo_reference()
    335         if self.zone == other.zone or\
    336                self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE:
     386        if (self.zone == other.zone or
     387            self.zone == DEFAULT_ZONE and
     388            other.zone == DEFAULT_ZONE):
    337389            pass
    338390        elif self.zone == DEFAULT_ZONE:
    339391            self.zone = other.zone
    340392        elif other.zone == DEFAULT_ZONE:
    341             other.zone = self.zone           
    342         else:   
    343             msg = 'Geospatial data must be in the same '+\
    344                   'ZONE to allow reconciliation. I got zone %d and %d'\
    345                   %(self.zone, other.zone)
     393            other.zone = self.zone
     394        else:
     395            msg = ('Geospatial data must be in the same '
     396                   'ZONE to allow reconciliation. I got zone %d and %d'
     397                   % (self.zone, other.zone))
    346398            raise ANUGAError, msg
    347    
     399
    348400    #def easting_northing2geo_reffed_point(self, x, y):
    349401    #    return [x-self.xllcorner, y - self.xllcorner]
    … …  
    352404    #    return [x-self.xllcorner, y - self.xllcorner]
    353405
     406    ##
     407    # @brief Get origin of this geo_reference.
     408    # @return (zone, xllcorner, yllcorner).
    354409    def get_origin(self):
    355410        return (self.zone, self.xllcorner, self.yllcorner)
    356    
     411
     412    ##
     413    # @brief Get a string representation of this geo_reference instance.
    357414    def __repr__(self):
    358         return "(zone=%i easting=%f, northing=%f)" %(self.zone, self.xllcorner, self.yllcorner)
    359    
    360     def __cmp__(self,other):
     415        return ('(zone=%i easting=%f, northing=%f)'
     416                % (self.zone, self.xllcorner, self.yllcorner))
     417
     418    ##
     419    # @brief Compare two geo_reference instances.
     420    # @param self This geo_reference instance.
     421    # @param other Another geo_reference instance to compare against.
     422    # @return 0 if instances have the same attributes, else 1.
     423    # @note Attributes are: zone, xllcorner, yllcorner.
     424    def __cmp__(self, other):
    361425        # FIXME (DSG) add a tolerence
    362426        if other is None:
    … …  
    371435        return cmp
    372436
     437
     438##
     439# @brief Write a geo_reference to a NetCDF file (usually SWW).
     440# @param origin A georef instance or parameters to create a georef instance.
     441# @param outfile Path to file to write.
     442# @return A normalized geo_reference.
    373443def write_NetCDF_georeference(origin, outfile):
    374     """
    375     Write georeferrence info to a netcdf file, usually sww.
    376 
    377     The origin can be a georef instance or parrameters for a geo_ref instance
     444    """Write georeference info to a netcdf file, usually sww.
     445
     446    The origin can be a georef instance or parameters for a geo_ref instance
    378447
    379448    outfile is the name of the file to be written to.
    380449    """
     450
    381451    geo_ref = ensure_geo_reference(origin)
    382452    geo_ref.write_NetCDF(outfile)
    383453    return geo_ref
    384454
     455
     456##
     457# @brief Convert an object to a georeference instance.
     458# @param origin A georef instance or (zone, xllcorner, yllcorner)
     459# @return A georef object, or None if 'origin' was None.
    385460def ensure_geo_reference(origin):
    386461    """
    … …  
    391466    effect code logic
    392467    """
    393     if isinstance(origin, Geo_reference):
     468
     469    if isinstance(origin, Geo_reference):
    394470        geo_ref = origin
    395471    elif origin is None:
    396472        geo_ref = None
    397473    else:
    398         geo_ref = apply(Geo_reference, origin)       
     474        geo_ref = apply(Geo_reference, origin)
     475
    399476    return geo_ref
    400477
    401    
     478
    402479#-----------------------------------------------------------------------
    403480

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

    r6304 r6360  
    326326        self.failUnless(isinstance(new_points, num.ndarray), 'failed')
    327327        self.failUnless(type(new_points) == type(points), 'failed')
    328         msg = ('Second call of .get_absolute() returned %s\nexpected %s'
     328        msg = ('Second call of .get_absolute() returned\n%s\nexpected\n%s'
    329329               % (str(new_points), str(expected_new_points)))
    330330        self.failUnless(num.alltrue(expected_new_points == new_points), msg)
    … …  
    507507if __name__ == "__main__":
    508508
    509     #suite = unittest.makeSuite(geo_referenceTestCase,'test')
    510     suite = unittest.makeSuite(geo_referenceTestCase,'test_get_absolute')
     509    suite = unittest.makeSuite(geo_referenceTestCase,'test')
     510    #suite = unittest.makeSuite(geo_referenceTestCase,'test_get_absolute')
    511511    runner = unittest.TextTestRunner() #verbosity=2)
    512512    runner.run(suite)

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

    r2253 r6360  
    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 
Note: See TracChangeset for help on using the changeset viewer.