OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
oblique_mercator.c
Go to the documentation of this file.
1 /*******************************************************************************
2 Name: OBLIQUE MERCATOR (HOTINE)
3 
4 Purpose: Provides the transformations between Easting/Northing and longitude/
5  latitude for the Oblique Mercator projection. The Easting and Northing
6  values are in meters. The longitude and latitude are in radians.
7 
8 Algorithm References
9 
10 1. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
11  Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
12  State Government Printing Office, Washington D.C., 1987.
13 
14 2. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
15  U.S. Geological Survey Professional Paper 1453 , United State Government
16  Printing Office, Washington D.C., 1989.
17 *******************************************************************************/
18 #include <stdlib.h>
19 #include "oli_cproj.h"
20 #include "gctp.h"
21 #include "oli_local.h"
22 
23 /* structure to hold the setup data relevant to this projection */
24 struct om_proj
25 {
26  double r_major; /* major axis */
27  double r_minor; /* minor axis */
28  double scale_factor; /* scale factor */
29  double lon_origin; /* center longitude */
30  double lat_origin; /* center latitude */
31  double lon1; /* first point to define central line */
32  double lat1; /* first point to define central line */
33  double lon2; /* second point to define central line */
34  double lat2; /* second point to define central line */
35  double azimuth; /* azimuth east of north */
36  double false_northing; /* y offset in meters */
37  double false_easting; /* x offset in meters */
38  double e; /* eccentricity constants */
39  double es;
40  double sin_p20; /* sin value */
41  double cos_p20; /* cos values */
42  double bl;
43  double al;
44  double d;
45  double el;
46  double u;
47  double singam;
48  double cosgam;
49  double sinaz;
50  double cosaz;
51  int mode; /* format type 0 = A, 1 = B */
52 };
53 
54 /*****************************************************************************
55 Name: print_info
56 
57 Purpose: Prints a summary of information about this projection.
58 
59 Returns:
60  nothing
61 
62 *****************************************************************************/
63 static void print_info
64 (
65  const TRANSFORMATION *trans
66 )
67 {
68  struct om_proj *cache_ptr = (struct om_proj *)trans->cache;
69 
70  gctp_print_title("OBLIQUE MERCATOR (HOTINE)");
71  gctp_print_radius2(cache_ptr->r_major, cache_ptr->r_minor);
73  "Scale Factor at C. Meridian: ");
74  gctp_print_offsetp(cache_ptr->false_easting, cache_ptr->false_northing);
75 
76  if (cache_ptr->mode == 0)
77  {
78  /* Format A parameters */
79  gctp_print_genrpt(cache_ptr->lon1 * R2D,
80  "Longitude of First Point: ");
81  gctp_print_genrpt(cache_ptr->lat1 * R2D,
82  "Latitude of First Point: ");
83  gctp_print_genrpt(cache_ptr->lon2 * R2D,
84  "Longitude of Second Point: ");
85  gctp_print_genrpt(cache_ptr->lat2 * R2D,
86  "Latitude of Second Point: ");
87  }
88  else
89  {
90  /* Format B parameters */
91  gctp_print_genrpt(cache_ptr->azimuth * R2D,
92  "Azimuth of Central Line: ");
93  gctp_print_cenlon(cache_ptr->lon_origin);
94  gctp_print_cenlat(cache_ptr->lat_origin);
95 
96  }
97 }
98 
99 /*******************************************************************************
100 Name: common_init
101 
102 Purpose: Initialization routine for initializing the projection information
103  that is common to both the forward and inverse transformations.
104 
105 Returns:
106  GCTP_SUCCESS or GCTP_ERROR
107 
108 *******************************************************************************/
109 static int common_init
110 (
111  TRANSFORMATION *trans /* I/O: transformation to initialize */
112 )
113 {
114  double r_major; /* major axis */
115  double r_minor; /* minor axis */
116  double radius;
117  double lon_origin; /* center longitude */
118  double lat_origin; /* center latitude */
119  double false_easting; /* x offset in meters */
120  double false_northing; /* y offset in meters */
121  double scale_factor; /* scale factor */
122  double azimuth; /* azimuth east of north */
123  double lon1; /* first point to define central line */
124  double lat1; /* first point to define central line */
125  double lon2; /* second point to define central line */
126  double lat2; /* second point to define central line */
127  int mode; /* format A/B indicator */
128  double temp; /* temporary variable */
129  double sinphi; /* sin value */
130  double con,com;
131  double ts;
132  double ts1,ts2;
133  double h,l;
134  double j,p,dlon;
135  double f = 0.0;
136  double g,gama;
137 
138  const GCTP_PROJECTION *proj = &trans->proj;
139  struct om_proj *cache = NULL;
140  int spheroid = proj->spheroid;
141 
142  gctp_get_spheroid(spheroid, proj->parameters, &r_major, &r_minor, &radius);
143  scale_factor = proj->parameters[2];
144  false_easting = proj->parameters[6];
145  false_northing = proj->parameters[7];
146 
147  if (gctp_dms2degrees(proj->parameters[5], &lat_origin) != GCTP_SUCCESS)
148  {
149  GCTP_PRINT_ERROR("Error converting center latitude in parameter 4 "
150  "from DMS to degrees: %f", proj->parameters[5]);
151  return GCTP_ERROR;
152  }
153  lat_origin *= 3600 * S2R;
154 
155  if (proj->parameters[12] != 0)
156  mode = 1;
157  else
158  mode = 0;
159 
160  if (mode == 0)
161  {
162  if (gctp_dms2degrees(proj->parameters[8], &lon1) != GCTP_SUCCESS)
163  {
164  GCTP_PRINT_ERROR("Error converting first point Long1 in parameter "
165  "8 from DMS to degrees: %f", proj->parameters[8]);
166  return GCTP_ERROR;
167  }
168  lon1 *= 3600 * S2R;
169 
170  if (gctp_dms2degrees(proj->parameters[9], &lat1) != GCTP_SUCCESS)
171  {
172  GCTP_PRINT_ERROR("Error converting first point Lat1 in parameter "
173  "9 from DMS to degrees: %f", proj->parameters[9]);
174  return GCTP_ERROR;
175  }
176  lat1 *= 3600 * S2R;
177 
178  if (gctp_dms2degrees(proj->parameters[10], &lon2) != GCTP_SUCCESS)
179  {
180  GCTP_PRINT_ERROR("Error converting first point Long2 in parameter "
181  "10 from DMS to degrees: %f", proj->parameters[10]);
182  return GCTP_ERROR;
183  }
184  lon2 *= 3600 * S2R;
185 
186  if (gctp_dms2degrees(proj->parameters[11], &lat2) != GCTP_SUCCESS)
187  {
188  GCTP_PRINT_ERROR("Error converting first point Lat2 in parameter "
189  "11 from DMS to degrees: %f", proj->parameters[11]);
190  return GCTP_ERROR;
191  }
192  lat2 *= 3600 * S2R;
193 
194  azimuth = 0.0;
195  lon_origin = 0.0;
196  }
197  else
198  {
199  if (gctp_dms2degrees(proj->parameters[3], &azimuth) != GCTP_SUCCESS)
200  {
201  GCTP_PRINT_ERROR("Error converting center azimuth angle in "
202  "parameter 3 from DMS to degrees: %f", proj->parameters[3]);
203  return GCTP_ERROR;
204  }
205  azimuth *= 3600 * S2R;
206 
207  if (gctp_dms2degrees(proj->parameters[4], &lon_origin) != GCTP_SUCCESS)
208  {
209  GCTP_PRINT_ERROR("Error converting center azimuth point in "
210  "parameter 4 from DMS to degrees: %f", proj->parameters[4]);
211  return GCTP_ERROR;
212  }
213  lon_origin *= 3600 * S2R;
214 
215  lon1 = 0.0;
216  lat1 = 0.0;
217  lon2 = 0.0;
218  lat2 = 0.0;
219  }
220 
221  /* Allocate a structure for the cached info */
222  cache = malloc(sizeof(*cache));
223  if (!cache)
224  {
225  GCTP_PRINT_ERROR("Error allocating memory for cache buffer");
226  return GCTP_ERROR;
227  }
228  trans->cache = cache;
229 
230  /* Save the information to the cache */
231  cache->r_major = r_major;
232  cache->r_minor = r_minor;
233  cache->false_northing = false_northing;
234  cache->false_easting = false_easting;
235  cache->scale_factor = scale_factor;
236  cache->lat_origin = lat_origin;
237  cache->lon1 = lon1;
238  cache->lat1 = lat1;
239  cache->lon2 = lon2;
240  cache->lat2 = lat2;
241  cache->mode = mode;
242 
243  temp = r_minor / r_major;
244  cache->es = 1.0 - SQUARE(temp);
245  cache->e = sqrt(cache->es);
246 
247  sincos(lat_origin, &cache->sin_p20, &cache->cos_p20);
248  con = 1.0 - cache->es * cache->sin_p20 * cache->sin_p20;
249  com = sqrt(1.0 - cache->es);
250  cache->bl = sqrt(1.0 + cache->es * pow(cache->cos_p20, 4.0)
251  / (1.0 - cache->es));
252  cache->al = r_major * cache->bl * cache->scale_factor * com / con;
253  if (fabs(lat_origin) < EPSLN)
254  {
255  ts = 1.0;
256  cache->d = 1.0;
257  cache->el = 1.0;
258  }
259  else
260  {
261  ts = gctp_calc_small_t(cache->e, lat_origin, cache->sin_p20);
262  con = sqrt(con);
263  cache->d = cache->bl * com / (cache->cos_p20 * con);
264  if ((cache->d * cache->d - 1.0) > 0.0)
265  {
266  if (lat_origin >= 0.0)
267  f = cache->d + sqrt(cache->d * cache->d - 1.0);
268  else
269  f = cache->d - sqrt(cache->d * cache->d - 1.0);
270  }
271  else
272  f = cache->d;
273  cache->el = f * pow(ts, cache->bl);
274  }
275 
276  if (mode != 0)
277  {
278  g = .5 * (f - 1.0/f);
279  gama = asinz(sin(azimuth) / cache->d);
280  lon_origin = lon_origin - asinz(g * tan(gama))/cache->bl;
281 
282  con = fabs(lat_origin);
283  if ((con > EPSLN) && (fabs(con - HALF_PI) > EPSLN))
284  {
285  sincos(gama, &cache->singam, &cache->cosgam);
286  sincos(azimuth, &cache->sinaz, &cache->cosaz);
287  if (lat_origin >= 0)
288  {
289  cache->u = (cache->al / cache->bl)
290  * atan(sqrt(cache->d * cache->d - 1.0) / cache->cosaz);
291  }
292  else
293  {
294  cache->u = -(cache->al / cache->bl)
295  * atan(sqrt(cache->d * cache->d - 1.0) / cache->cosaz);
296  }
297  }
298  else
299  {
300  GCTP_PRINT_ERROR("Input data error");
301  free(cache);
302  trans->cache = NULL;
303  return GCTP_ERROR;
304  }
305  }
306  else
307  {
308  sinphi = sin(lat1);
309  ts1 = gctp_calc_small_t(cache->e, lat1, sinphi);
310  sinphi = sin(lat2);
311  ts2 = gctp_calc_small_t(cache->e, lat2, sinphi);
312  h = pow(ts1, cache->bl);
313  l = pow(ts2, cache->bl);
314  f = cache->el / h;
315  g = .5 * (f - 1.0/f);
316  j = (cache->el * cache->el - l * h)/(cache->el * cache->el + l * h);
317  p = (l - h) / (l + h);
318  dlon = lon1 - lon2;
319  if (dlon < -PI)
320  lon2 = lon2 - 2.0 * PI;
321  if (dlon > PI)
322  lon2 = lon2 + 2.0 * PI;
323  dlon = lon1 - lon2;
324  lon_origin = .5 * (lon1 + lon2)
325  - atan(j * tan(.5 * cache->bl * dlon) / p) / cache->bl;
326  dlon = adjust_lon(lon1 - lon_origin);
327  gama = atan(sin(cache->bl * dlon) / g);
328  azimuth = asinz(cache->d * sin(gama));
329 
330  if (fabs(lat1 - lat2) <= EPSLN)
331  {
332  GCTP_PRINT_ERROR("Input data error");
333  free(cache);
334  trans->cache = NULL;
335  return GCTP_ERROR;
336  }
337  else
338  con = fabs(lat1);
339  if ((con <= EPSLN) || (fabs(con - HALF_PI) <= EPSLN))
340  {
341  GCTP_PRINT_ERROR("Input data error");
342  free(cache);
343  trans->cache = NULL;
344  return GCTP_ERROR;
345  }
346  else if (fabs(fabs(lat_origin) - HALF_PI) <= EPSLN)
347  {
348  GCTP_PRINT_ERROR("Input data error");
349  free(cache);
350  trans->cache = NULL;
351  return GCTP_ERROR;
352  }
353 
354  sincos(gama, &cache->singam, &cache->cosgam);
355  sincos(azimuth, &cache->sinaz, &cache->cosaz);
356  if (lat_origin >= 0)
357  {
358  cache->u = (cache->al / cache->bl)
359  * atan(sqrt(cache->d * cache->d - 1.0)/cache->cosaz);
360  }
361  else
362  {
363  cache->u = -(cache->al / cache->bl)
364  * atan(sqrt(cache->d * cache->d - 1.0)/cache->cosaz);
365  }
366  }
367 
368  cache->lon_origin = lon_origin;
369  cache->azimuth = azimuth;
370 
371  trans->print_info = print_info;
372 
373  return GCTP_SUCCESS;
374 }
375 
376 /*****************************************************************************
377 Name: inverse_transform
378 
379 Purpose: Transforms X,Y to lat,long
380 
381 Returns:
382  GCTP_SUCCESS or GCTP_ERROR
383 
384 *****************************************************************************/
385 static int inverse_transform
386 (
387  const TRANSFORMATION *trans, /* I: transformation information */
388  double x, /* I: X projection coordinate */
389  double y, /* I: Y projection coordinate */
390  double *lon, /* O: Longitude */
391  double *lat /* O: Latitude */
392 )
393 {
394  struct om_proj *cache_ptr = (struct om_proj *)trans->cache;
395  double theta; /* angle */
396  double t; /* temporary values */
397  double con; /* cone constant, small m */
398  double vs,us,q,s,ts1;
399  double vl,ul;
400 
401  x -= cache_ptr->false_easting;
402  y -= cache_ptr->false_northing;
403  vs = x * cache_ptr->cosaz - y * cache_ptr->sinaz;
404  us = y * cache_ptr->cosaz + x * cache_ptr->sinaz;
405  us = us + cache_ptr->u;
406  q = exp(-cache_ptr->bl * vs / cache_ptr->al);
407  s = .5 * (q - 1.0/q);
408  t = .5 * (q + 1.0/q);
409  vl = sin(cache_ptr->bl * us / cache_ptr->al);
410  ul = (vl * cache_ptr->cosgam + s * cache_ptr->singam)/t;
411  if (fabs(fabs(ul) - 1.0) <= EPSLN)
412  {
413  *lon = cache_ptr->lon_origin;
414  if (ul >= 0.0)
415  *lat = HALF_PI;
416  else
417  *lat = -HALF_PI;
418  }
419  else
420  {
421  con = 1.0 / cache_ptr->bl;
422  ts1 = pow((cache_ptr->el / sqrt((1.0 + ul) / (1.0 - ul))), con);
423  if (gctp_calc_phi2(cache_ptr->e, ts1, lat) != GCTP_SUCCESS)
424  return GCTP_ERROR;
425  con = cos(cache_ptr->bl * us /cache_ptr->al);
426  theta = cache_ptr->lon_origin - atan2((s * cache_ptr->cosgam - vl
427  * cache_ptr->singam), con) / cache_ptr->bl;
428  *lon = adjust_lon(theta);
429  }
430 
431  return GCTP_SUCCESS;
432 }
433 
434 /*****************************************************************************
435 Name: forward_transform
436 
437 Purpose: Transforms lat,long to oblique mercator X,Y
438 
439 Returns:
440  GCTP_SUCCESS or GCTP_ERROR
441 
442 *****************************************************************************/
443 static int forward_transform
444 (
445  const TRANSFORMATION *trans, /* I: transformation information */
446  double lon, /* I: Longitude */
447  double lat, /* I: Latitude */
448  double *x, /* O: X projection coordinate */
449  double *y /* O: Y projection coordinate */
450 )
451 {
452  struct om_proj *cache_ptr = (struct om_proj *)trans->cache;
453  double sin_phi; /* sin value */
454  double t; /* temporary values */
455  double con; /* cone constant, small m */
456  double q,us,vl;
457  double ul,vs;
458  double s;
459  double dlon;
460  double ts1;
461 
462  sin_phi = sin(lat);
463  dlon = adjust_lon(lon - cache_ptr->lon_origin);
464  vl = sin(cache_ptr->bl * dlon);
465  if (fabs(fabs(lat) - HALF_PI) > EPSLN)
466  {
467  ts1 = gctp_calc_small_t(cache_ptr->e, lat, sin_phi);
468  q = cache_ptr->el / (pow(ts1, cache_ptr->bl));
469  s = .5 * (q - 1.0 / q);
470  t = .5 * (q + 1.0/ q);
471  ul = (s * cache_ptr->singam - vl * cache_ptr->cosgam) / t;
472  con = cos(cache_ptr->bl * dlon);
473  if (fabs(con) < .0000001)
474  {
475  us = cache_ptr->al * dlon;
476  }
477  else
478  {
479  us = cache_ptr->al * atan((s * cache_ptr->cosgam + vl
480  * cache_ptr->singam) / con)/cache_ptr->bl;
481  if (con < 0)
482  us = us + PI * cache_ptr->al / cache_ptr->bl;
483  }
484  }
485  else
486  {
487  if (lat >= 0)
488  ul = cache_ptr->singam;
489  else
490  ul = -cache_ptr->singam;
491  us = cache_ptr->al * lat / cache_ptr->bl;
492  }
493  if (fabs(fabs(ul) - 1.0) <= EPSLN)
494  {
495  GCTP_PRINT_ERROR("Point projects into infinity");
496  return GCTP_ERROR;
497  }
498  vs = .5 * cache_ptr->al * log((1.0 - ul)/(1.0 + ul)) / cache_ptr->bl;
499  us = us - cache_ptr->u;
500  *x = cache_ptr->false_easting + vs * cache_ptr->cosaz
501  + us * cache_ptr->sinaz;
502  *y = cache_ptr->false_northing + us * cache_ptr->cosaz
503  - vs * cache_ptr->sinaz;
504 
505  return GCTP_SUCCESS;
506 }
507 
508 
509 /*****************************************************************************
510 Name: gctp_om_inverse_init
511 
512 Purpose: Initializes the inverse oblique mercator transformation
513 
514 Returns:
515  GCTP_SUCCESS or GCTP_ERROR
516 
517 *****************************************************************************/
519 (
520  TRANSFORMATION *trans
521 )
522 {
523  /* Call the common routine used for the forward and inverse init */
524  if (common_init(trans) != GCTP_SUCCESS)
525  {
527  "Error initializing oblique mercator inverse projection");
528  return GCTP_ERROR;
529  }
530 
531  trans->transform = inverse_transform;
532 
533  return GCTP_SUCCESS;
534 }
535 
536 /*****************************************************************************
537 Name: gctp_om_forward_init
538 
539 Purpose: Initializes the forward oblique mercator transformation
540 
541 Returns:
542  GCTP_SUCCESS or GCTP_ERROR
543 
544 *****************************************************************************/
546 (
547  TRANSFORMATION *trans
548 )
549 {
550  /* Call the common routine used for the forward and inverse init */
551  if (common_init(trans) != GCTP_SUCCESS)
552  {
554  "Error initializing oblique mercator forward projection");
555  return GCTP_ERROR;
556  }
557 
558  trans->transform = forward_transform;
559 
560  return GCTP_SUCCESS;
561 }
data_t q
Definition: decode_rs.h:74
#define SQUARE(x)
Definition: proj_define.h:99
data_t t[NROOTS+1]
Definition: decode_rs.h:77
double lon_origin
double lat_origin
double sinaz
int j
Definition: decode_rs.h:73
double singam
void gctp_print_title(const char *proj_name)
Definition: gctp_report.c:14
#define GCTP_ERROR
Definition: gctp.h:81
double sin_p20
#define NULL
Definition: decode_rs.h:63
void gctp_print_cenlon(double A)
Definition: gctp_report.c:40
README for MOD_PR02AQUA(AQUA) Version to set to For disabling creating and output data sets when in night mode
Definition: README.txt:96
#define GCTP_PRINT_ERROR(format,...)
Definition: oli_local.h:81
float h[MODELMAX]
Definition: atrem_corl1.h:131
float * lat
double cosgam
double adjust_lon(double x)
Definition: proj_cproj.c:349
#define PI
Definition: l3_get_org.c:6
double r_minor
double precision function f(R1)
Definition: tmd.lp.f:1454
def cache(filename, recache=False)
Definition: utils.py:145
int gctp_calc_phi2(double eccent, double ts, double *phi2)
Definition: gctp_utility.c:209
int gctp_om_forward_init(TRANSFORMATION *trans)
#define HALF_PI
Definition: proj_define.h:84
double gctp_calc_small_t(double eccent, double phi, double sinphi)
Definition: gctp_utility.c:277
double r_major
void gctp_print_offsetp(double A, double B)
Definition: gctp_report.c:91
double lat2
#define GCTP_SUCCESS
Definition: gctp.h:82
double false_easting
int gctp_dms2degrees(double angle, double *degrees)
void gctp_print_genrpt(double A, const char *S)
Definition: gctp_report.c:117
#define gctp_get_spheroid
Definition: oli_local.h:7
double scale_factor
#define sincos
Definition: proj_define.h:108
int gctp_om_inverse_init(TRANSFORMATION *trans)
void gctp_print_radius2(double radius1, double radius2)
Definition: gctp_report.c:30
#define fabs(a)
Definition: misc.h:93
double lat1
#define S2R
Definition: proj_define.h:92
#define R2D
Definition: proj_define.h:87
double azimuth
float * lon
data_t s[NROOTS]
Definition: decode_rs.h:75
double cos_p20
double cosaz
void radius(double A)
Definition: proj_report.c:132
double lon1
double false_northing
double lon2
double asinz(double con)
Definition: proj_cproj.c:67
float p[MODELMAX]
Definition: atrem_corl1.h:131
void gctp_print_cenlat(double A)
Definition: gctp_report.c:57
#define EPSLN
Definition: proj_define.h:86