1
/***************************************************************************/
2
/* RSC IDENTIFIER: ALBERS
6
* This component provides conversions between Geodetic coordinates
7
* (latitude and longitude in radians) and Albers Equal Area Conic
8
* projection coordinates (easting and northing in meters) defined
9
* by two standard parallels.
13
* This component checks parameters for valid values. If an invalid value
14
* is found the error code is combined with the current error code using
15
* the bitwise or. This combining allows multiple error codes to be
16
* returned. The possible error codes are:
18
* ALBERS_NO_ERROR : No errors occurred in function
19
* ALBERS_LAT_ERROR : Latitude outside of valid range
21
* ALBERS_LON_ERROR : Longitude outside of valid range
22
* (-180 to 360 degrees)
23
* ALBERS_EASTING_ERROR : Easting outside of valid range
24
* (depends on ellipsoid and projection
26
* ALBERS_NORTHING_ERROR : Northing outside of valid range
27
* (depends on ellipsoid and projection
29
* ALBERS_FIRST_STDP_ERROR : First standard parallel outside of valid
30
* range (-90 to 90 degrees)
31
* ALBERS_SECOND_STDP_ERROR : Second standard parallel outside of valid
32
* range (-90 to 90 degrees)
33
* ALBERS_ORIGIN_LAT_ERROR : Origin latitude outside of valid range
35
* ALBERS_CENT_MER_ERROR : Central meridian outside of valid range
36
* (-180 to 360 degrees)
37
* ALBERS_A_ERROR : Semi-major axis less than or equal to zero
38
* ALBERS_INV_F_ERROR : Inverse flattening outside of valid range
40
* ALBERS_HEMISPHERE_ERROR : Standard parallels cannot be opposite
42
* ALBERS_FIRST_SECOND_ERROR : The 1st & 2nd standard parallels cannot
48
* ALBERS is intended for reuse by any application that performs an Albers
49
* Equal Area Conic projection or its inverse.
53
* Further information on ALBERS can be found in the Reuse Manual.
55
* ALBERS originated from: U.S. Army Topographic Engineering Center
56
* Geospatial Information Division
58
* Alexandria, VA 22310-3864
62
* None apply to this component.
66
* ALBERS has no restrictions.
70
* ALBERS was tested and certified in the following environments:
72
* 1. Solaris 2.5 with GCC, version 2.8.1
73
* 2. MSDOS with MS Visual C++, version 6
79
* 07-09-99 Original Code
86
/***************************************************************************/
95
* math.h - Standard C math library
96
* albers.h - Is for prototype error checking
100
/***************************************************************************/
105
#define PI 3.14159265358979323e0 /* PI */
106
#define PI_OVER_2 ( PI / 2.0)
107
#define TWO_PI (2.0 * PI)
108
#define ES_SIN(sinlat) (es * sinlat)
109
#define ONE_MINUS_SQR(x) (1.0 - x * x)
110
#define ALBERS_M(clat,oneminussqressin) (clat / sqrt(oneminussqressin))
111
#define ALBERS_Q(slat,oneminussqressin,essin) (one_MINUS_es2)*(slat / (oneminussqressin)- \
112
(1 / (two_es)) *log((1 - essin) / (1 + essin)))
114
/***************************************************************************/
119
/* Ellipsoid Parameters, default to WGS 84 */
120
static double Albers_a = 6378137.0; /* Semi-major axis of ellipsoid in meters */
121
static double Albers_f = 1 / 298.257223563; /* Flattening of ellipsoid */
122
static double es = 0.08181919084262188000; /* Eccentricity of ellipsoid */
123
static double es2 = 0.0066943799901413800; /* Eccentricity squared */
124
static double C =1.4896626908850; /* constant c */
125
static double rho0 = 6388749.3391064; /* height above ellipsoid */
126
static double n = .70443998701755; /* ratio between meridians */
127
static double Albers_a_OVER_n = 9054194.9882824; /* Albers_a / n */
128
static double one_MINUS_es2 = .99330562000986; /* 1 - es2 */
129
static double two_es = .16363838168524; /* 2 * es */
131
/* Albers Projection Parameters */
132
static double Albers_Origin_Lat = (45 * PI / 180); /* Latitude of origin in radians */
133
static double Albers_Origin_Long = 0.0; /* Longitude of origin in radians */
134
static double Albers_Std_Parallel_1 = (40 * PI / 180);
135
static double Albers_Std_Parallel_2 = (50 * PI / 180);
136
static double Albers_False_Easting = 0.0;
137
static double Albers_False_Northing = 0.0;
139
static double Albers_Delta_Northing = 40000000;
140
static double Albers_Delta_Easting = 40000000;
142
* These state variables are for optimization purposes. The only function
143
* that should modify them is Set_Albers_Parameters.
147
/***************************************************************************/
151
long Set_Albers_Parameters(double a,
153
double Origin_Latitude,
154
double Central_Meridian,
155
double Std_Parallel_1,
156
double Std_Parallel_2,
157
double False_Easting,
158
double False_Northing)
160
{ /* BEGIN Set_Albers_Parameters */
162
* The function Set_Albers_Parameters receives the ellipsoid parameters and
163
* projection parameters as inputs, and sets the corresponding state
164
* variables. If any errors occur, the error code(s) are returned by the function,
165
* otherwise ALBERS_NO_ERROR is returned.
167
* a : Semi-major axis of ellipsoid, in meters (input)
168
* f : Flattening of ellipsoid (input)
169
* Origin_Latitude : Latitude in radians at which the (input)
170
* point scale factor is 1.0
171
* Central_Meridian : Longitude in radians at the center of (input)
173
* Std_Parallel_1 : First standard parallel (input)
174
* Std_Parallel_2 : Second standard parallel (input)
175
* False_Easting : A coordinate value in meters assigned to the
176
* central meridian of the projection. (input)
177
* False_Northing : A coordinate value in meters assigned to the
178
* origin latitude of the projection (input)
181
double sin_lat, sin_lat_1, cos_lat;
182
double m1, m2, SQRm1;
184
double es_sin, one_MINUS_SQRes_sin;
186
double inv_f = 1 / f;
187
long Error_Code = ALBERS_NO_ERROR;
190
{ /* Semi-major axis must be greater than zero */
191
Error_Code |= ALBERS_A_ERROR;
193
if ((inv_f < 250) || (inv_f > 350))
194
{ /* Inverse flattening must be between 250 and 350 */
195
Error_Code |= ALBERS_INV_F_ERROR;
197
if ((Origin_Latitude < -PI_OVER_2) || (Origin_Latitude > PI_OVER_2))
198
{ /* origin latitude out of range */
199
Error_Code |= ALBERS_ORIGIN_LAT_ERROR;
201
if ((Central_Meridian < -PI) || (Central_Meridian > TWO_PI))
202
{ /* origin longitude out of range */
203
Error_Code |= ALBERS_CENT_MER_ERROR;
205
if ((Std_Parallel_1 < -PI_OVER_2) || (Std_Parallel_1 > PI_OVER_2))
206
{ /* First Standard Parallel out of range */
207
Error_Code |= ALBERS_FIRST_STDP_ERROR;
209
if ((Std_Parallel_2 < -PI_OVER_2) || (Std_Parallel_2 > PI_OVER_2))
210
{ /* Second Standard Parallel out of range */
211
Error_Code |= ALBERS_SECOND_STDP_ERROR;
213
if ((Std_Parallel_1 == 0.0) && (Std_Parallel_2 == 0.0))
214
{ /* First & Second Standard Parallels equal 0 */
215
Error_Code |= ALBERS_FIRST_SECOND_ERROR;
217
if (Std_Parallel_1 == -Std_Parallel_2)
218
{ /* Parallels are opposite latitudes */
219
Error_Code |= ALBERS_HEMISPHERE_ERROR;
226
Albers_Origin_Lat = Origin_Latitude;
227
Albers_Std_Parallel_1 = Std_Parallel_1;
228
Albers_Std_Parallel_2 = Std_Parallel_2;
229
if (Central_Meridian > PI)
230
Central_Meridian -= TWO_PI;
231
Albers_Origin_Long = Central_Meridian;
232
Albers_False_Easting = False_Easting;
233
Albers_False_Northing = False_Northing;
235
es2 = 2 * Albers_f - Albers_f * Albers_f;
237
one_MINUS_es2 = 1 - es2;
240
sin_lat = sin(Albers_Origin_Lat);
241
es_sin = ES_SIN(sin_lat);
242
one_MINUS_SQRes_sin = ONE_MINUS_SQR(es_sin);
243
q0 = ALBERS_Q(sin_lat, one_MINUS_SQRes_sin, es_sin);
245
sin_lat_1 = sin(Albers_Std_Parallel_1);
246
cos_lat = cos(Albers_Std_Parallel_1);
247
es_sin = ES_SIN(sin_lat_1);
248
one_MINUS_SQRes_sin = ONE_MINUS_SQR(es_sin);
249
m1 = ALBERS_M(cos_lat, one_MINUS_SQRes_sin);
250
q1 = ALBERS_Q(sin_lat_1, one_MINUS_SQRes_sin, es_sin);
253
if (fabs(Albers_Std_Parallel_1 - Albers_Std_Parallel_2) > 1.0e-10)
255
sin_lat = sin(Albers_Std_Parallel_2);
256
cos_lat = cos(Albers_Std_Parallel_2);
257
es_sin = ES_SIN(sin_lat);
258
one_MINUS_SQRes_sin = ONE_MINUS_SQR(es_sin);
259
m2 = ALBERS_M(cos_lat, one_MINUS_SQRes_sin);
260
q2 = ALBERS_Q(sin_lat, one_MINUS_SQRes_sin, es_sin);
261
n = (SQRm1 - m2 * m2) / (q2 - q1);
267
Albers_a_OVER_n = Albers_a / n;
272
rho0 = Albers_a_OVER_n * sqrt(C - nq0);
275
} /* END OF if(!Error_Code) */
277
} /* END OF Set_Albers_Parameters */
280
void Get_Albers_Parameters(double *a,
282
double *Origin_Latitude,
283
double *Central_Meridian,
284
double *Std_Parallel_1,
285
double *Std_Parallel_2,
286
double *False_Easting,
287
double *False_Northing)
289
{ /* BEGIN Get_Albers_Parameters */
291
* The function Get_Albers_Parameters returns the current ellipsoid
292
* parameters, and Albers projection parameters.
294
* a : Semi-major axis of ellipsoid, in meters (output)
295
* f : Flattening of ellipsoid (output)
296
* Origin_Latitude : Latitude in radians at which the (output)
297
* point scale factor is 1.0
298
* Origin_Longitude : Longitude in radians at the center of (output)
300
* Std_Parallel_1 : First standard parallel (output)
301
* Std_Parallel_2 : Second standard parallel (output)
302
* False_Easting : A coordinate value in meters assigned to the
303
* central meridian of the projection. (output)
304
* False_Northing : A coordinate value in meters assigned to the
305
* origin latitude of the projection (output)
310
*Origin_Latitude = Albers_Origin_Lat;
311
*Std_Parallel_1 = Albers_Std_Parallel_1;
312
*Std_Parallel_2 = Albers_Std_Parallel_2;
313
*Central_Meridian = Albers_Origin_Long;
314
*False_Easting = Albers_False_Easting;
315
*False_Northing = Albers_False_Northing;
317
} /* END OF Get_Albers_Parameters */
320
long Convert_Geodetic_To_Albers (double Latitude,
325
{ /* BEGIN Convert_Geodetic_To_Albers */
327
* The function Convert_Geodetic_To_Albers converts geodetic (latitude and
328
* longitude) coordinates to Albers projection (easting and northing)
329
* coordinates, according to the current ellipsoid and Albers projection
330
* parameters. If any errors occur, the error code(s) are returned by the
331
* function, otherwise ALBERS_NO_ERROR is returned.
333
* Latitude : Latitude (phi) in radians (input)
334
* Longitude : Longitude (lambda) in radians (input)
335
* Easting : Easting (X) in meters (output)
336
* Northing : Northing (Y) in meters (output)
339
double dlam; /* Longitude - Central Meridan */
340
double sin_lat, cos_lat;
341
double es_sin, one_MINUS_SQRes_sin;
346
long Error_Code = ALBERS_NO_ERROR;
348
if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
349
{ /* Latitude out of range */
350
Error_Code |= ALBERS_LAT_ERROR;
352
if ((Longitude < -PI) || (Longitude > TWO_PI))
353
{ /* Longitude out of range */
354
Error_Code|= ALBERS_LON_ERROR;
360
dlam = Longitude - Albers_Origin_Long;
369
sin_lat = sin(Latitude);
370
cos_lat = cos(Latitude);
371
es_sin = ES_SIN(sin_lat);
372
one_MINUS_SQRes_sin = ONE_MINUS_SQR(es_sin);
373
q = ALBERS_Q(sin_lat, one_MINUS_SQRes_sin, es_sin);
378
rho = Albers_a_OVER_n * sqrt(C - nq);
382
*Easting = rho * sin(theta) + Albers_False_Easting;
383
*Northing = rho0 - rho * cos(theta) + Albers_False_Northing;
386
} /* END OF Convert_Geodetic_To_Albers */
389
long Convert_Albers_To_Geodetic(double Easting,
393
{ /* BEGIN Convert_Albers_To_Geodetic */
395
* The function Convert_Albers_To_Geodetic converts Albers projection
396
* (easting and northing) coordinates to geodetic (latitude and longitude)
397
* coordinates, according to the current ellipsoid and Albers projection
398
* coordinates. If any errors occur, the error code(s) are returned by the
399
* function, otherwise ALBERS_NO_ERROR is returned.
401
* Easting : Easting (X) in meters (input)
402
* Northing : Northing (Y) in meters (input)
403
* Latitude : Latitude (phi) in radians (output)
404
* Longitude : Longitude (lambda) in radians (output)
407
double rho0_MINUS_dy;
408
double q, qconst, q_OVER_2;
410
double PHI, Delta_PHI = 1.0;
412
double es_sin, one_MINUS_SQRes_sin;
415
double tolerance = 4.85e-10; /* approximately 1/1000th of
416
an arc second or 1/10th meter */
417
long Error_Code = ALBERS_NO_ERROR;
419
if ((Easting < (Albers_False_Easting - Albers_Delta_Easting))
420
|| (Easting > Albers_False_Easting + Albers_Delta_Easting))
421
{ /* Easting out of range */
422
Error_Code |= ALBERS_EASTING_ERROR;
424
if ((Northing < (Albers_False_Northing - Albers_Delta_Northing))
425
|| (Northing > Albers_False_Northing + Albers_Delta_Northing))
426
{ /* Northing out of range */
427
Error_Code |= ALBERS_NORTHING_ERROR;
432
dy = Northing - Albers_False_Northing;
433
dx = Easting - Albers_False_Easting;
434
rho0_MINUS_dy = rho0 - dy;
435
rho = sqrt(dx * dx + rho0_MINUS_dy * rho0_MINUS_dy);
442
rho0_MINUS_dy *= -1.0;
446
theta = atan2(dx, rho0_MINUS_dy);
448
q = (C - (rho_n * rho_n) / (Albers_a * Albers_a)) / n;
449
qconst = 1 - ((one_MINUS_es2) / (two_es)) * log((1.0 - es) / (1.0 + es));
450
if (fabs(fabs(qconst) - fabs(q)) > 1.0e-6)
454
*Latitude = PI_OVER_2;
455
else if (q_OVER_2 < -1.0)
456
*Latitude = -PI_OVER_2;
459
PHI = asin(q_OVER_2);
464
while ((fabs(Delta_PHI) > tolerance) && count)
467
es_sin = ES_SIN(sin_phi);
468
one_MINUS_SQRes_sin = ONE_MINUS_SQR(es_sin);
469
Delta_PHI = (one_MINUS_SQRes_sin * one_MINUS_SQRes_sin) / (2.0 * cos(PHI)) *
470
(q / (one_MINUS_es2) - sin_phi / one_MINUS_SQRes_sin +
471
(log((1.0 - es_sin) / (1.0 + es_sin)) / (two_es)));
477
return Error_Code |= ALBERS_NORTHING_ERROR;
482
if (*Latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
483
*Latitude = PI_OVER_2;
484
else if (*Latitude < -PI_OVER_2)
485
*Latitude = -PI_OVER_2;
492
*Latitude = PI_OVER_2;
494
*Latitude = -PI_OVER_2;
496
*Longitude = Albers_Origin_Long + theta / n;
499
*Longitude -= TWO_PI;
500
if (*Longitude < -PI)
501
*Longitude += TWO_PI;
503
if (*Longitude > PI) /* force distorted values to 180, -180 degrees */
505
else if (*Longitude < -PI)
510
} /* END OF Convert_Albers_To_Geodetic */
added
removed
dt_cc/albers/albers.c
