Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

Close

[d89be7]: src / _op_inv.cc Maximize Restore History

Download this file

_op_inv.cc    161 lines (160 with data), 6.1 kB

  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
/* -*- coding: utf-8 -*- */
/* Copyright (C) 2009 José Luis García Pallero, <jgpallero@gmail.com>
*
* This file is part of OctPROJ.
*
* OctPROJ is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this software; see the file COPYING. If not, see
* <http://www.gnu.org/licenses/>.
*/
/******************************************************************************/
/******************************************************************************/
#define HELPTEXT "\
-*- texinfo -*-\n\
@deftypefn{Loadable Function}{[@var{lon},@var{lat}] =}\
_op_inv(@var{X},@var{Y},@var{params})\n\
\n\
@cindex Performs inverse projection step.\n\
\n\
This function unprojects cartesian projected coordinates (in a defined\n\
cartographic projection) into geodetic coordinates using the PROJ.4 function \
pj_inv().\n\
\n\
@var{X} is a column vector containing the X projected coordinates.\n\
@var{Y} is a column vector containing the Y projected coordinates.\n\
@var{params} is a text string containing the projection parameters in PROJ.4 \
format.\n\
\n\
The coordinate vectors @var{X} and @var{Y} must be both scalars or both\n\
column vectors (of the same size).\n\
\n\
@var{lon} is a column vector containing the geodetic longitude, in radians.\n\
@var{lat} is a column vector containing the geodetic latitude, in radians.\n\
\n\
If a projection error occurs, the resultant coordinates for the affected\n\
points have both Inf value and a warning message is emitted (one for each\n\
erroneous point).\n\
@seealso{_op_fwd}\n\
@end deftypefn"
/******************************************************************************/
/******************************************************************************/
#include<octave/oct.h>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include"projwrap.h"
/******************************************************************************/
/******************************************************************************/
#define ERRORTEXT 1000
/******************************************************************************/
/******************************************************************************/
DEFUN_DLD(_op_inv,args,,HELPTEXT)
{
//error message
char errorText[ERRORTEXT]="_op_inv:\n\t";
//output list
octave_value_list outputList;
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
//checking input arguments
if(args.length()!=3)
{
//error text
sprintf(&errorText[strlen(errorText)],
"Incorrect number of input arguments\n\t"
"See help _op_inv");
//error message
error(errorText);
}
else
{
//error code
int idErr=0;
//error in projection
int projectionError=0;
//auxiliar message
char auxErrorText[ERRORTEXT];
//loop index
size_t i=0;
//projected coordinates
ColumnVector X=args(0).column_vector_value();
ColumnVector Y=args(1).column_vector_value();
//projection parameters
std::string params=args(2).string_value();
//number of elements
size_t nElem=static_cast<size_t>(X.rows());
//geodetic coordinates
ColumnVector lonOut(nElem);
ColumnVector latOut(nElem);
//computation vectors
double* lon=NULL;
double* lat=NULL;
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
//copy input data
for(i=0;i<nElem;i++)
{
lonOut(i) = X(i);
latOut(i) = Y(i);
}
//pointers to output data
lon = const_cast<double*>(lonOut.data());
lat = const_cast<double*>(latOut.data());
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
//projection
idErr = proj_inv(lon,lat,nElem,params.c_str(),&errorText[strlen(errorText)],
&projectionError);
//error checking
if(idErr)
{
//type of error
if(projectionError)
{
//warning message
warning(errorText);
//positions of projection errors
for(i=0;i<nElem;i++)
{
//testing for HUGE_VAL
if((lon[i]==HUGE_VAL)&&(lat[i]==HUGE_VAL))
{
//error text
sprintf(auxErrorText,
"Projection error in point %d (index starts at 1)",
i+1);
//warning
warning(auxErrorText);
}
}
}
else
{
//error message
error(errorText);
//exit
return outputList;
}
}
////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////
//output parameters list
outputList(0) = lonOut;
outputList(1) = latOut;
}
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
//exit
return outputList;
}