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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
|
/***********************************************************************
* GRASS 5.0 gis library
* get_datum.c, get datum parameters from location database
*
* Andreas Lange, andreas.lange@rhein-main.de
* version 0.9
* modified Jul 13 2000
*
***********************************************************************/
#include "gis.h"
#include "glocale.h"
#include <unistd.h>
#include <ctype.h>
#include <string.h>
#include <stdlib.h>
static int getvalue(const char *, double *);
static char *PERMANENT = "PERMANENT";
/* G_get_datum_parameters
*
* This routine returns the datum parameters from the database.
* If the PROJECTION_FILE exists in the PERMANENT mapset, read info from
* that file, otherwise return WGS 84 values.
*
* Returns: 1 ok, 0 default values used, -1 on internal error.
*/
/*!
* \brief get datum parameters from database
*
* This function sets the datum parameters for the map datum of the current database.
* These are the semi-major axis <b>a</b> (in meters), the eccentricity squared
* <b>e2</b> and the inverse flattening <b>f</b> of the spheroid associated
* with the database and the x shift <b>dx</b>, the y shift <b>dy</b> and the
* z shift <b>dz</b> of the map datum associated with the database. If there is
* no map datum explicitely associated with the actual database, the standard
* values for the WGS84 spheroid and map datum are set. The funcion returns 1 on
* success, 0 if the default WGS84 parameters are set. If an error occurs, the
* function dies with a diagnostic message (HINT: to change, very bad practice to
* die in a library function!).
*
* \param a
* \param e2
* \param f
* \param dx
* \param dy
* \param dz
* \return int
*/
int
G_get_datum_parameters (double *a, double *e2, double *f,
double *dx, double *dy, double *dz)
{
int in_stat, get_parms = 0;
char *str, *dat, *ellps, err[1024], ipath[1024];
struct Key_Value *proj_keys;
G__file_name (ipath, "", PROJECTION_FILE, PERMANENT);
/* no PROJ_INFO file existing, return default wgs84 values */
if (access(ipath,0) !=0)
{
*a = 6378137.0 ;
*e2 = .006694385 ;
*f = 298.257223563;
*dx = 0.0;
*dy = 0.0;
*dz = 0.0;
return 0;
}
/* read in key values */
proj_keys = G_read_key_value_file(ipath, &in_stat);
if (in_stat !=0)
{
sprintf (err, _("Unable to open file %s in %s"),PROJECTION_FILE,PERMANENT);
G_fatal_error (err);
}
/* if datum key is existing, process values */
if ((dat = G_find_key_value("datum",proj_keys))!=NULL) {
str = G_find_key_value("a",proj_keys);
get_parms += getvalue(str, a);
str = G_find_key_value("es",proj_keys);
get_parms += getvalue(str, e2);
str = G_find_key_value("f",proj_keys);
get_parms += getvalue(str, f);
str = G_find_key_value("dx",proj_keys);
get_parms += getvalue(str, dx);
str = G_find_key_value("dy",proj_keys);
get_parms += getvalue(str, dy);
str = G_find_key_value("dz",proj_keys);
get_parms += getvalue(str, dz);
if (get_parms) {
if (G_datum_shift(G_get_datum_by_name(dat), dx, dy, dz) == 0) {
sprintf(err, _("Error reading datum shift parameters for %s from table"), dat);
G_fatal_error(err);
return -1;
}
/* get ellipsoid parameters from PROJ_INFO */
ellps = G_find_key_value("ellps",proj_keys);
if (ellps!=NULL) {
if (G_get_spheroid_by_name(ellps, a, e2, f) == 0) {
sprintf(err, _("Error reading ellipsoid parameters for %s from table"), ellps);
G_fatal_error(err);
return -1;
}
} else {
sprintf(err, _("No ellipsoid field %s in file %s in %s"), ellps, PROJECTION_FILE,PERMANENT);
G_fatal_error (err);
return -1;
}
}
return 1;
} else {
/* no datum key in PROJ_INFO, supply wgs84 params */
*a = 6378137.0 ;
*e2 = .006694385 ;
*f = 298.257223563;
*dx = 0.0;
*dy = 0.0;
*dz = 0.0;
return 0;
}
}
/* placeholder, not yet implemented */
/*!
* \brief get datum parameters from database
*
* This is a placeholder as the 7
* parameter datum shift support is not implemented yet.
*
* \param a
* \param e2
* \param f
* \param dx
* \param dy
* \param dz
* \param rx
* \param ry
* \param rz
* \param m
* \return int
*/
int
G_get_datum_parameters7(double *a, double *e2, double *f,
double *dx, double *dy, double *dz,
double *rx, double *ry, double *rz, double *m)
{
return -1;
}
static int
getvalue (const char *key, double *value)
{
char err[512];
if (key!=NULL) {
if(sscanf(key,"%lf",value)!=1) {
sprintf (err, _("invalid value: field %s in file %s in %s")
,key,PROJECTION_FILE,PERMANENT);
G_fatal_error (err);
} else {
return 0;
}
}
return 1;
}
|
