proj4-sh.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 (defpackage proj
  19   (:use :common-lisp :cffi)
  20   (:export cs2cs degrees-to-radians radians-to-degrees version)
  21   (:documentation
  22    "Interface to the PROJ.4 cartographic projection library.
  23 This is an alternative to the CFFI interface defined in proj4.lisp which doesn't work well with certain versions of proj4, see below.
  24
  25 ;;; (proj:version)
  26 ;;; \"Rel. 4.7.1, 23 September 2009\"
  27 ;;; (proj:cs2cs (list (proj:degrees-to-radians 9.28838684)
  28 ;;;                   (proj:degrees-to-radians 53.19274138)
  29 ;;;                   68.969)
  30 ;;;             :source-cs \"+proj=longlat +datum=WGS84 +no_defs\"
  31 ;;;             :destination-cs \"+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs\")
  32 ;;; (519267.42716064927 5893750.589875963 68.969)
  33 ;;; echo 9.28838684 -53.19274138 68.969 | cs2cs -f %.9f +proj=longlat +datum=WGS84 +no_defs +to  +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
  34 ;;; 519267.427160649    5893750.589875964 68.969000000
  35
  36 ;;; (proj:version)
  37 ;;; \"Rel. 4.6.1, 21 August 2008\"
  38 ;;; (proj:cs2cs (list (proj:degrees-to-radians 9.28838684)
  39 ;;;                   (proj:degrees-to-radians 53.19274138)
  40 ;;;                   68.969)
  41 ;;;             :source-cs \" +proj=longlat +datum=WGS84 +no_defs\"
  42 ;;;             :destination-cs \"+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs\")
  43 ;;; (519267.46293972526 5893728.74020727 68.969)
  44 ;;; echo 9.28838684 -53.19274138 68.969 | cs2cs -f %.9f +proj=longlat +datum=WGS84 +no_defs +to  +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
  45 ;;; 519267.427160649    5893750.589875964 68.969000000"))
  46
  47 (in-package :proj)
  48
  49 (defun version ()
  50   "PROJ.4 version."
  51   (multiple-value-bind (standard-output error-output exit-status)
  52       (trivial-shell:shell-command "cs2cs")
  53     (declare (ignore standard-output) (ignore exit-status))
  54     (car (cl-utilities:split-sequence #\Newline error-output))))
  55
  56 (defun degrees-to-radians (degrees)
  57   "Convert degrees to radians."
  58   (* degrees  (/ pi 180)))
  59
  60 (defun radians-to-degrees (radians)
  61   "Convert radians to degrees."
  62   (* radians  (/ 180 pi)))
  63
  64 (defun cs2cs (point &key
  65               (source-cs "+proj=latlong +datum=WGS84")
  66               (destination-cs "+proj=latlong +datum=WGS84"))
  67   "Transform point (a list of (x y z)) from source-cs to
  68 destination-cs.  Geographic coordinates are in radians."
  69   (flet ((degrees-p (number) (<= (abs number) 360)) ;kludges
  70          (radians-p (number) (<= (abs number) (* 2 pi))))
  71     (let ((command (format
  72                     nil
  73                     "cs2cs -f %.9f ~A +no_defs +to ~A +no_defs"
  74                     source-cs destination-cs))
  75           (horizontal-coordinates (butlast point))
  76           (height (third point))
  77           (*read-default-float-format* 'double-float))
  78       (when (every #'radians-p horizontal-coordinates)
  79         (setf horizontal-coordinates
  80               (mapcar #'proj:radians-to-degrees horizontal-coordinates)))
  81       (multiple-value-bind (standard-output error-output exit-status)
  82           (trivial-shell:shell-command
  83            command
  84            :input (format nil "~{~S ~}~S"
  85                           horizontal-coordinates
  86                           height))
  87         (unless (zerop exit-status)
  88           (error "Attempt to call `~A´ returned ~D: ~A"
  89                  command exit-status error-output))
  90         (with-input-from-string (stream standard-output)
  91           (let ((result
  92                  (loop
  93                     for number = (read stream nil)
  94                     repeat 3
  95                     while number
  96                     do (assert (numberp number)
  97                                () "cs2cs didn't return numbers.")
  98                     collect number)))
  99             (if (degrees-p (first result))
 100                 (setf (first result)
 101                       (proj:degrees-to-radians (first result))))
 102             (if (degrees-p (second result))
 103                 (setf (second result)
 104                       (proj:degrees-to-radians (second result))))
 105             result))))))
 106