OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
lambert_conformal_conic.c
Go to the documentation of this file.
1 /*******************************************************************************
2 Name: LAMBERT CONFORMAL CONIC
3 
4 Purpose: Provides the transformations between Easting/Northing and longitude/
5  latitude for the Lambert Conformal Conic projection. The Easting and
6  Northing 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 *******************************************************************************/
17 #include <stdlib.h>
18 #include <math.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 lamcc_proj
25 {
26  double r_major; /* major axis */
27  double r_minor; /* minor axis */
28  double es; /* eccentricity squared */
29  double e; /* eccentricity */
30  double lat1; /* first standard parallel */
31  double lat2; /* second standard parallel */
32  double center_lon; /* center longitude */
33  double lat_origin; /* latitude of origin */
34  double ns; /* ratio of angle between meridian */
35  double f0; /* flattening of ellipsoid */
36  double rh; /* height above ellipsoid */
37  double false_easting; /* x offset in meters */
38  double false_northing; /* y offset in meters */
39 };
40 
41 /*****************************************************************************
42 Name: print_info
43 
44 Purpose: Prints a summary of information about this projection.
45 
46 Returns:
47  nothing
48 
49 *****************************************************************************/
50 static void print_info
51 (
52  const TRANSFORMATION *trans
53 )
54 {
55  struct lamcc_proj *cache_ptr = (struct lamcc_proj *)trans->cache;
56 
57  gctp_print_title("LAMBERT CONFORMAL CONIC");
58  gctp_print_radius2(cache_ptr->r_major, cache_ptr->r_minor);
59  gctp_print_stanparl(cache_ptr->lat1, cache_ptr->lat2);
60  gctp_print_cenlonmer(cache_ptr->center_lon);
61  gctp_print_origin(cache_ptr->lat_origin);
62  gctp_print_offsetp(cache_ptr->false_easting, cache_ptr->false_northing);
63 }
64 
65 /*******************************************************************************
66 Name: common_init
67 
68 Purpose: Initialization routine for initializing the projection information
69  that is common to both the forward and inverse transformations.
70 
71 Returns:
72  GCTP_SUCCESS or GCTP_ERROR
73 
74 *******************************************************************************/
75 static int common_init
76 (
77  TRANSFORMATION *trans /* I/O: transformation to initialize */
78 )
79 {
80  double r_major; /* major axis */
81  double r_minor; /* minor axis */
82  double radius; /* earth radius */
83  double lat1; /* first standard parallel */
84  double lat2; /* second standard parallel */
85  double center_lon; /* center longitude */
86  double lat_origin; /* center latitude */
87  double false_easting; /* x offset in meters */
88  double false_northing; /* y offset in meters */
89  double sin_po; /* sin value */
90  double cos_po; /* cos value */
91  double con; /* temporary variable */
92  double ms1; /* small m 1 */
93  double ms2; /* small m 2 */
94  double temp; /* temporary variable */
95  double ts0; /* small t 0 */
96  double ts1; /* small t 1 */
97  double ts2; /* small t 2 */
98 
99  const GCTP_PROJECTION *proj = &trans->proj;
100  struct lamcc_proj *cache = NULL;
101  int spheroid = proj->spheroid;
102 
103  gctp_get_spheroid(spheroid, proj->parameters, &r_major, &r_minor, &radius);
104  false_easting = proj->parameters[6];
105  false_northing = proj->parameters[7];
106 
107  if (gctp_dms2degrees(proj->parameters[2], &lat1) != GCTP_SUCCESS)
108  {
109  GCTP_PRINT_ERROR("Error converting standard parallel 1 in parameter 2 "
110  "from DMS to degrees: %f", proj->parameters[2]);
111  return GCTP_ERROR;
112  }
113  lat1 *= 3600 * S2R;
114 
115  if (gctp_dms2degrees(proj->parameters[3], &lat2) != GCTP_SUCCESS)
116  {
117  GCTP_PRINT_ERROR("Error converting standard parallel 2 in parameter 3 "
118  "from DMS to degrees: %f", proj->parameters[3]);
119  return GCTP_ERROR;
120  }
121  lat2 *= 3600 * S2R;
122 
123  if (gctp_dms2degrees(proj->parameters[4], &center_lon) != GCTP_SUCCESS)
124  {
125  GCTP_PRINT_ERROR("Error converting center longitude in parameter 4 "
126  "from DMS to degrees: %f", proj->parameters[4]);
127  return GCTP_ERROR;
128  }
129  center_lon *= 3600 * S2R;
130 
131  if (gctp_dms2degrees(proj->parameters[5], &lat_origin) != GCTP_SUCCESS)
132  {
133  GCTP_PRINT_ERROR("Error converting center latitude in parameter 4 "
134  "from DMS to degrees: %f", proj->parameters[5]);
135  return GCTP_ERROR;
136  }
137  lat_origin *= 3600 * S2R;
138 
139  /* Standard Parallels cannot be equal and on opposite sides of the
140  equator */
141  if (fabs(lat1 + lat2) < EPSLN)
142  {
143  GCTP_PRINT_ERROR("Equal latitudes for Standard Parallels on opposite "
144  "sides of equator");
145  return GCTP_ERROR;
146  }
147 
148  /* Allocate a structure for the cached info */
149  cache = malloc(sizeof(*cache));
150  if (!cache)
151  {
152  GCTP_PRINT_ERROR("Error allocating memory for cache buffer");
153  return GCTP_ERROR;
154  }
155  trans->cache = cache;
156 
157  /* Save the information to the cache */
158  cache->r_major = r_major;
159  cache->r_minor = r_minor;
160  cache->false_northing = false_northing;
161  cache->false_easting = false_easting;
162  cache->lat1 = lat1;
163  cache->lat2 = lat2;
164  cache->center_lon = center_lon;
165  cache->lat_origin = lat_origin;
166 
167  temp = r_minor / r_major;
168  cache->es = 1.0 - SQUARE(temp);
169  cache->e = sqrt(cache->es);
170 
171  sincos(lat1, &sin_po, &cos_po);
172  con = sin_po;
173  ms1 = gctp_calc_small_radius(cache->e, sin_po, cos_po);
174  ts1 = gctp_calc_small_t(cache->e, lat1, sin_po);
175  sincos(lat2, &sin_po, &cos_po);
176  ms2 = gctp_calc_small_radius(cache->e, sin_po, cos_po);
177  ts2 = gctp_calc_small_t(cache->e, lat2, sin_po);
178  sin_po = sin(lat_origin);
179  ts0 = gctp_calc_small_t(cache->e, lat_origin, sin_po);
180 
181  if (fabs(lat1 - lat2) > EPSLN)
182  cache->ns = log(ms1/ms2)/ log(ts1/ts2);
183  else
184  cache->ns = con;
185  cache->f0 = ms1 / (cache->ns * pow(ts1, cache->ns));
186  cache->rh = r_major * cache->f0 * pow(ts0, cache->ns);
187 
188  trans->print_info = print_info;
189 
190  return GCTP_SUCCESS;
191 }
192 
193 /*****************************************************************************
194 Name: inverse_transform
195 
196 Purpose: Transforms X,Y to lat,long
197 
198 Returns:
199  GCTP_SUCCESS or GCTP_ERROR
200 
201 *****************************************************************************/
202 static int inverse_transform
203 (
204  const TRANSFORMATION *trans, /* I: transformation information */
205  double x, /* I: X projection coordinate */
206  double y, /* I: Y projection coordinate */
207  double *lon, /* O: Longitude */
208  double *lat /* O: Latitude */
209 )
210 {
211  struct lamcc_proj *cache_ptr = (struct lamcc_proj *)trans->cache;
212  double rh1; /* height above ellipsoid */
213  double con; /* sign variable */
214  double ts; /* small t */
215  double theta; /* angle */
216 
217  x -= cache_ptr->false_easting;
218  y = cache_ptr->rh - y + cache_ptr->false_northing;
219  if (cache_ptr->ns > 0)
220  {
221  rh1 = sqrt (x * x + y * y);
222  con = 1.0;
223  }
224  else
225  {
226  rh1 = -sqrt (x * x + y * y);
227  con = -1.0;
228  }
229  theta = 0.0;
230  if (rh1 != 0)
231  theta = atan2((con * x),(con * y));
232 
233  if ((rh1 != 0) || (cache_ptr->ns > 0.0))
234  {
235  con = 1.0/cache_ptr->ns;
236  ts = pow((rh1/(cache_ptr->r_major * cache_ptr->f0)),con);
237  if (gctp_calc_phi2(cache_ptr->e, ts, lat) != GCTP_SUCCESS)
238  return GCTP_ERROR;
239  }
240  else
241  *lat = -HALF_PI;
242 
243  *lon = adjust_lon(theta/cache_ptr->ns + cache_ptr->center_lon);
244 
245  return GCTP_SUCCESS;
246 }
247 
248 /*****************************************************************************
249 Name: forward_transform
250 
251 Purpose: Transforms lat,long to X,Y
252 
253 Returns:
254  GCTP_SUCCESS or GCTP_ERROR
255 
256 *****************************************************************************/
257 static int forward_transform
258 (
259  const TRANSFORMATION *trans, /* I: transformation information */
260  double lon, /* I: Longitude */
261  double lat, /* I: Latitude */
262  double *x, /* O: X projection coordinate */
263  double *y /* O: Y projection coordinate */
264 )
265 {
266  struct lamcc_proj *cache_ptr = (struct lamcc_proj *)trans->cache;
267  double con; /* temporary angle variable */
268  double rh1; /* height above ellipsoid */
269  double sinphi; /* sin value */
270  double theta; /* angle */
271  double ts; /* small value t */
272 
273  con = fabs( fabs(lat) - HALF_PI);
274  if (con > EPSLN)
275  {
276  sinphi = sin(lat);
277  ts = gctp_calc_small_t(cache_ptr->e, lat, sinphi);
278  rh1 = cache_ptr->r_major * cache_ptr->f0 * pow(ts, cache_ptr->ns);
279  }
280  else
281  {
282  con = lat * cache_ptr->ns;
283  if (con <= 0)
284  {
285  GCTP_PRINT_ERROR("Point can not be projected");
286  return GCTP_ERROR;
287  }
288  rh1 = 0;
289  }
290  theta = cache_ptr->ns * adjust_lon(lon - cache_ptr->center_lon);
291  *x = rh1 * sin(theta) + cache_ptr->false_easting;
292  *y = cache_ptr->rh - rh1 * cos(theta) + cache_ptr->false_northing;
293 
294  return GCTP_SUCCESS;
295 }
296 
297 /*****************************************************************************
298 Name: gctp_lamcc_inverse_init
299 
300 Purpose: Initializes the inverse Lambert Conformal Conic transformation
301 
302 Returns:
303  GCTP_SUCCESS or GCTP_ERROR
304 
305 *****************************************************************************/
307 (
308  TRANSFORMATION *trans
309 )
310 {
311  /* Call the common routine used for the forward and inverse init */
312  if (common_init(trans) != GCTP_SUCCESS)
313  {
315  "Error initializing lambert conformal conic inverse projection");
316  return GCTP_ERROR;
317  }
318 
319  trans->transform = inverse_transform;
320 
321  return GCTP_SUCCESS;
322 }
323 
324 /*****************************************************************************
325 Name: gctp_lamcc_forward_init
326 
327 Purpose: Initializes the forward Lambert Conformal Conic transformation
328 
329 Returns:
330  GCTP_SUCCESS or GCTP_ERROR
331 
332 *****************************************************************************/
334 (
335  TRANSFORMATION *trans
336 )
337 {
338  /* Call the common routine used for the forward and inverse init */
339  if (common_init(trans) != GCTP_SUCCESS)
340  {
342  "Error initializing lambert conformal conic forward projection");
343  return GCTP_ERROR;
344  }
345 
346  trans->transform = forward_transform;
347 
348  return GCTP_SUCCESS;
349 }
#define SQUARE(x)
Definition: proj_define.h:99
void gctp_print_origin(double A)
Definition: gctp_report.c:65
int gctp_lamcc_forward_init(TRANSFORMATION *trans)
void gctp_print_title(const char *proj_name)
Definition: gctp_report.c:14
#define GCTP_ERROR
Definition: gctp.h:81
#define NULL
Definition: decode_rs.h:63
#define GCTP_PRINT_ERROR(format,...)
Definition: oli_local.h:81
float * lat
character(len=1000) if
Definition: names.f90:13
double adjust_lon(double x)
Definition: proj_cproj.c:349
int gctp_lamcc_inverse_init(TRANSFORMATION *trans)
def cache(filename, recache=False)
Definition: utils.py:145
void gctp_print_stanparl(double A, double B)
Definition: gctp_report.c:73
int gctp_calc_phi2(double eccent, double ts, double *phi2)
Definition: gctp_utility.c:209
#define HALF_PI
Definition: proj_define.h:84
double gctp_calc_small_t(double eccent, double phi, double sinphi)
Definition: gctp_utility.c:277
void gctp_print_offsetp(double A, double B)
Definition: gctp_report.c:91
#define GCTP_SUCCESS
Definition: gctp.h:82
int gctp_dms2degrees(double angle, double *degrees)
#define gctp_get_spheroid
Definition: oli_local.h:7
#define sincos
Definition: proj_define.h:108
void gctp_print_radius2(double radius1, double radius2)
Definition: gctp_report.c:30
void gctp_print_cenlonmer(double A)
Definition: gctp_report.c:48
double gctp_calc_small_radius(double eccent, double sinphi, double cosphi)
Definition: gctp_utility.c:253
#define fabs(a)
Definition: misc.h:93
#define S2R
Definition: proj_define.h:92
float * lon
void radius(double A)
Definition: proj_report.c:132
#define EPSLN
Definition: proj_define.h:86