1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
|
from pyproj import Geod
import commands, cPickle
g = Geod(ellps='clrk66')
lat1pt = 42.+(15./60.)
lon1pt = -71.-(7./60.)
lat2pt = 45.+(31./60.)
lon2pt = -123.-(41./60.)
print lat1pt,lon1pt,lat2pt,lon2pt
print 'inverse transform'
print 'from proj.4 invgeod:'
print commands.getoutput('echo "42d15\'N 71d07\'W 45d31\'N 123d41\'W" | geod +ellps=clrk66 -I -f "%.3f"')
print 'from pyproj._Geod._inv:'
az12,az21,dist = g.inv(lon1pt,lat1pt,lon2pt,lat2pt)
print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
print 'forward transform'
print 'from proj.4 geod:'
print commands.getoutput('echo "42d15\'N 71d07\'W -66d31\'50.141 4164192.708" | geod +ellps=clrk66 -f "%.3f"')
endlon,endlat,backaz = g.fwd(lon1pt,lat1pt,az12,dist)
print 'from pyproj._Geod._fwd:'
print "%6.3f %6.3f %13.3f" % (endlat,endlon,backaz)
print 'intermediate points:'
print 'from geod with +lat_1,+lon_1,+lat_2,+lon_2,+n_S:'
points = '+lon_1=%s +lat_1=%s +lon_2=%s +lat_2=%s' % (lon1pt,lat1pt,lon2pt,lat2pt,)
print points
print commands.getoutput('geod +ellps=clrk66 -f "%.3f" +n_S=5 '+points)
print 'from pyproj._Geod._npts:'
npts = 4
lonlats = g.npts(lon1pt,lat1pt,lon2pt,lat2pt,npts)
lonprev = lon1pt
latprev = lat1pt
print dist/(npts+1)
print '%6.3f %7.3f' % (lat1pt, lon1pt)
for lon, lat in lonlats:
az12,az21,dist = g.inv(lonprev,latprev,lon,lat)
print '%6.3f %7.3f %11.3f' % (lat, lon, dist)
latprev = lat; lonprev = lon
az12,az21,dist = g.inv(lonprev,latprev,lon2pt,lat2pt)
print '%6.3f %7.3f %11.3f' % (lat2pt, lon2pt, dist)
# should raise exception (antipodal point)
try:
print 'testing antipodal point, should raise exception ...'
az12,az21,dist = g.inv(0.,90.,0.,-90.,0.)
except ValueError:
print 'OK'
else:
print 'not OK, no exception raised!'
# should raise exception (equatorial arc)
try:
print 'testing equatorial arc, should raise execption ...'
az12,az21,dist = g.inv(180.,0.0,200.,0.)
endlon,endlat,backaz = g.fwd(180.,0.0,az12,0.)
except ValueError:
print 'OK'
else:
print 'not OK, no exception raised!'
# specify the lat/lons of some cities.
boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.)
portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.)
newyork_lat = 40.+(47./60.); newyork_lon = -73.-(58./60.)
london_lat = 51.+(32./60.); london_lon = -(5./60.)
g1 = Geod(ellps='clrk66')
cPickle.dump(g1,open('geod1.pickle','wb'),-1)
g2 = Geod(ellps='WGS84')
cPickle.dump(g2,open('geod2.pickle','wb'),-1)
az12,az21,dist = g1.inv(boston_lon,boston_lat,portland_lon,portland_lat)
print "distance between boston and portland, clrk66:"
print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
print "distance between boston and portland, WGS84:"
az12,az21,dist = g2.inv(boston_lon,boston_lat,portland_lon,portland_lat)
print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
print "testing pickling of Geod instance"
g3 = cPickle.load(open('geod1.pickle','rb'))
g4 = cPickle.load(open('geod2.pickle','rb'))
az12,az21,dist = g3.inv(boston_lon,boston_lat,portland_lon,portland_lat)
print "distance between boston and portland, clrk66 (from pickle):"
print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
az12,az21,dist = g4.inv(boston_lon,boston_lat,portland_lon,portland_lat)
print "distance between boston and portland, WGS84 (from pickle):"
print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
g3 = Geod('+ellps=clrk66') # proj4 style init string
print 'inverse transform'
print 'from proj.4 invgeod:'
print commands.getoutput('echo "42d15\'N 71d07\'W 45d31\'N 123d41\'W" | geod +ellps=clrk66 -I -f "%.3f"')
print 'from pyproj._Geod._inv:'
az12,az21,dist = g3.inv(lon1pt,lat1pt,lon2pt,lat2pt)
print "%7.3f %6.3f %12.3f" % (az12,az21,dist)
|
