Changeset 875


Ignore:
Timestamp:
Feb 14, 2005, 8:42:23 AM (20 years ago)
Author:
ole
Message:

Cleaned up and made test a proper unit test.
Still thinking about what interface should look like??

Location:
inundation/ga/storm_surge/pyvolution/coordinate_transforms
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/coordinate_transforms/redfearn.py

    r868 r875  
     1"""Implementation of Redfearn's formula to compute UTM projections from latitude and longitude
    12
     3"""
    24
    35def degminsec2decimal_degrees(dd,mm,ss):
    … …  
    3840   
    3941
    40     #Convert to radians
    41     phi = lat*pi/180
    42     lam = lon*pi/180
    4342
    4443    #GDA Specifications
    45     a = 6378137.0         #Semi major axis
     44    a = 6378137.0                       #Semi major axis
    4645    inverse_flattening = 298.257222101  #1/f
    47     zone_width = 6              #Degrees
     46    K0 = 0.9996                         #Central scale factor   
     47    zone_width = 6                      #Degrees
     48
     49    longitude_of_central_meridian_zone0 = -183   
     50    longitude_of_western_edge_zone0 = -186   
    4851    false_easting = 500000
    49 
     52   
    5053    if lat < 0:
    5154        false_northing = 10000000  #Southern hemisphere
    5255    else:
    5356        false_northing = 0         #Northern hemisphere)
    54     longitude_of_central_meridian_zone1 = -177  #Needed???
    55     longitude_of_central_meridian_zone0 = -183   
    56     longitude_of_western_edge_zone0 = -186
    57     K0 = 0.9996                #Central scale factor
     57       
     58
     59
     60
    5861   
    5962    #Derived constants
    … …  
    7679    G = a*(1-n)*(1-n2)*(1+9*n2/4+225*n4/64)*pi/180
    7780
    78     sinphi = sin(phi)   #E9
     81
     82    phi = lat*pi/180     #Convert latitude to radians
     83
     84    sinphi = sin(phi)   
    7985    sin2phi = sin(2*phi)
    8086    sin4phi = sin(4*phi)
    8187    sin6phi = sin(6*phi)
    8288
    83     cosphi = cos(phi)   #E22
     89    cosphi = cos(phi)
    8490    cosphi2 = cosphi*cosphi
    8591    cosphi3 = cosphi*cosphi2
    … …  
    119125    m = term1 + term2 + term3 + term4 #OK
    120126
    121 
    122     zone_r = (lon - longitude_of_western_edge_zone0)/zone_width
    123     zone = int(zone_r)
     127    #Zone
     128    zone = int((lon - longitude_of_western_edge_zone0)/zone_width)
    124129    central_meridian = zone*zone_width+longitude_of_central_meridian_zone0
    125130
    … …  
    131136    omega6 = omega3*omega3
    132137    omega7 = omega*omega6
    133     omega8 = omega4*omega4   
     138    omega8 = omega4*omega4
     139   
    134140   
    135141    #Northing
    … …  
    147153    term3 = nu*cosphi5*(4*psi3*(1-6*t2)+psi2*(1+8*t2)-2*psi*t2+t4)*omega5/120
    148154    term4 = nu*cosphi7*(61-479*t2+179*t4-t6)*omega7/5040
    149     sum = term1 + term2 + term3 + term4
    150     easting = false_easting + K0*sum
     155    easting = false_easting + K0*(term1 + term2 + term3 + term4)
    151156   
    152157    return zone, easting, northing

  • inundation/ga/storm_surge/pyvolution/coordinate_transforms/test_redfearn.py

    r868 r875  
    11
    2 from redfearn import *
    3 from Numeric import allclose
    4 
     2#Test of redfearns formula. Tests can be verified at
     3#
    54#http://www.cellspark.com/UTM.html
    65#http://www.ga.gov.au/nmd/geodesy/datums/redfearn_geo_to_grid.jsp
    76
    87
    9 #TEST 1
    10 
    11 #latitude:  -37 39' 10.15610"
    12 #Longitude: 143 55' 35.38390"
    13 #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
    14 #Zone:   54   
    15 #Easting:  758173.797  Northing: 5828674.340
    16 #Latitude:   -37  39 ' 10.15610 ''  Longitude: 143  55 ' 35.38390 ''
    17 #Grid Convergence:  1  47 ' 19.36 ''  Point Scale: 1.00042107
    18 
    19 lat = degminsec2decimal_degrees(-37,39,10.15610)
    20 lon = degminsec2decimal_degrees(143,55,35.38390)
    21 assert allclose(lat, -37.65282114)
    22 assert allclose(lon, 143.9264955)
    23 
    24 dd,mm,ss = decimal_degrees2degminsec(-37.65282114)
    25 assert dd==-37
    26 assert mm==39
    27 assert allclose(ss, 10.15610)
    28 
    29 dd,mm,ss = decimal_degrees2degminsec(143.9264955)
    30 assert dd==143
    31 assert mm==55
    32 assert allclose(ss, 35.38390)
    33 
    34 #print dd,mm,ss
    35 
    36 zone, easting, northing = redfearn(lat,lon)
    37 
    38 assert zone == 54
    39 assert allclose(easting, 758173.797)
    40 assert allclose(northing, 5828674.340)
     8import unittest
     9from redfearn import *
     10from Numeric import allclose
    4111
    4212
     13#-------------------------------------------------------------
    4314
    44 #TEST 2
     15class TestCase(unittest.TestCase):
    4516
    46 #Latitude:  -37 57 03.7203
    47 #Longitude: 144 25 29.5244
    48 #Zone:   55   
    49 #Easting:  273741.297  Northing: 5796489.777
    50 #Latitude:   -37  57 ' 3.72030 ''  Longitude: 144  25 ' 29.52440 ''
    51 #Grid Convergence:  -1  35 ' 3.65 ''  Point Scale: 1.00023056
     17    def test_decimal_degrees_conversion(Self):
     18        lat = degminsec2decimal_degrees(-37,39,10.15610)
     19        lon = degminsec2decimal_degrees(143,55,35.38390)
     20        assert allclose(lat, -37.65282114)
     21        assert allclose(lon, 143.9264955)
    5222
    53 lat = degminsec2decimal_degrees(-37,57,03.7203)
    54 lon = degminsec2decimal_degrees(144,25,29.5244)
    55 print lat, lon
     23        dd,mm,ss = decimal_degrees2degminsec(-37.65282114)
     24        assert dd==-37
     25        assert mm==39
     26        assert allclose(ss, 10.15610)
    5627
    57 zone, easting, northing = redfearn(lat,lon)
    58 
    59 assert zone == 55
    60 assert allclose(easting, 273741.297)
    61 assert allclose(northing, 5796489.777)
     28        dd,mm,ss = decimal_degrees2degminsec(143.9264955)
     29        assert dd==143
     30        assert mm==55
     31        assert allclose(ss, 35.38390)
    6232
    6333
     34    def test_UTM_1(self):
     35        #latitude:  -37 39' 10.15610"
     36        #Longitude: 143 55' 35.38390"
     37        #Site Name:    GDA-MGA: (UTM with GRS80 ellipsoid)
     38        #Zone:   54   
     39        #Easting:  758173.797  Northing: 5828674.340
     40        #Latitude:   -37  39 ' 10.15610 ''  Longitude: 143  55 ' 35.38390 ''
     41        #Grid Convergence:  1  47 ' 19.36 ''  Point Scale: 1.00042107
    6442
    65 #Test 3
    66 lat = degminsec2decimal_degrees(-60,0,0)
    67 lon = degminsec2decimal_degrees(130,0,0)
    68 
    69 zone, easting, northing = redfearn(lat,lon)
    70 #print zone, easting, northing
    71 
    72 assert zone == 52
    73 assert allclose(easting, 555776.267)
    74 assert allclose(northing, 3348167.264)
     43        lat = degminsec2decimal_degrees(-37,39,10.15610)
     44        lon = degminsec2decimal_degrees(143,55,35.38390)
     45        assert allclose(lat, -37.65282114)
     46        assert allclose(lon, 143.9264955)
    7547
    7648
    77 #Test 4 (Kobenhavn)
    78 lat = 55.70248
    79 dd,mm,ss = decimal_degrees2degminsec(lat)
    80 print dd, mm, ss
     49        zone, easting, northing = redfearn(lat,lon)
    8150
    82 lon = 12.58364
    83 dd,mm,ss = decimal_degrees2degminsec(lon)
    84 print dd, mm, ss
    85 
    86 zone, easting, northing = redfearn(lat,lon)
    87 print zone, easting, northing
    88 
    89 assert zone == 33
    90 assert allclose(easting, 348157.631)
    91 #assert allclose(northing, 16175612.993) #With false_norting 10000000
    92 assert allclose(northing, 6175612.993)
     51        assert zone == 54
     52        assert allclose(easting, 758173.797)
     53        assert allclose(northing, 5828674.340)
    9354
    9455
     56    def test_UTM_2(self):
     57        #TEST 2
     58
     59        #Latitude:  -37 57 03.7203
     60        #Longitude: 144 25 29.5244
     61        #Zone:   55   
     62        #Easting:  273741.297  Northing: 5796489.777
     63        #Latitude:   -37  57 ' 3.72030 ''  Longitude: 144  25 ' 29.52440 ''
     64        #Grid Convergence:  -1  35 ' 3.65 ''  Point Scale: 1.00023056
     65
     66        lat = degminsec2decimal_degrees(-37,57,03.7203)
     67        lon = degminsec2decimal_degrees(144,25,29.5244)
     68        #print lat, lon
     69
     70        zone, easting, northing = redfearn(lat,lon)
     71
     72        assert zone == 55
     73        assert allclose(easting, 273741.297)
     74        assert allclose(northing, 5796489.777)
     75
     76       
     77    def test_UTM_3(self):
     78        #Test 3
     79        lat = degminsec2decimal_degrees(-60,0,0)
     80        lon = degminsec2decimal_degrees(130,0,0)
     81
     82        zone, easting, northing = redfearn(lat,lon)
     83        #print zone, easting, northing
     84
     85        assert zone == 52
     86        assert allclose(easting, 555776.267)
     87        assert allclose(northing, 3348167.264)
     88
     89
     90    def test_UTM_4(self):
     91
     92        #Test 4 (Kobenhavn, Northern hemisphere)
     93        lat = 55.70248
     94        dd,mm,ss = decimal_degrees2degminsec(lat)
     95
     96        lon = 12.58364
     97        dd,mm,ss = decimal_degrees2degminsec(lon)
     98
     99        zone, easting, northing = redfearn(lat,lon)
     100
     101
     102        assert zone == 33
     103        assert allclose(easting, 348157.631)
     104        assert allclose(northing, 6175612.993)
     105
     106
     107#-------------------------------------------------------------
     108if __name__ == "__main__":
     109
     110    mysuite = unittest.makeSuite(TestCase,'test')
     111    runner = unittest.TextTestRunner()
     112    runner.run(mysuite)
Note: See TracChangeset for help on using the changeset viewer.