1\documentclass[10pt,a4paper]{article} 2\usepackage{tocbibind} 3\usepackage{hyperref} 4\hypersetup{colorlinks,citecolor=red,linkcolor=red,urlcolor=red} 5 6\newcommand{\octproj}{\texttt{OctPROJ}} 7\newcommand{\proj}{\texttt{PROJ}} 8\newcommand{\octave}{GNU Octave} 9 10\title{\octproj\\Calling \proj{} from \octave\footnote{This document is 11 distributed under the terms of the GNU Free Documentation License. 12 Please, see \url{http://www.gnu.org/licenses/}.}} 13\author{Jos\'e Luis Garc\'ia Pallero\footnote{ETSI en Topograf\'ia, Geodesia y 14 Cartograf\'ia, Universidad Polit\'ecnica de Madrid. 15 \texttt{jlg.pallero@upm.es}, \texttt{jgpallero@gmail.com}.}} 16\date{May 16, 2020 (version 2.0.1)\\ 17 May 9, 2020 (version 2.0.0)\\ 18 April 2, 2019 (version 1.1.6)\\ 19 June 16, 2015 (version 1.1.5)\\ 20 February 13, 2015 (version 1.1.4)\\ 21 June 20, 2013 (version 1.1.3)\\ 22 October 1, 2012 (version 1.1.1)\\ 23 April 13, 2012 (version 1.1.0)\\ 24 May 13, 2011 (version 1.0.2)\\ 25 November 29, 2010 (version 1.0.1)\\ 26 February 9, 2010 (version 1.0.0)} 27 28\begin{document} 29\maketitle 30% \tableofcontents 31 32\begin{abstract} 33This is a small introduction to using the \octproj{} package. In this text, you 34can overview the basic usage of the functions in 35\octave\footnote{\url{http://www.octave.org}.}. If you need a detailed 36description about the options, available projections and coordinate reference 37systems, please visit the \proj{} website\footnote{\url{https://proj.org/}.}. 38Starting with \octproj{} 2.0.0 the code was upgraded to \proj{} version 6.3.0, 39so older versions of the library could not work. 40\end{abstract} 41 42\tableofcontents 43 44\section{Overview} 45 46\octproj{} allows you to perform cartographic projections and coordinate 47reference systems transformation in \octave{} using the \proj{} library. You can 48take the power of \proj{} using the facilities that \octave{} provides, without 49know the internals of the \proj{} C API. You can use the conversion programs 50coming with \proj{} distribution, but for use them in \octave{} you must make 51system calls. With \octproj{}, \proj{} can be integrated in \octave{} scripts in 52a natural way. 53 54\section{Installation} 55 56As several \octave{} packages, \octproj{} installation consists in compiling the 57C++ kernel sources (see section \ref{op-kw}), link them against \octave{} 58library to generate \texttt{*.oct} functions and copy this \texttt{*.oct} 59executables and other \texttt{*.m} functions into a working directory. 60 61The automagic procedure can be easily done by running the command: 62 63\begin{verbatim} 64octave:1> pkg install octproj-x.x.x.tar.gz 65\end{verbatim} 66where \texttt{x.x.x} is the version number. 67 68After that, the functions and documentation are installed in your machine and 69you are ready for use the package. 70 71\section{Kernel wrapper} 72\label{op-kw} 73 74The main functions (the functions which make the actual computations) are 75written in separate files: \texttt{projwrap.h} and \texttt{projwrap.c}. 76 77The files contain three functions: 78\begin{itemize} 79\item \texttt{proj\_fwd}: forward computation of geodetic to projected 80 coordinates. 81\item \texttt{proj\_inv}: inverse computation of projected to geodetic 82 coordinates. 83\item \texttt{proj\_transform}: general transformations. It is possible to make 84 forward, inverse and other transformations using only one function (see 85 \proj{} documentation). 86\end{itemize} 87 88\subsection{Error handling} 89 90Error handling in the kernel wrapper is based on error codes from \proj. 91Functions in \texttt{projwrap} return the \proj{} error code and the \proj{} 92text string error message, which can be catched in order to work in this case. 93 94The functions can stop the execution in presence of errors depending on the 95nature of the error. 96\begin{itemize} 97\item If exist an error involving a general parameter of the projection, the 98 execution stops. 99\item If an error is found due to a wrong or out of domain input coordinate, the 100 execution don't stops, but the error code is activated and the output 101 coordinate corresponding to the error position have a special value. 102\end{itemize} 103 104\section{\octave{} functions} 105 106Two types of functions are programmed for \octave: \texttt{*.oct} functions and 107\texttt{*.m} functions. 108 109\subsection{\texttt{*.oct} functions} 110\label{op-of} 111 112These functions are linked with \texttt{projwrap}. You can use it, but it is not 113recommended because the input arguments are more strict (always column vectors) 114than \texttt{*.m} functions and don't check for errors. 115 116The functions are: 117\begin{itemize} 118\item \texttt{\_op\_geod2geoc}: transformation between geodetic and cartesian 119 tridimensional geocentric coordinates. 120\item \texttt{\_op\_geoc2geod}: transformation between cartesian tridimensional 121 geocentric and geodetic coordinates. 122\item \texttt{\_op\_fwd}: forward projection (calls the function 123 \texttt{proj\_trans\_generic} with \texttt{PJ\_FWD} option). 124\item \texttt{\_op\_inv}: inverse projection (calls the function 125 \texttt{proj\_trans\_generic} with \texttt{PJ\_INV} option). 126\item \texttt{\_op\_transform}: general transformation (calls 127 \texttt{proj\_trans\_generic} with \texttt{PJ\_FWD} option). 128\end{itemize} 129 130\subsection{\texttt{*.m} functions} 131 132These functions make the computations by calling the \texttt{*.oct} functions. 133You must call these functions because you can use various types of input 134(scalars, vectors or matrices) and checking of input arguments (data type, 135dimensions) is performed. 136 137The functions are the same as in section \ref{op-of} (without the \texttt{\_} at 138the beginning of the name): 139\begin{itemize} 140\item \texttt{op\_geod2geoc}: calls \texttt{\_op\_geod2geoc}. 141\item \texttt{op\_geoc2geod}: calls \texttt{\_op\_geoc2geod}. 142\item \texttt{op\_fwd}: calls \texttt{\_op\_fwd}. 143\item \texttt{op\_inv}: calls \texttt{\_op\_inv}. 144\item \texttt{op\_transform}: calls \texttt{\_op\_transform}. 145\end{itemize} 146 147\subsection{Error handling} 148 149\texttt{*.oct} and \texttt{*.m} functions can emit errors or warnings, some due 150to errors in input arguments and other due to errors in functions from 151\texttt{projwrap} kernel execution (see section \ref{op-kw}). 152 153Errors due to wrong input arguments (data types, dimensions, etc.) can be only 154given for \texttt{*.m} functions and this is the reason because the use of these 155functions are recommended. In this case the execution is aborted and nothing is 156stored in output arguments. 157 158Errors due to the execution of \texttt{projwrap} kernel can be emitted for both 159\texttt{*.oct} and \texttt{*.m} functions. If the error is due to an erroneous 160projection parameter the execution is aborted and nothing is stored in output 161arguments; but if the error is due to a wrong or out of domain input coordinate, 162a warning is emitted and the execution has a normal end. 163 164\section{Examples} 165 166\subsection{Geodetic to geocentric and vice versa} 167 168\begin{verbatim} 169lon=-6*pi/180;lat=43*pi/180;h=1000; 170[x,y,z]=op_geod2geoc(lon,lat,h,6378388,1/297) 171x = 4647300.723262573 172y = -488450.9885681378 173z = 4328259.364257743 174 175[lon,lat,h]=op_geoc2geod(x,y,z,6378388,1/297); 176lon*=180/pi,lat*=180/pi,h 177lon = -6 178lat = 42.99999999999999 179h = 1000 180\end{verbatim} 181 182\subsection{Forward and inverse projection} 183 184\begin{verbatim} 185lon=-6*pi/180;lat=43*pi/180; 186[x,y]=op_fwd(lon,lat,'+proj=utm +lon_0=3w +ellps=GRS80') 187x = 255466.9805506805 188y = 4765182.932683380 189 190[lon,lat]=op_inv(x,y,'+proj=utm +lon_0=3w +ellps=GRS80'); 191lon*=180/pi,lat*=180/pi 192lon = -5.999999999999999 193lat = 43 194\end{verbatim} 195 196\subsection{Forward and inverse transformation: \texttt{op\_transform}} 197 198This function takes into account some of the new \proj{} capabilities. The main 199of them are: 200\begin{itemize} 201\item Projection and coordinate reference systems can be defined not only via 202 the classical \texttt{+proj=} format, but also using 203 EPSG\footnote{\url{http://www.epsg-registry.org/}.} codes. 204\item Time coordinate for dynamical coordinate reference systems can be used. 205\end{itemize} 206 207\subsubsection{With altitude and time} 208 209\begin{verbatim} 210lon=-6;lat=43;h=1000;t=1; 211[x,y,h,t]=op_transform(lon,lat,h,t,... 212 '+proj=latlong +ellps=GRS80',... 213 '+proj=utm +lon_0=3w +ellps=GRS80') 214x = 255466.9805506804 215y = 4765182.932683380 216h = 1000 217t = 1 218 219[lon,lat,h,t]=op_transform(x,y,h,t,... 220 '+proj=utm +lon_0=3w +ellps=GRS80',... 221 '+proj=latlong +ellps=GRS80') 222lon = -6 223lat = 43 224h = 1000 225t = 1 226\end{verbatim} 227 228Note that in this case, where the \texttt{+proj=latlong} was used, the input 229geodetic coordinates are provided in degrees instead of radians. This is a 230feature of \proj{} when \texttt{+proj=latlong} is used in coordinate reference 231system transformation, but not in cartographic projection, i.e., when it is used 232in \texttt{op\_transform}, but not in \texttt{op\_fwd} nor \texttt{op\_inv}. 233 234\subsubsection{With altitude} 235 236\begin{verbatim} 237lon=-6;lat=43;h=1000; 238[x,y,h]=op_transform(lon,lat,h,... 239 '+proj=latlong +ellps=GRS80',... 240 '+proj=utm +lon_0=3w +ellps=GRS80') 241x = 255466.9805506804 242y = 4765182.932683380 243h = 1000 244 245[lon,lat,h]=op_transform(x,y,h,... 246 '+proj=utm +lon_0=3w +ellps=GRS80',... 247 '+proj=latlong +ellps=GRS80') 248lon = -6 249lat = 43 250h = 1000 251\end{verbatim} 252 253\subsubsection{Without altitude} 254 255\begin{verbatim} 256lon=-6;lat=43; 257[x,y]=op_transform(lon,lat,'+proj=latlong +ellps=GRS80',... 258 '+proj=utm +lon_0=3w +ellps=GRS80') 259x = 255466.9805506804 260y = 4765182.932683380 261 262[lon,lat]=op_transform(x,y,'+proj=utm +lon_0=3w +ellps=GRS80',... 263 '+proj=latlong +ellps=GRS80') 264lon = -6 265lat = 43 266\end{verbatim} 267 268In all the preceding examples the input coordinates can be scalars, vectors or 269matrices, and height and/or time can be provided as empty matrices. 270 271\subsubsection{Using EPSG codes} 272 273Supose a transformation from UTM zone 30 coordinates to UTM zone 31 in the 274ETRS89 coordinate reference system. Using the classical definition the operation 275is performed as 276 277\begin{verbatim} 278x1=[600000;700000];y1=[4800000;4900000];z1=100; 279[x2,y2,z2]=op_transform(x1,y1,z1,... 280 '+proj=utm +lon_0=3w +ellps=GRS80',... 281 '+proj=utm +lon_0=3e +ellps=GRS80') 282x2 = 283 284 113677.9843623374 285 220776.6200746754 286 287y2 = 288 289 4810303.384473779 290 4902896.002979981 291 292z2 = 293 294 100 295 100 296\end{verbatim} 297 298The same transformation can be carried out using the corresponding EPSG 299codes, i.e., \texttt{EPSG:25830} and 300\texttt{EPSG:25831}\footnote{See 301\url{https://spatialreference.org/ref/epsg/25830/}, 302and \url{https://spatialreference.org/ref/epsg/25831/}.}: 303 304\begin{verbatim} 305x1=[600000;700000];y1=[4800000;4900000];z1=100; 306[x2,y2,z2]=op_transform(x1,y1,z1,'EPSG:25830','EPSG:25831') 307x2 = 308 309 113677.9843623374 310 220776.6200746754 311 312y2 = 313 314 4810303.384473779 315 4902896.002979981 316 317z2 = 318 319 100 320 100 321\end{verbatim} 322 323\subsection{Error due to an erroneous parameter} 324 325\begin{verbatim} 326lon=-6*pi/180;lat=43*pi/180; 327[x,y]=op_fwd(lon,lat,'+proj=utm +lon_0=3w +ellps=GRS8') 328proj_create: Error -9: unknown elliptical parameter name 329error: 330 In function op_fwd: 331 In function _op_fwd: 332 Error in projection parameters 333 unknown elliptical parameter name 334 +proj=utm +lon_0=3w +ellps=GRS8 335error: called from 336 op_fwd at line 96 column 5 337\end{verbatim} 338 339\subsection{Error due to latitude too big} 340 341\begin{verbatim} 342lon=[-6*pi/180;-6*pi/180];lat=[43*pi/180;43]; 343[x,y]=_op_fwd(lon,lat,'+proj=utm +lon_0=3w +ellps=GRS80') 344warning: _op_fwd: 345 346warning: Projection error in point 2 (index starts at 1) 347x = 348 349 255466.98055 350 Inf 351 352y = 353 354 4765182.93268 355 Inf 356\end{verbatim} 357 358\section{Notes} 359 360Apart from \url{http://octave.sourceforge.net/octproj/index.html}, an up to date 361version of \octproj{} can be downloaded from 362\url{https://bitbucket.org/jgpallero/octproj/}. 363 364\begin{thebibliography}{99} 365\bibitem{eat-om} \textsc{Eaton}, John W.; \textsc{Bateman}, David; 366 \textsc{Hauberg}, S\o{}ren; and \textsc{Wehbring}, Rik; 367 GNU Octave. A high-level interactive language for numerical 368 computations; Edition 5 for Octave version 5.2.0, January 2020; 369 \url{https://www.gnu.org/software/octave/octave.pdf}; 370 Permanently updated at 371 \url{https://www.gnu.org/software/octave/support.html}. 372\bibitem{projman} \textsc{Evenden}, Gerald I.; Cartographic Projection 373 Procedures for the UNIX Environment---A User's Manual; USGS 374 Open-File Report 90-284; 2003; 375 \url{ftp://ftp.remotesensing.org/proj/OF90-284.pdf}. 376\bibitem{projir1} \textsc{Evenden}, Gerald I.; Cartographic Projection 377 Procedures Release 4, Interim Report; 2003; 378 \url{ftp://ftp.remotesensing.org/proj/proj.4.3.pdf}. 379\bibitem{projir2} \textsc{Evenden}, Gerald I.; Cartographic Projection 380 Procedures Release 4, Second Interim Report; 2003; 381 \url{ftp://ftp.remotesensing.org/proj/proj.4.3.I2.pdf}. 382\bibitem{sny-wm} \textsc{Snyder}, John Parr; Map Projections: A Working Manual; 383 USGS series, Professional Paper 1395; Geological Survey 384 (U. S.), 1987; 385 \url{http://pubs.er.usgs.gov/usgspubs/pp/pp1395}. 386\end{thebibliography} 387 388\end{document} 389 390%Copyright (C) 2009-2020, José Luis García Pallero, <jgpallero@gmail.com> 391%This document is distributed under the terms of the GNU Free Documentation 392%License. Please, see http://www.gnu.org/licenses/ 393