File: _geod.pyx

package info (click to toggle)
python-pyproj 1.8.9-1
  • links: PTS, VCS
  • area: main
  • in suites: wheezy
  • size: 9,192 kB
  • sloc: ansic: 10,983; python: 510; makefile: 6
file content (175 lines) | stat: -rw-r--r-- 7,667 bytes parent folder | download | duplicates (2)
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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
# Make changes to this file, not the c-wrappers that Pyrex generates.

include "_pyproj.pxi"

cdef class Geod:
    cdef GEODESIC_T geodesic_t
    cdef public object geodstring
    cdef public object proj_version
    cdef char *geodinitstring

    def __cinit__(self, geodstring):
        cdef GEODESIC_T GEOD_T
        # setup geod initialization string.
        self.geodstring = geodstring
        bytestr = _strencode(geodstring)
        self.geodinitstring = bytestr
        # initialize projection
        self.geodesic_t = GEOD_init_plus(self.geodinitstring, &GEOD_T)[0]
        if pj_errno != 0:
            raise RuntimeError(pj_strerrno(pj_errno))
        self.proj_version = PJ_VERSION/100.

    def __reduce__(self):
        """special method that allows pyproj.Geod instance to be pickled"""
        return (self.__class__,(self.geodstring,))

    def _fwd(self, object lons, object lats, object az, object dist, radians=False):
        """
 forward transformation - determine longitude, latitude and back azimuth 
 of a terminus point given an initial point longitude and latitude, plus
 forward azimuth and distance.
 if radians=True, lons/lats are radians instead of degrees.
        """
        cdef Py_ssize_t buflenlons, buflenlats, buflenaz, buflend, ndim, i
        cdef double *lonsdata, *latsdata, *azdata, *distdata
        cdef void *londata, *latdata, *azdat, *distdat
        # if buffer api is supported, get pointer to data buffers.
        if PyObject_AsWriteBuffer(lons, &londata, &buflenlons) <> 0:
            raise RuntimeError
        if PyObject_AsWriteBuffer(lats, &latdata, &buflenlats) <> 0:
            raise RuntimeError
        if PyObject_AsWriteBuffer(az, &azdat, &buflenaz) <> 0:
            raise RuntimeError
        if PyObject_AsWriteBuffer(dist, &distdat, &buflend) <> 0:
            raise RuntimeError
        # process data in buffer
        if not buflenlons == buflenlats == buflenaz == buflend:
            raise RuntimeError("Buffer lengths not the same")
        ndim = buflenlons//_doublesize
        lonsdata = <double *>londata
        latsdata = <double *>latdata
        azdata = <double *>azdat
        distdata = <double *>distdat
        for i from 0 <= i < ndim:
            if radians:
                self.geodesic_t.p1.v = lonsdata[i]
                self.geodesic_t.p1.u = latsdata[i]
                self.geodesic_t.ALPHA12 = azdata[i]
                self.geodesic_t.DIST = distdata[i]
            else:
                self.geodesic_t.p1.v = _dg2rad*lonsdata[i]
                self.geodesic_t.p1.u = _dg2rad*latsdata[i]
                self.geodesic_t.ALPHA12 = _dg2rad*azdata[i]
                self.geodesic_t.DIST = distdata[i]
            geod_pre(&self.geodesic_t)
            if pj_errno != 0:
                raise RuntimeError(pj_strerrno(pj_errno))
            geod_for(&self.geodesic_t)
            if pj_errno != 0:
                raise RuntimeError(pj_strerrno(pj_errno))
            # check for NaN.
            if self.geodesic_t.ALPHA21 != self.geodesic_t.ALPHA21:
                raise ValueError('undefined inverse geodesic (may be an antipodal point)')
            if radians:
                lonsdata[i] = self.geodesic_t.p2.v
                latsdata[i] = self.geodesic_t.p2.u
                azdata[i] = self.geodesic_t.ALPHA21
            else:
                lonsdata[i] = _rad2dg*self.geodesic_t.p2.v
                latsdata[i] = _rad2dg*self.geodesic_t.p2.u
                azdata[i] = _rad2dg*self.geodesic_t.ALPHA21

    def _inv(self, object lons1, object lats1, object lons2, object lats2, radians=False):
        """
 inverse transformation - return forward and back azimuths, plus distance
 between an initial and terminus lat/lon pair.
 if radians=True, lons/lats are radians instead of degrees.
        """
        cdef Py_ssize_t buflenlons, buflenlats, buflenaz, buflend, ndim, i
        cdef double *lonsdata, *latsdata, *azdata, *distdata
        cdef void *londata, *latdata, *azdat, *distdat
        # if buffer api is supported, get pointer to data buffers.
        if PyObject_AsWriteBuffer(lons1, &londata, &buflenlons) <> 0:
            raise RuntimeError
        if PyObject_AsWriteBuffer(lats1, &latdata, &buflenlats) <> 0:
            raise RuntimeError
        if PyObject_AsWriteBuffer(lons2, &azdat, &buflenaz) <> 0:
            raise RuntimeError
        if PyObject_AsWriteBuffer(lats2, &distdat, &buflend) <> 0:
            raise RuntimeError
        # process data in buffer
        if not buflenlons == buflenlats == buflenaz == buflend:
            raise RuntimeError("Buffer lengths not the same")
        ndim = buflenlons//_doublesize
        lonsdata = <double *>londata
        latsdata = <double *>latdata
        azdata = <double *>azdat
        distdata = <double *>distdat
        errmsg = 'undefined inverse geodesic (may be an antipodal point)'
        for i from 0 <= i < ndim:
            if radians:
                self.geodesic_t.p1.v = lonsdata[i]
                self.geodesic_t.p1.u = latsdata[i]
                self.geodesic_t.p2.v = azdata[i]
                self.geodesic_t.p2.u = distdata[i]
            else:
                self.geodesic_t.p1.v = _dg2rad*lonsdata[i]
                self.geodesic_t.p1.u = _dg2rad*latsdata[i]
                self.geodesic_t.p2.v = _dg2rad*azdata[i]
                self.geodesic_t.p2.u = _dg2rad*distdata[i]
            geod_inv(&self.geodesic_t)
            if self.geodesic_t.DIST != self.geodesic_t.DIST: # check for NaN
                raise ValueError('undefined inverse geodesic (may be an antipodal point)')
            if pj_errno != 0:
                raise RuntimeError(pj_strerrno(pj_errno))
            if radians:
                lonsdata[i] = self.geodesic_t.ALPHA12
                latsdata[i] = self.geodesic_t.ALPHA21
            else:
                lonsdata[i] = _rad2dg*self.geodesic_t.ALPHA12
                latsdata[i] = _rad2dg*self.geodesic_t.ALPHA21
            azdata[i] = self.geodesic_t.DIST

    def _npts(self, double lon1, double lat1, double lon2, double lat2, int npts, radians=False):
        """
 given initial and terminus lat/lon, find npts intermediate points."""
        cdef int i
        cdef double del_s
        if radians:
            self.geodesic_t.p1.v = lon1
            self.geodesic_t.p1.u = lat1
            self.geodesic_t.p2.v = lon2
            self.geodesic_t.p2.u = lat2
        else:
            self.geodesic_t.p1.v = _dg2rad*lon1
            self.geodesic_t.p1.u = _dg2rad*lat1
            self.geodesic_t.p2.v = _dg2rad*lon2
            self.geodesic_t.p2.u = _dg2rad*lat2
        # do inverse computation to set azimuths, distance.
        geod_inv(&self.geodesic_t)
        # set up some constants needed for forward computation.
        geod_pre(&self.geodesic_t)
        # distance increment.
        del_s = self.geodesic_t.DIST/(npts+1)
        # initialize output tuples.
        lats = ()
        lons = ()
        # loop over intermediate points, compute lat/lons.
        for i from 1 <= i < npts+1:
            self.geodesic_t.DIST = i*del_s
            geod_for(&self.geodesic_t)
            if radians:
                lats = lats + (self.geodesic_t.p2.u,)
                lons = lons + (self.geodesic_t.p2.v,)
            else:
                lats = lats + (_rad2dg*self.geodesic_t.p2.u,)
                lons = lons + (_rad2dg*self.geodesic_t.p2.v,)
        return lons, lats   

cdef _strencode(pystr,encoding='ascii'):
    # encode a string into bytes.  If already bytes, do nothing.
    try:
        return pystr.encode(encoding)
    except AttributeError:
        return pystr # already bytes?