_proj.pyx

#cimport c_numpy
#c_numpy.import_array()

import math

cdef double _dg2rad, _rad2dg

_dg2rad = math.radians(1.)
_rad2dg = math.degrees(1.)
_doublesize = sizeof(double)
__version__ = "1.9.6"

cdef extern from "math.h":
    cdef enum:
        HUGE_VAL
        FP_NAN

cdef extern from "geodesic.h":
  struct geod_geodesic:
        pass
  struct geod_geodesicline:
        pass
  void geod_init(geod_geodesic* g, double a, double f)
  void geod_direct(geod_geodesic* g,\
              double lat1, double lon1, double azi1, double s12,\
              double* plat2, double* plon2, double* pazi2)
  void geod_inverse(geod_geodesic* g,\
               double lat1, double lon1, double lat2, double lon2,\
               double* ps12, double* pazi1, double* pazi2)
  void geod_lineinit(geod_geodesicline* l,\
               geod_geodesic* g,\
               double lat1, double lon1, double azi1, unsigned caps)
  void geod_position(geod_geodesicline* l, double s12,\
               double* plat2, double* plon2, double* pazi2);
  cdef enum:
      GEODESIC_VERSION_MAJOR
      GEODESIC_VERSION_MINOR
      GEODESIC_VERSION_PATCH

# define part of the struct PJconsts from projects.h
ctypedef void (*c_func_type)()

ctypedef struct PJconsts:
    void *ctx
    c_func_type fwd
    c_func_type inv
    # ignore all other components of this struct, we don't need them

cdef extern from "proj_api.h":
    ctypedef struct projUV:
        double u
        double v
    ctypedef void *projPJ
    ctypedef void *projCtx
    projPJ pj_init_plus(char *)
    projPJ pj_init_plus_ctx(projCtx, char *)
    projUV pj_fwd(projUV, projPJ)
    projUV pj_inv(projUV, projPJ)
    projPJ pj_latlong_from_proj(projPJ)
    int pj_transform(projPJ src, projPJ dst, long point_count, int point_offset,
                     double *x, double *y, double *z)
    int pj_is_latlong(projPJ)
    char *pj_get_def( projPJ pj, int options)
    int pj_is_geocent(projPJ)
    char *pj_strerrno(int)
    void pj_ctx_free( projCtx )
    int pj_ctx_get_errno( projCtx )
    projCtx pj_ctx_alloc()
    projCtx pj_get_default_ctx()
    void pj_free(projPJ)
    # pj_dalloc: useful for deallocating stuff like i.e. pj_get_def char buffer
    void pj_dalloc(void *)
    void pj_set_searchpath ( int count, char **path )
    void pj_dalloc(void *)
    cdef enum:
        PJ_VERSION

cdef extern from "Python.h":
    int PyObject_AsWriteBuffer(object, void **rbuf, Py_ssize_t *len)


# version number strings for proj.4 and Geod
if PJ_VERSION > 499:
# proj.4 Version 4.10.0 and later: PJ_VERSION=MMMNNNPP later where MMM, NNN, PP
# are the major, minor, and patch numbers
    proj_version_str = "{0}.{1}.{2}".format(PJ_VERSION // 10**5 % 1000,
                           PJ_VERSION // 10**2 % 1000, PJ_VERSION % 100)
else:
#  before proj.4 version 4.10.0: PJ_VERSION=MNP where M, N, and P are the major,
#   minor, and patch numbers;
    proj_version_str = "{0}.{1}.{2}".format(PJ_VERSION // 100 % 10,
                                 PJ_VERSION // 10 % 10, PJ_VERSION % 10)

geodesic_version_str = "{0}.{1}.{2}".format(GEODESIC_VERSION_MAJOR,
                         GEODESIC_VERSION_MINOR, GEODESIC_VERSION_PATCH)


def set_datapath(datapath):
    # convert the datapath into bytes if it is not already
    if not isinstance(datapath, bytes):
        datapath = datapath.encode('utf-8')
    cdef const char *searchpath = datapath
    pj_set_searchpath(1, &searchpath)

# deprecated: used in _proj.Proj.to_latlong
def _createproj(projstring):
    return Proj(projstring)

cdef class Proj:
    cdef projPJ projpj
    cdef projCtx projctx
    cdef public object proj_version
    cdef char *pjinitstring
    cdef public object srs

    def __cinit__(self, projstring):
        # setup proj initialization string.
        cdef int err
        self.srs = projstring
        bytestr = _strencode(projstring)
        self.pjinitstring = bytestr
        # initialize projection
        self.projctx = pj_ctx_alloc()
        self.projpj = pj_init_plus_ctx(self.projctx, self.pjinitstring)
        err = pj_ctx_get_errno(self.projctx)
        if err != 0:
             raise RuntimeError(pj_strerrno(err))
        self.proj_version = PJ_VERSION/100.

    def __dealloc__(self):
        """destroy projection definition"""
        pj_free(self.projpj)
        pj_ctx_free(self.projctx)

    def to_latlong_def(self):
        """return the definition string of the geographic (lat/lon)
        coordinate version of the current projection"""
        # This is a little hacky way of getting a latlong proj object
        # Maybe instead of this function the __cinit__ function can take a
        # Proj object and a type (where type = "geographic") as the libproj
        # java wrapper
        cdef projPJ llpj
        cdef char *cstring_def
        cdef int err

        llpj = pj_latlong_from_proj(self.projpj) # create temp proj
        if llpj is not NULL:
            cstring_def = pj_get_def(llpj, 0) # get definition c string
            pj_free(llpj) # deallocate temp proj
            if cstring_def is not NULL:
                try:
                    pystring_def = <bytes>cstring_def # copy to python string
                    return pystring_def
                finally:
                    pj_dalloc(cstring_def) # deallocate c string

        raise RuntimeError("could not create latlong definition")


    # deprecated : using in transform raised a TypeError in release 1.9.5.1
    # reported in issue #53, resolved in #73.
    def to_latlong(self):
        """return a new Proj instance which is the geographic (lat/lon)
        coordinate version of the current projection"""
        cdef projPJ llpj
        llpj = pj_latlong_from_proj(self.projpj)
        initstring = pj_get_def(llpj, 0) # this leaks the c char buffer
        pj_free(llpj)
        return _createproj(initstring)

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

    def _fwd(self, object lons, object lats, radians=False, errcheck=False):
        """
 forward transformation - lons,lats to x,y (done in place).
 if radians=True, lons/lats are radians instead of degrees.
 if errcheck=True, an exception is raised if the forward transformation is invalid.
 if errcheck=False and the forward transformation is invalid, no exception is
 raised and 1.e30 is returned.
        """
        cdef projUV projxyout, projlonlatin
        cdef Py_ssize_t buflenx, bufleny, ndim, i
        cdef double u, v
        cdef double *lonsdata
        cdef double *latsdata
        cdef void *londata
        cdef void *latdata
        cdef int err
        # if buffer api is supported, get pointer to data buffers.
        if PyObject_AsWriteBuffer(lons, &londata, &buflenx) <> 0:
            raise RuntimeError
        if PyObject_AsWriteBuffer(lats, &latdata, &bufleny) <> 0:
            raise RuntimeError
        # process data in buffer
        if buflenx != bufleny:
            raise RuntimeError("Buffer lengths not the same")
        ndim = buflenx//_doublesize
        lonsdata = <double *>londata
        latsdata = <double *>latdata
        for i from 0 <= i < ndim:
            # if inputs are nan's, return big number.
            if lonsdata[i] != lonsdata[i] or latsdata[i] != latsdata[i]:
                lonsdata[i]=1.e30; latsdata[i]=1.e30
                if errcheck:
                    raise RuntimeError('projection undefined')
                continue
            if radians:
                projlonlatin.u = lonsdata[i]
                projlonlatin.v = latsdata[i]
            else:
                projlonlatin.u = _dg2rad*lonsdata[i]
                projlonlatin.v = _dg2rad*latsdata[i]
            projxyout = pj_fwd(projlonlatin,self.projpj)
            if errcheck:
                err = pj_ctx_get_errno(self.projctx)
                if err != 0:
                     raise RuntimeError(pj_strerrno(err))
            # since HUGE_VAL can be 'inf',
            # change it to a real (but very large) number.
            # also check for NaNs.
            if projxyout.u == HUGE_VAL or\
               projxyout.u != projxyout.u:
                if errcheck:
                    raise RuntimeError('projection undefined')
                lonsdata[i] = 1.e30
            else:
                lonsdata[i] = projxyout.u
            if projxyout.v == HUGE_VAL or\
               projxyout.u != projxyout.u:
                if errcheck:
                    raise RuntimeError('projection undefined')
                latsdata[i] = 1.e30
            else:
                latsdata[i] = projxyout.v

    def _inv(self, object x, object y, radians=False, errcheck=False):
        """
 inverse transformation - x,y to lons,lats (done in place).
 if radians=True, lons/lats are radians instead of degrees.
 if errcheck=True, an exception is raised if the inverse transformation is invalid.
 if errcheck=False and the inverse transformation is invalid, no exception is
 raised and 1.e30 is returned.
        """
        cdef projUV projxyin, projlonlatout
        cdef Py_ssize_t buflenx, bufleny, ndim, i
        cdef double u, v
        cdef void *xdata
        cdef void *ydata
        cdef double *xdatab
        cdef double *ydatab
        # if buffer api is supported, get pointer to data buffers.
        if PyObject_AsWriteBuffer(x, &xdata, &buflenx) <> 0:
            raise RuntimeError
        if PyObject_AsWriteBuffer(y, &ydata, &bufleny) <> 0:
            raise RuntimeError
        # process data in buffer
        # (for numpy/regular python arrays).
        if buflenx != bufleny:
            raise RuntimeError("Buffer lengths not the same")
        ndim = buflenx//_doublesize
        xdatab = <double *>xdata
        ydatab = <double *>ydata
        for i from 0 <= i < ndim:
            # if inputs are nan's, return big number.
            if xdatab[i] != xdatab[i] or ydatab[i] != ydatab[i]:
                xdatab[i]=1.e30; ydatab[i]=1.e30
                if errcheck:
                    raise RuntimeError('projection undefined')
                continue
            projxyin.u = xdatab[i]
            projxyin.v = ydatab[i]
            projpj2 = <PJconsts *> self.projpj
            if (projpj2.inv == NULL):
                raise RuntimeError('inverse projection undefined')
            projlonlatout = pj_inv(projxyin,self.projpj)
            if errcheck:
                err = pj_ctx_get_errno(self.projctx)
                if err != 0:
                     raise RuntimeError(pj_strerrno(err))
            # since HUGE_VAL can be 'inf',
            # change it to a real (but very large) number.
            # also check for NaNs.
            if projlonlatout.u == HUGE_VAL or \
               projlonlatout.u != projlonlatout.u:
                if errcheck:
                    raise RuntimeError('projection undefined')
                xdatab[i] = 1.e30
            elif radians:
                xdatab[i] = projlonlatout.u
            else:
                xdatab[i] = _rad2dg*projlonlatout.u
            if projlonlatout.v == HUGE_VAL or \
               projlonlatout.v != projlonlatout.v:
                if errcheck:
                    raise RuntimeError('projection undefined')
                ydatab[i] = 1.e30
            elif radians:
                ydatab[i] = projlonlatout.v
            else:
                ydatab[i] = _rad2dg*projlonlatout.v

#   def _fwdn(self, c_numpy.ndarray lonlat, radians=False, errcheck=False):
#       """
#forward transformation - lons,lats to x,y (done in place).
#Uses ndarray of shape ...,2.
#if radians=True, lons/lats are radians instead of degrees.
#if errcheck=True, an exception is raised if the forward
#                transformation is invalid.
#if errcheck=False and the forward transformation is
#                invalid, no exception is
#                raised and 1.e30 is returned.
#       """
#       cdef projUV projxyout, projlonlatin
#       cdef projUV *llptr
#       cdef int err
#       cdef Py_ssize_t npts, i
#       npts = c_numpy.PyArray_SIZE(lonlat)//2
#       llptr = <projUV *>lonlat.data
#       for i from 0 <= i < npts:
#           if radians:
#               projlonlatin = llptr[i]
#           else:
#               projlonlatin.u = _dg2rad*llptr[i].u
#               projlonlatin.v = _dg2rad*llptr[i].v
#           projxyout = pj_fwd(projlonlatin,self.projpj)

#           if errcheck:
#               err = pj_ctx_get_errno(self.projctx)
#               if err != 0:
#                    raise RuntimeError(pj_strerrno(err))
#           # since HUGE_VAL can be 'inf',
#           # change it to a real (but very large) number.
#           if projxyout.u == HUGE_VAL:
#               llptr[i].u = 1.e30
#           else:
#               llptr[i].u = projxyout.u
#           if projxyout.v == HUGE_VAL:
#               llptr[i].u = 1.e30
#           else:
#               llptr[i].v = projxyout.v

#   def _invn(self, c_numpy.ndarray xy, radians=False, errcheck=False):
#       """
#inverse transformation - x,y to lons,lats (done in place).
#Uses ndarray of shape ...,2.
#if radians=True, lons/lats are radians instead of degrees.
#if errcheck=True, an exception is raised if the inverse transformation is invalid.
#if errcheck=False and the inverse transformation is invalid, no exception is
#raised and 1.e30 is returned.
#       """
#       cdef projUV projxyin, projlonlatout
#       cdef projUV *llptr
#       cdef Py_ssize_t npts, i
#       npts = c_numpy.PyArray_SIZE(xy)//2
#       llptr = <projUV *>xy.data

#       for i from 0 <= i < npts:
#           projxyin = llptr[i]
#           projlonlatout = pj_inv(projxyin, self.projpj)
#           if errcheck:
#               err = pj_ctx_get_errno(self.projctx)
#               if err != 0:
#                    raise RuntimeError(pj_strerrno(err))
#           # since HUGE_VAL can be 'inf',
#           # change it to a real (but very large) number.
#           if projlonlatout.u == HUGE_VAL:
#               llptr[i].u = 1.e30
#           elif radians:
#               llptr[i].u = projlonlatout.u
#           else:
#               llptr[i].u = _rad2dg*projlonlatout.u
#           if projlonlatout.v == HUGE_VAL:
#               llptr[i].v = 1.e30
#           elif radians:
#               llptr[i].v = projlonlatout.v
#           else:
#               llptr[i].v = _rad2dg*projlonlatout.v

    def is_latlong(self):
        # returns True if projection in geographic (lon/lat) coordinates
        cdef int i
        i = pj_is_latlong(self.projpj)
        if i:
            return True
        else:
            return False

    def is_geocent(self):
        # returns True if projection in geocentric (x/y) coordinates
        cdef int i
        i = pj_is_geocent(self.projpj)
        if i:
            return True
        else:
            return False

    def definition_string(self):
        # Returns projection definition
        cdef bytes py_def_string
        cdef char* c_def_string
        c_def_string = pj_get_def(self.projpj, 0)
        py_def_string = c_def_string
        pj_dalloc(c_def_string)
        return py_def_string

    def __repr__(self):
        return "{modname}.{classname}({srs!r}, preserve_units=True)".format(
                   modname=self.__module__, classname=self.__class__.__name__,
                   srs=self.srs)

def _transform(Proj p1, Proj p2, inx, iny, inz, radians):
    # private function to call pj_transform
    cdef void *xdata
    cdef void *ydata
    cdef void *zdata
    cdef double *xx
    cdef double *yy
    cdef double *zz
    cdef Py_ssize_t buflenx, bufleny, buflenz, npts, i
    cdef int err
    if PyObject_AsWriteBuffer(inx, &xdata, &buflenx) <> 0:
        raise RuntimeError
    if PyObject_AsWriteBuffer(iny, &ydata, &bufleny) <> 0:
        raise RuntimeError
    if inz is not None:
        if PyObject_AsWriteBuffer(inz, &zdata, &buflenz) <> 0:
            raise RuntimeError
    else:
        buflenz = bufleny
    if not (buflenx == bufleny == buflenz):
        raise RuntimeError('x,y and z must be same size')
    xx = <double *>xdata
    yy = <double *>ydata
    if inz is not None:
        zz = <double *>zdata
    npts = buflenx//8
    if not radians and p1.is_latlong():
        for i from 0 <= i < npts:
            xx[i] = xx[i]*_dg2rad
            yy[i] = yy[i]*_dg2rad
    if inz is not None:
        err = pj_transform(p1.projpj, p2.projpj, npts, 0, xx, yy, zz)
    else:
        err = pj_transform(p1.projpj, p2.projpj, npts, 0, xx, yy, NULL)
    if err != 0:
        raise RuntimeError(pj_strerrno(err))
    if not radians and p2.is_latlong():
        for i from 0 <= i < npts:
            xx[i] = xx[i]*_rad2dg
            yy[i] = yy[i]*_rad2dg

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?

def _transform_sequence(Proj p1, Proj p2, Py_ssize_t stride, inseq, bint radians, bint switch):
    # private function to itransform function
    cdef:
        void *buffer
        double *coords
        double *x
        double *y
        double *z
        Py_ssize_t buflen, npts, i, j
        int err

    if stride < 2:
        raise RuntimeError("coordinates must contain at least 2 values")
    if PyObject_AsWriteBuffer(inseq, &buffer, &buflen) <> 0:
        raise RuntimeError("object does not provide the python buffer writeable interface")

    coords = <double*>buffer
    npts = buflen // (stride * sizeof(double))

    if not radians and p1.is_latlong():
        for i from 0 <= i < npts:
            j = stride*i
            coords[j] *= _dg2rad
            coords[j+1] *= _dg2rad

    if not switch:
        x = coords
        y = coords + 1
    else:
        x = coords + 1
        y = coords

    if stride == 2:
        z = NULL
    else:
        z = coords + 2

    err = pj_transform(p1.projpj, p2.projpj, npts, stride, x, y, z)

    if err != 0:
        raise RuntimeError(pj_strerrno(err))

    if not radians and p2.is_latlong():
        for i from 0 <= i < npts:
            j = stride*i
            coords[j] *= _dg2rad
            coords[j+1] *= _dg2rad


cdef class Geod:
    cdef geod_geodesic _geod_geodesic
    cdef public object initstring

    def __cinit__(self, a, f, sphere, b, es):
        geod_init(&self._geod_geodesic, <double> a, <double> f)
        self.a = a
        self.f = f
        if isinstance(a, float) and a.is_integer():
            # convert 'a' only for initstring
            a = int(a)
        if f == 0.0:
            f = 0
        self.initstring = '+a=%s +f=%s' % (a, f)
        self.sphere = sphere
        self.b = b
        self.es = es

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

    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 lat1,lon1,az1,s12,plon2,plat2,pazi2
        cdef double *lonsdata
        cdef double *latsdata
        cdef double *azdata
        cdef double *distdata
        cdef void *londata
        cdef void *latdata
        cdef void *azdat
        cdef void *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 not radians:
                lon1 = lonsdata[i]
                lat1 = latsdata[i]
                az1 = azdata[i]
                s12 = distdata[i]
            else:
                lon1 = _rad2dg*lonsdata[i]
                lat1 = _rad2dg*latsdata[i]
                az1 = _rad2dg*azdata[i]
                s12 = distdata[i]
            geod_direct(&self._geod_geodesic, lat1, lon1, az1, s12,\
                   &plat2, &plon2, &pazi2)
            # back azimuth needs to be flipped 180 degrees
            # to match what proj4 geod utility produces.
            if pazi2 > 0:
                pazi2 = pazi2-180.
            elif pazi2 <= 0:
                pazi2 = pazi2+180.
            if not radians:
                lonsdata[i] = plon2
                latsdata[i] = plat2
                azdata[i] = pazi2
            else:
                lonsdata[i] = _dg2rad*plon2
                latsdata[i] = _dg2rad*plat2
                azdata[i] = _dg2rad*pazi2

    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 double lat1,lon1,lat2,lon2,pazi1,pazi2,ps12
        cdef Py_ssize_t buflenlons, buflenlats, buflenaz, buflend, ndim, i
        cdef double *lonsdata
        cdef double *latsdata
        cdef double *azdata
        cdef double *distdata
        cdef void *londata
        cdef void *latdata
        cdef void *azdat
        cdef void *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
        for i from 0 <= i < ndim:
            if radians:
                lon1 = _rad2dg*lonsdata[i]
                lat1 = _rad2dg*latsdata[i]
                lon2 = _rad2dg*azdata[i]
                lat2 = _rad2dg*distdata[i]
            else:
                lon1 = lonsdata[i]
                lat1 = latsdata[i]
                lon2 = azdata[i]
                lat2 = distdata[i]
            geod_inverse(&self._geod_geodesic, lat1, lon1, lat2, lon2,
                    &ps12, &pazi1, &pazi2)
            # back azimuth needs to be flipped 180 degrees
            # to match what proj4 geod utility produces.
            if pazi2 > 0:
                pazi2 = pazi2-180.
            elif pazi2 <= 0:
                pazi2 = pazi2+180.
            if radians:
                lonsdata[i] = _dg2rad*pazi1
                latsdata[i] = _dg2rad*pazi2
            else:
                lonsdata[i] = pazi1
                latsdata[i] = pazi2
            azdata[i] = ps12

    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,ps12,pazi1,pazi2,s12,plon2,plat2
        cdef geod_geodesicline line
        if radians:
            lon1 = _rad2dg*lon1
            lat1 = _rad2dg*lat1
            lon2 = _rad2dg*lon2
            lat2 = _rad2dg*lat2
        # do inverse computation to set azimuths, distance.
        # in proj 4.9.3 and later the next two steps can be replace by a call
        # to geod_inverseline with del_s = line.s13/(npts+1)
        geod_inverse(&self._geod_geodesic, lat1, lon1,  lat2, lon2,
                &ps12, &pazi1, &pazi2)
        geod_lineinit(&line, &self._geod_geodesic, lat1, lon1, pazi1, 0u)
        # distance increment.
        del_s = ps12/(npts+1)
        # initialize output tuples.
        lats = ()
        lons = ()
        # loop over intermediate points, compute lat/lons.
        for i from 1 <= i < npts+1:
            s12 = i*del_s
            geod_position(&line, s12, &plat2, &plon2, &pazi2);
            if radians:
                lats = lats + (_dg2rad*plat2,)
                lons = lons + (_dg2rad*plon2,)
            else:
                lats = lats + (plat2,)
                lons = lons + (plon2,)
        return lons, lats

    def __repr__(self):
        return "{modname}.{classname}({init!r})".format(modname=self.__module__,
                                              classname=self.__class__.__name__,
                                              init=self.initstring)