proj4.lisp

   1 ;;; PHOROS -- Photogrammetric Road Survey
   2 ;;; Copyright (C) 2011 Bert Burgemeister
   3 ;;;
   4 ;;; This program is free software; you can redistribute it and/or modify
   5 ;;; it under the terms of the GNU General Public License as published by
   6 ;;; the Free Software Foundation; either version 2 of the License, or
   7 ;;; (at your option) any later version.
   8 ;;;
   9 ;;; This program is distributed in the hope that it will be useful,
  10 ;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 ;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 ;;; GNU General Public License for more details.
  13 ;;;
  14 ;;; You should have received a copy of the GNU General Public License along
  15 ;;; with this program; if not, write to the Free Software Foundation, Inc.,
  16 ;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  17
  18
  19 (defpackage proj
  20   (:use :common-lisp :cffi)
  21   (:export cs2cs degrees-to-radians radians-to-degrees version)
  22   (:documentation "Interface to the PROJ.4 cartographic projection library."))
  23
  24 (in-package :proj)
  25
  26 (load-foreign-library '(:default "libproj"))
  27
  28 (cffi:defcfun "pj_transform" :int
  29   (src :pointer)
  30   (dst :pointer)
  31   (point_count :long)
  32   (point_offset :int)
  33   (x :pointer)
  34   (y :pointer)
  35   (z :pointer))
  36
  37 (cffi:defcfun "pj_free" :void
  38   (arg0 :pointer))
  39
  40 (cffi:defcfun "pj_init_plus" :pointer
  41   (arg0 :string))
  42
  43 (cffi:defcfun "pj_get_release" :string)
  44
  45 (defun version ()
  46   "PROJ.4 version.  Last working one was \"Rel. 4.7.1, 23 September 2009\""
  47   (pj-get-release))
  48
  49 (defun degrees-to-radians (degrees)
  50   "Convert degrees into radians."
  51   (* degrees  (/ pi 180)))
  52
  53 (defun radians-to-degrees (radians)
  54   "Convert radians into degrees."
  55   (* radians  (/ 180 pi)))
  56
  57 (defun cs2cs (point &key
  58               (source-cs "+proj=latlong +datum=WGS84")
  59               (destination-cs "+proj=latlong +datum=WGS84"))
  60   "Transform point (a list of (x y z)) from source-cs to
  61 destination-cs.  Geographic coordinates are in radians."
  62   (with-foreign-objects ((x :double) (y :double) (z :double)
  63                          (proj-source :pointer) (proj-destination :pointer))
  64     (setf proj-source (pj-init-plus source-cs)
  65           proj-destination (pj-init-plus destination-cs)
  66           (mem-ref x :double) (coerce (first point) 'double-float)
  67           (mem-ref y :double) (coerce (second point) 'double-float)
  68           (mem-ref z :double) (coerce (third point) 'double-float))
  69     (pj-transform proj-source proj-destination 1 1 x y z)
  70     (pj-free proj-source) (pj-free proj-destination)
  71     (list (mem-ref x :double) (mem-ref y :double) (mem-ref z :double))))