OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
gctp_create_transformation.c
Go to the documentation of this file.
1 /****************************************************************************
2 Name: gctp_create_transformation
3 
4 Purpose: Given input and output projections creates a projection
5  transformation object (GTCP_TRANSFORMATION) that is used when calling
6  the gctp_transform to do the actual transformation.
7 
8 Note: The GCTP_TRANSFORMATION includes the input and output units specified.
9  If the same transformation is needed with different input or output
10  units it requires the creation of a different GCTP_TRANSFORMATION.
11 
12 ****************************************************************************/
13 #include <stdlib.h>
14 #include <stdio.h>
15 #include "oli_cproj.h"
16 #include "gctp.h"
17 #include "oli_local.h"
18 
19 /* Define a type for a projection initialization function pointer. */
20 typedef int (*INIT_ROUTINE)(TRANSFORMATION *trans);
21 
22 /* Define a lookup table for the forward transform init routines
23  (lat/lon to x/y) */
24 static INIT_ROUTINE forward_init[GCTP_MAX_PROJ_CODE + 1] =
25  {
26  gctp_geo_init, /* 0 = Geographic */
27  gctp_utm_forward_init, /* 1 = Universal Transverse Mercator (UTM) */
28  gctp_state_plane_forward_init, /* 2 = State Plane Coordinates */
29  NULL, /* 3 = Albers Conical Equal Area */
30  gctp_lamcc_forward_init, /* 4 = Lambert Conformal Conic */
31  NULL, /* 5 = Mercator */
32  gctp_ps_forward_init, /* 6 = Polar Stereographic */
33  gctp_poly_forward_init,/* 7 = Polyconic */
34  NULL, /* 8 = Equidistant Conic */
35  gctp_tm_forward_init, /* 9 = Transverse Mercator */
36  NULL, /* 10 = Stereographic */
37  NULL, /* 11 = Lambert Azimuthal Equal Area */
38  NULL, /* 12 = Azimuthal Equidistant */
39  NULL, /* 13 = Gnomonic */
40  NULL, /* 14 = Orthographic */
41  NULL, /* 15 = General Vertical Near-Side Perspective */
42  NULL, /* 16 = Sinusiodal */
43  NULL, /* 17 = Equirectangular */
44  NULL, /* 18 = Miller Cylindrical */
45  NULL, /* 19 = Van der Grinten */
46  gctp_om_forward_init, /* 20 = (Hotine) Oblique Mercator */
47  NULL, /* 21 = Robinson */
48  gctp_som_forward_init, /* 22 = Space Oblique Mercator (SOM) */
49  NULL, /* 23 = Alaska Conformal */
50  NULL, /* 24 = Interrupted Goode Homolosine */
51  NULL, /* 25 = Mollweide */
52  NULL, /* 26 = Interrupted Mollweide */
53  NULL, /* 27 = Hammer */
54  NULL, /* 28 = Wagner IV */
55  NULL, /* 29 = Wagner VII */
56  NULL, /* 30 = Oblated Equal Area */
57  NULL, /* 31 = Integerized Sinusiodal */
58  };
59 
60 /* Define a lookup table for the inverse transform init routines
61  (x/y to lat/lon) */
62 static INIT_ROUTINE inverse_init[GCTP_MAX_PROJ_CODE + 1] =
63  {
64  gctp_geo_init, /* 0 = Geographic */
65  gctp_utm_inverse_init, /* 1 = Universal Transverse Mercator (UTM) */
66  gctp_state_plane_inverse_init, /* 2 = State Plane Coordinates */
67  NULL, /* 3 = Albers Conical Equal Area */
68  gctp_lamcc_inverse_init, /* 4 = Lambert Conformal Conic */
69  NULL, /* 5 = Mercator */
70  gctp_ps_inverse_init, /* 6 = Polar Stereographic */
71  gctp_poly_inverse_init,/* 7 = Polyconic */
72  NULL, /* 8 = Equidistant Conic */
73  gctp_tm_inverse_init, /* 9 = Transverse Mercator */
74  NULL, /* 10 = Stereographic */
75  NULL, /* 11 = Lambert Azimuthal Equal Area */
76  NULL, /* 12 = Azimuthal Equidistant */
77  NULL, /* 13 = Gnomonic */
78  NULL, /* 14 = Orthographic */
79  NULL, /* 15 = General Vertical Near-Side Perspective */
80  NULL, /* 16 = Sinusiodal */
81  NULL, /* 17 = Equirectangular */
82  NULL, /* 18 = Miller Cylindrical */
83  NULL, /* 19 = Van der Grinten */
84  gctp_om_inverse_init, /* 20 = (Hotine) Oblique Mercator */
85  NULL, /* 21 = Robinson */
86  gctp_som_inverse_init, /* 22 = Space Oblique Mercator (SOM) */
87  NULL, /* 23 = Alaska Conformal */
88  NULL, /* 24 = Interrupted Goode Homolosine */
89  NULL, /* 25 = Mollweide */
90  NULL, /* 26 = Interrupted Mollweide */
91  NULL, /* 27 = Hammer */
92  NULL, /* 28 = Wagner IV */
93  NULL, /* 29 = Wagner VII */
94  NULL, /* 30 = Oblated Equal Area */
95  NULL, /* 31 = Integerized Sinusiodal */
96  };
97 
98 /* Flag to only allow threadsafe transforms to be created. This is "temporary"
99  and can be removed when all the projection transformations have been updated
100  to be threadsafe. */
101 static int only_threadsafe = 0;
102 
103 /* Routine to get the conversion factor for converting between the input and
104  output units. */
105 static int get_unit_conversion_factor
106 (
107  int in_units,
108  int out_units,
109  double *factor
110 )
111 {
112  /* Define an array of conversion factors between the in_units (first index
113  in the array) and the out_units (second index of the array) */
114  static const double factors[6][6] =
115  {{1.0, 0.0, 0.0, 206264.8062470963, 57.295779513082323, 0.0},
116  {0.0, 1.0, .3048006096012192, 0.0, 0.0, 1.000002000004},
117  {0.0, 3.280833333333333, 1.0, 0.0, 0.0, 3.280839895013124},
118  {.484813681109536e-5, 0.0, 0.0, 1.0, .27777777777778e-3, 0.0},
119  {.01745329251994329, 0.0, 0.0, 3600, 1.0, 0.0},
120  {0.0, .999998, .3048, 0.0, 0.0, 1.0}};
121 
122  if ((out_units >= 0) && (out_units <= MAXUNIT) && (in_units >= 0)
123  && (in_units <= MAXUNIT))
124  {
125  *factor = factors[in_units][out_units];
126 
127  /* Angle units can not be converted to length units */
128  if (*factor == 0.0)
129  {
130  GCTP_PRINT_ERROR("Incompatable unit codes: %d and %d",
131  in_units, out_units);
132  return GCTP_ERROR;
133  }
134  }
135  else
136  {
137  GCTP_PRINT_ERROR("Illegal source or target unit code: %d and %d",
138  in_units, out_units);
139  return GCTP_ERROR;
140  }
141 
142  return GCTP_SUCCESS;
143 }
144 
145 /****************************************************************************
146 Name: gctp_create_transformation
147 
148 Purpose: Given input and output projections creates a projection
149  transformation.
150 
151 Returns:
152  A pointer to the transformation created, or NULL if there is an error.
153 
154 Notes:
155  - gctp_destroy_transformation should be called for transformations that
156  are created to free any resources allocated by the creation process.
157 
158 ****************************************************************************/
159 GCTP_TRANSFORMATION *gctp_create_transformation
160 (
161  const GCTP_PROJECTION *input_projection,
162  const GCTP_PROJECTION *output_projection
163 )
164 {
165  GCTP_TRANSFORMATION *trans;
166  INIT_ROUTINE inverse_init_func;
167  INIT_ROUTINE forward_init_func;
168  int units;
169 
170  /* Verify the proj codes are valid */
171  if (input_projection->proj_code < 0
172  || input_projection->proj_code > GCTP_MAX_PROJ_CODE)
173  {
174  GCTP_PRINT_ERROR("Invalid input projection code: %d",
175  input_projection->proj_code);
176  return NULL;
177  }
178 
179  if (output_projection->proj_code < 0
180  || output_projection->proj_code > GCTP_MAX_PROJ_CODE)
181  {
182  GCTP_PRINT_ERROR("Invalid output projection code: %d",
183  input_projection->proj_code);
184  return NULL;
185  }
186 
187  /* Allocate memory for the transformation object */
188  trans = malloc(sizeof(*trans));
189  if (!trans)
190  {
191  GCTP_PRINT_ERROR("Error allocating memory for transformation object");
192  return NULL;
193  }
194 
195  /* Copy the projection information to the transformation object */
196  trans->inverse.proj = *input_projection;
197  trans->forward.proj = *output_projection;
198 
199  /* default to no forward or inverse transformation routines */
200  trans->forward.transform = NULL;
201  trans->forward.destroy = NULL;
202  trans->forward.cache = NULL;
203  trans->forward.print_info = NULL;
204 
205  trans->inverse.transform = NULL;
206  trans->inverse.destroy = NULL;
207  trans->inverse.cache = NULL;
208  trans->inverse.print_info = NULL;
209 
210  /* Get the inverse and forward init routine pointers */
211  inverse_init_func = inverse_init[input_projection->proj_code];
212  forward_init_func = forward_init[output_projection->proj_code];
213 
214  /* If either of the init routines are not available, fall back to using
215  GCTP for the transformations */
216  /* TODO - remove this when all transformations are supported with the
217  new interface */
218  trans->use_gctp = 0;
219  if (!inverse_init_func || !forward_init_func)
220  {
221  /* If only threadsafe transforms are allowed, consider trying to use
222  the old gctp interface an error since it isn't threadsafe */
223  if (only_threadsafe)
224  {
225  GCTP_PRINT_ERROR("Error: trying to use a projection "
226  "transformation that isn't threadsafe");
227  free(trans);
228  return NULL;
229  }
230  trans->use_gctp = 1;
231  return trans;
232  }
233 
234  /* Create the inverse transformation */
235  if (inverse_init_func)
236  {
237  if (inverse_init_func(&trans->inverse) != GCTP_SUCCESS)
238  {
239  GCTP_PRINT_ERROR("Error initializing inverse transformation");
240  free(trans);
241  return NULL;
242  }
243  }
244  else
245  {
246  GCTP_PRINT_ERROR("Unsupported projection code: %d",
247  input_projection->proj_code);
248  free(trans);
249  return NULL;
250  }
251 
252  /* Create the forward transformation */
253  if (forward_init_func)
254  {
255  if (forward_init_func(&trans->forward) != GCTP_SUCCESS)
256  {
257  GCTP_PRINT_ERROR("Error initializing forward transformation");
259  return NULL;
260  }
261  }
262  else
263  {
264  GCTP_PRINT_ERROR("Unsupported projection code: %d",
265  output_projection->proj_code);
267  return NULL;
268  }
269 
270  /* find the factor unit conversions--all transformations are in radians
271  or meters */
272  if (input_projection->proj_code == GEO)
273  units = RADIAN;
274  else
275  units = METER;
276  if (get_unit_conversion_factor(input_projection->units, units,
277  &trans->inverse.unit_conversion_factor) != GCTP_SUCCESS)
278  {
279  GCTP_PRINT_ERROR("Error getting the input unit conversion factor");
281  return NULL;
282  }
283 
284  if (output_projection->proj_code == GEO)
285  units = RADIAN;
286  else
287  units = METER;
288  if (get_unit_conversion_factor(units, output_projection->units,
289  &trans->forward.unit_conversion_factor) != GCTP_SUCCESS)
290  {
291  GCTP_PRINT_ERROR("Error getting the output unit conversion factor");
293  return NULL;
294  }
295 
296  /* Return the transformation */
297  return trans;
298 }
299 
300 
301 /****************************************************************************
302 Name: gctp_destroy_transformation
303 
304 Purpose: Releases the resources allocated when a transformation is created.
305  After this is called, the transformation is no longer valid.
306 
307 Returns:
308  nothing
309 
310 ****************************************************************************/
312 (
313  GCTP_TRANSFORMATION *trans
314 )
315 {
316  if (trans)
317  {
318  if (trans->forward.destroy)
319  trans->forward.destroy(&trans->forward);
320  if (trans->inverse.destroy)
321  trans->inverse.destroy(&trans->inverse);
322  free(trans->forward.cache);
323  trans->forward.cache = NULL;
324  free(trans->inverse.cache);
325  trans->inverse.cache = NULL;
326  free(trans);
327  }
328 }
329 
330 /****************************************************************************
331 Name: gctp_only_allow_threadsafe_transforms
332 
333 Purpose: Many of the transformations haven't been made threadsafe yet. Until
334  they are, there needs to be a way to allow threaded applications to
335  indicate that they are threaded and therefore only threadsafe
336  transformations should be allowed. This routine should be called from any
337  application that uses threading. After calling this, non-threadsafe
338  projections will cause the transformation creation to return an error.
339 
340 Returns:
341  nothing
342 
343 ****************************************************************************/
345 {
346  only_threadsafe = 1;
347 }
int gctp_lamcc_forward_init(TRANSFORMATION *trans)
void gctp_destroy_transformation(GCTP_TRANSFORMATION *trans)
#define GCTP_ERROR
Definition: gctp.h:81
#define NULL
Definition: decode_rs.h:63
#define GCTP_PRINT_ERROR(format,...)
Definition: oli_local.h:81
int gctp_geo_init(TRANSFORMATION *trans)
Definition: geographic.c:38
int gctp_ps_inverse_init(TRANSFORMATION *trans)
#define GEO
Definition: proj_define.h:36
int gctp_lamcc_inverse_init(TRANSFORMATION *trans)
void gctp_only_allow_threadsafe_transforms()
int gctp_utm_inverse_init(TRANSFORMATION *trans)
Definition: tm.c:487
int(* INIT_ROUTINE)(TRANSFORMATION *trans)
int gctp_om_forward_init(TRANSFORMATION *trans)
int gctp_state_plane_forward_init(TRANSFORMATION *trans)
Definition: state_plane.c:429
#define MAXUNIT
Definition: proj_define.h:75
int gctp_tm_inverse_init(TRANSFORMATION *trans)
Definition: tm.c:539
#define GCTP_SUCCESS
Definition: gctp.h:82
#define RADIAN
Definition: gctp.h:89
int gctp_om_inverse_init(TRANSFORMATION *trans)
int gctp_tm_forward_init(TRANSFORMATION *trans)
Definition: tm.c:565
GCTP_TRANSFORMATION * gctp_create_transformation(const GCTP_PROJECTION *input_projection, const GCTP_PROJECTION *output_projection)
int gctp_poly_inverse_init(TRANSFORMATION *trans)
Definition: polyconic.c:307
int gctp_som_inverse_init(TRANSFORMATION *trans)
Definition: som.c:573
int gctp_som_forward_init(TRANSFORMATION *trans)
Definition: som.c:599
int gctp_ps_forward_init(TRANSFORMATION *trans)
int gctp_poly_forward_init(TRANSFORMATION *trans)
Definition: polyconic.c:334
int gctp_state_plane_inverse_init(TRANSFORMATION *trans)
Definition: state_plane.c:388
#define METER
Definition: gctp.h:91
#define GCTP_MAX_PROJ_CODE
Definition: oli_local.h:10
int gctp_utm_forward_init(TRANSFORMATION *trans)
Definition: tm.c:513