gridnav.c
References to this file elsewhere.
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include "gridnav.h"
5
6 #define MISSING -999
7 #define PI 3.1415927
8 #define RAD_TO_DEG 57.29577951
9 #define EARTH_RAD 6371.200
10
11 int fill_proj_parms(GridNav *gridnav);
12
13 /*****************************************************************************
14 * Fills up the gridnav structure using arguments input to the function
15 *
16 * input:
17 * central_lat - the central latitude of the projection, in degrees
18 * central_lon - the central longitude of the projection, in degrees
19 * projection - the projection number (defined by #define's in gridnav.h)
20 * truelat1 - the first "true" latitude. Only has an effect with
21 * polar stereographic and lambert projections
22 * truelat2 - the second true latitude. Only has an effect with lambert
23 * projection
24 * num_columns - number of columns in grid
25 * num_rows - number of rows in grid
26 * dx,dy - east-west (dx) and north-south (dy) distance between grid
27 * points, in degrees for latlon (i.e., cylindrical
28 * equidistant) projection, in km for all other projections
29 * (including mercator).
30 * lat_origin - latitude of grid point defined in origin_column, origin_row
31 * lon_origin - longitude of grid point defined in origin_column, origin_row
32 * origin_column - column for lat_origin, long_origin pair
33 * origin_row - row for lat_origin, long_origin pair
34 *
35 * output:
36 * gridnav - filled gridnav structure
37 *
38 *****************************************************************************/
39
40 int GRID_init(float central_lat, float central_lon, int projection,
41 float truelat1, float truelat2, int num_columns,
42 int num_rows, float dx, float dy, float lat_origin,
43 float lon_origin, float origin_column, float origin_row,
44 GridNav *gridnav)
45 {
46
47 int status = 1;
48
49 gridnav->proj.central_lon = central_lon;
50 gridnav->proj.central_lat = central_lat;
51 gridnav->proj.map_proj = projection;
52 gridnav->proj.truelat1 = truelat1;
53 gridnav->proj.truelat2 = truelat2;
54
55 gridnav->grid.num_columns = num_columns;
56 gridnav->grid.num_rows = num_rows;
57 gridnav->grid.dx = dx;
58 gridnav->grid.dy = dy;
59 gridnav->grid.lat_origin = lat_origin;
60 gridnav->grid.lon_origin = lon_origin;
61 gridnav->grid.origin_column = origin_column;
62 gridnav->grid.origin_row = origin_row;
63
64 fill_proj_parms(gridnav);
65
66 return status;
67 }
68
69 /*****************************************************************************
70 * Outputs the latitude and longitude of a grid point.
71 *
72 * input:
73 * *gridnav - a pointer to a filled gridnav structure. The gridnav
74 * structure should have been filled using one of the init
75 * functions.
76 * column - the column in the grid to retrieve.
77 * row - the row in the grid to retrieve.
78 *
79 * output:
80 * *lat - pointer to output latitude
81 * *lon - pointer to output longitude
82 *
83 *****************************************************************************/
84
85 int GRID_to_latlon(GridNav *gridnav, float column, float row, float *lat,
86 float *lon)
87 {
88 /*
89 * Calculate the latitude and longitude for an input grid point location
90 */
91 double X, Y, R;
92 int status = 1;
93
94 switch (gridnav->proj.map_proj)
95 {
96 case GRID_MERCATOR:
97 *lat = 2 * RAD_TO_DEG *
98 atan(exp((gridnav->proj_transform.parm2 + gridnav->grid.dx *
99 (row - gridnav->grid.origin_row)) /
100 gridnav->proj_transform.parm1 )) - 90;
101 *lon = gridnav->grid.lon_origin + RAD_TO_DEG * gridnav->grid.dx *
102 (column - gridnav->grid.origin_column) /
103 gridnav->proj_transform.parm1;
104 break;
105
106 case GRID_LAMCON:
107 X = (column - gridnav->grid.origin_column) * gridnav->grid.dx +
108 gridnav->proj_transform.parm6;
109 Y = (row - gridnav->grid.origin_row) * gridnav->grid.dy +
110 gridnav->proj_transform.parm3 + gridnav->proj_transform.parm7;
111 R = sqrt(X*X + Y*Y);
112 *lat = gridnav->proj_transform.parm5 * 90 -
113 2 * RAD_TO_DEG * atan(gridnav->proj_transform.parm4 *
114 pow(R, 1 /
115 gridnav->proj_transform.parm2));
116 *lon = gridnav->proj.central_lon +
117 RAD_TO_DEG / gridnav->proj_transform.parm2 *
118 atan(X / (gridnav->proj_transform.parm5 * -Y));
119 while (*lon > 180) *lon -= 360;
120 while (*lon <= -180) *lon += 360;
121 break;
122
123 case GRID_POLSTR:
124 X = (column - gridnav->grid.origin_column)*gridnav->grid.dx;
125 Y = (row - gridnav->grid.origin_row)*gridnav->grid.dy +
126 gridnav->proj_transform.parm3;
127 R = sqrt(X*X + Y*Y);
128 *lat = gridnav->proj_transform.parm5*90 - 2 * RAD_TO_DEG *
129 atan((gridnav->proj_transform.parm5*R/EARTH_RAD)/
130 (1+cos(gridnav->proj_transform.parm1)));
131 *lon = gridnav->grid.lon_origin + RAD_TO_DEG *
132 atan(X/(gridnav->proj_transform.parm5 * -Y));
133 while (*lon > 180) *lon -= 360;
134 while (*lon <= -180) *lon += 360;
135 break;
136
137 case GRID_LATLON:
138 *lat = (row - gridnav->grid.origin_row)*gridnav->grid.dy +
139 gridnav->grid.lat_origin;
140 *lon = (column - gridnav->grid.origin_column)*gridnav->grid.dx +
141 gridnav->grid.lon_origin;
142 break;
143
144 default:
145 /*
146 * Unsupported map projection: set lat-lon to no-data values.
147 */
148 fprintf(stderr,"GRID_to_latlon: Unsupport map projection type %d\n",
149 gridnav->proj.map_proj);
150 *lon = -9999;
151 *lat = -9999;
152 status = 0;
153 break;
154
155 } /* end of switch on map projection type */
156
157
158 return status;
159 }
160
161
162 /*****************************************************************************
163 * Outputs grid point indices, given lat/lon location.
164 *
165 * input:
166 * *gridnav - a pointer to a filled gridnav structure. The gridnav
167 * structure should have been filled using one of the init
168 * functions.
169 * lat - latitude for location
170 * lon - longitude for location
171 * output:
172 * *column - pointer to value for column
173 * *row - pointer to value for row
174 *
175 *****************************************************************************/
176
177 int GRID_from_latlon(GridNav *gridnav, float lat, float lon, float *column,
178 float *row)
179 {
180 double X, Y, Rs;
181 double lat_rad;
182 int status = 1;
183
184 switch (gridnav->proj.map_proj)
185 {
186 case GRID_MERCATOR:
187 X = gridnav->proj_transform.parm1 *
188 ((lon - gridnav->grid.lon_origin) / RAD_TO_DEG);
189 *column = gridnav->grid.origin_column + X / gridnav->grid.dx;
190 lat_rad = lat / RAD_TO_DEG;
191 Y = gridnav->proj_transform.parm1 * log((1 + sin(lat_rad)) /
192 cos(lat_rad));
193 *row = gridnav->grid.origin_row +
194 ((Y-gridnav->proj_transform.parm2) / gridnav->grid.dy);;
195 break;
196
197 case GRID_LAMCON:
198 Rs = (EARTH_RAD / gridnav->proj_transform.parm2) *
199 sin(gridnav->proj_transform.parm1) *
200 pow(tan((gridnav->proj_transform.parm5 * 90 - lat) /
201 (2 * RAD_TO_DEG))/
202 tan(gridnav->proj_transform.parm1 / 2),
203 gridnav->proj_transform.parm2);
204 *row = gridnav->grid.origin_row -
205 (1 / gridnav->grid.dy) *
206 (gridnav->proj_transform.parm3 +
207 Rs * cos(gridnav->proj_transform.parm2*
208 (lon - gridnav->proj.central_lon) / RAD_TO_DEG)) -
209 gridnav->proj_transform.parm7 / gridnav->grid.dy;
210 *column = gridnav->grid.origin_column +
211 ( gridnav->proj_transform.parm5*
212 ((Rs / gridnav->grid.dx)*
213 sin(gridnav->proj_transform.parm2 *
214 (lon - gridnav->proj.central_lon) / RAD_TO_DEG) -
215 gridnav->proj_transform.parm6 / gridnav->grid.dx));
216 break;
217
218 case GRID_POLSTR:
219 Rs = EARTH_RAD * sin((gridnav->proj_transform.parm5 * 90 - lat)
220 / RAD_TO_DEG)*
221 ((1 + cos(gridnav->proj_transform.parm1)) /
222 (1 + cos((gridnav->proj_transform.parm5 * 90 - lat)
223 / RAD_TO_DEG)) );
224 *row = gridnav->grid.origin_row -
225 (1 / gridnav->grid.dy) *
226 (gridnav->proj_transform.parm3 +
227 Rs * cos((lon - gridnav->grid.lon_origin) / RAD_TO_DEG));
228 *column = gridnav->grid.origin_column +
229 gridnav->proj_transform.parm5 *
230 ((Rs / gridnav->grid.dx) *
231 sin((lon - gridnav->grid.lon_origin) / RAD_TO_DEG));
232 break;
233
234 case GRID_LATLON:
235 *row = ((lat - gridnav->grid.lat_origin) / gridnav->grid.dy) +
236 gridnav->grid.origin_row;
237 /* If lon is negative, make it positive */
238 while (lon < 0) lon += 360.;
239 *column = ((lon - gridnav->grid.lon_origin) / gridnav->grid.dx) +
240 gridnav->grid.origin_column;
241 break;
242
243 default:
244 fprintf(stderr,"GRID_from_latlon: Unsupported map projection type %d\n",
245 gridnav->proj.map_proj);
246 *column = -9999;
247 *row = -9999;
248 status = 0;
249 break;
250
251 } /* End of switch on map projection type */
252
253 return status;
254 }
255
256 int fill_proj_parms(GridNav *gridnav)
257 {
258 double orig_lat_rad;
259 double R_orig;
260 int hemifactor;
261
262 switch (gridnav->proj.map_proj)
263 {
264 case GRID_MERCATOR:
265 gridnav->proj_transform.parm1 =
266 EARTH_RAD * cos(gridnav->proj.truelat1 / RAD_TO_DEG);
267 orig_lat_rad = (gridnav->grid.lat_origin) / RAD_TO_DEG;
268 gridnav->proj_transform.parm2 =
269 gridnav->proj_transform.parm1 *
270 log((1. + sin(orig_lat_rad)) / cos(orig_lat_rad));
271 gridnav->proj_transform.parm3 = MISSING;
272 gridnav->proj_transform.parm4 = MISSING;
273 gridnav->proj_transform.parm5 = MISSING;
274 break;
275 case GRID_LAMCON:
276 if (gridnav->proj.truelat1 >= 0)
277 {
278 hemifactor = 1;
279 }
280 else
281 {
282 hemifactor = -1;
283 }
284 /* This is Psi1 in MM5 speak */
285 gridnav->proj_transform.parm1 =
286 hemifactor*(PI/2 - fabs(gridnav->proj.truelat1) / RAD_TO_DEG);
287 /* This is Kappa in MM5 speak */
288 if (fabs(gridnav->proj.truelat1 - gridnav->proj.truelat2) < .01)
289 {
290 gridnav->proj_transform.parm2 =
291 sin(gridnav->proj.truelat1 / RAD_TO_DEG);
292 }
293 else
294 {
295 gridnav->proj_transform.parm2 =
296 (log10(cos(gridnav->proj.truelat1 / RAD_TO_DEG)) -
297 log10(cos(gridnav->proj.truelat2 / RAD_TO_DEG))) /
298 (log10(tan((45 - fabs(gridnav->proj.truelat1) / 2) / RAD_TO_DEG)) -
299 log10(tan((45 - fabs(gridnav->proj.truelat2) / 2) / RAD_TO_DEG)));
300 }
301 /* This is Yc in MM5 speak */
302 gridnav->proj_transform.parm3 =
303 -EARTH_RAD/gridnav->proj_transform.parm2 *
304 sin(gridnav->proj_transform.parm1) *
305 pow(
306 (tan((hemifactor * 90 - gridnav->grid.lat_origin) /
307 (RAD_TO_DEG * 2)) /
308 tan(gridnav->proj_transform.parm1 / 2)),
309 gridnav->proj_transform.parm2);
310 gridnav->proj_transform.parm4 =
311 tan(gridnav->proj_transform.parm1 / 2)*
312 pow(hemifactor * (gridnav->proj_transform.parm2)/
313 (EARTH_RAD * sin(gridnav->proj_transform.parm1)),
314 1 / gridnav->proj_transform.parm2);
315 gridnav->proj_transform.parm5 = hemifactor;
316 R_orig = (EARTH_RAD / gridnav->proj_transform.parm2) *
317 sin(gridnav->proj_transform.parm1) *
318 pow(tan((gridnav->proj_transform.parm5 * 90 -
319 gridnav->grid.lat_origin) /
320 (2 * RAD_TO_DEG))/
321 tan(gridnav->proj_transform.parm1 / 2),
322 gridnav->proj_transform.parm2);
323 /* X origin */
324 gridnav->proj_transform.parm6 =
325 R_orig * sin(gridnav->proj_transform.parm2 *
326 (gridnav->grid.lon_origin -
327 gridnav->proj.central_lon) / RAD_TO_DEG);
328
329 /* Y origin */
330 gridnav->proj_transform.parm7 = -(gridnav->proj_transform.parm3 +
331 R_orig * cos(gridnav->proj_transform.parm2 *
332 (gridnav->grid.lon_origin -
333 gridnav->proj.central_lon) / RAD_TO_DEG));
334 break;
335 case GRID_POLSTR:
336 if (gridnav->proj.truelat1 > 0)
337 {
338 hemifactor = 1;
339 }
340 else
341 {
342 hemifactor = -1;
343 }
344
345 /* This is Psi1 in MM5 speak */
346 gridnav->proj_transform.parm1 =
347 hemifactor * (PI/2 - fabs(gridnav->proj.truelat1) / RAD_TO_DEG);
348 gridnav->proj_transform.parm2 =
349 (1+log10(cos(gridnav->proj.truelat1 / RAD_TO_DEG))) /
350 ( -log10(tan(45 / RAD_TO_DEG - hemifactor *
351 gridnav->proj.truelat1 /
352 (2 * RAD_TO_DEG) )) );
353 /* This is Yc in MM5 speak */
354 gridnav->proj_transform.parm3 =
355 -EARTH_RAD * sin((hemifactor * 90 -
356 gridnav->grid.lat_origin) / RAD_TO_DEG)*
357 ( (1 + cos(gridnav->proj_transform.parm1))/
358 (1 + cos((hemifactor*90 - gridnav->grid.lat_origin) /
359 RAD_TO_DEG)) );
360 gridnav->proj_transform.parm4 = MISSING;
361 gridnav->proj_transform.parm5 = hemifactor;
362 break;
363 case GRID_LATLON:
364 gridnav->proj_transform.parm1 = MISSING;
365 gridnav->proj_transform.parm2 = MISSING;
366 gridnav->proj_transform.parm3 = MISSING;
367 gridnav->proj_transform.parm4 = MISSING;
368 gridnav->proj_transform.parm5 = MISSING;
369 break;
370
371 default:
372 fprintf(stderr,"GRID_init_mm5data: Invalid Projection\n");
373 return 0;
374 }
375 return 1;
376 }
377
378 /*****************************************************************************
379 * Rotates u and v components of wind, from being relative to earth, to
380 * relative to grid.
381 *
382 * input:
383 * *gridnav - a pointer to a filled gridnav structure. The gridnav
384 * structure should have been filled using one of the init
385 * functions.
386 * lon - longitude for location
387 * u_earth - the u component of the wind in earth (north-relative)
388 * coordinates.
389 * v_earth - the v component of the wind in earth (north-relative)
390 * coordinates.
391 * output:
392 * *u_grid - pointer to value for u in grid coordinates
393 * *v_grid - pointer to value for v in grid coordinates
394 *
395 *****************************************************************************/
396
397 int GRID_rotate_from_earth_coords(GridNav *gridnav, float lon, float u_earth,
398 float v_earth, float *u_grid, float *v_grid)
399 {
400 float speed, dir;
401 float dir_grid;
402
403 /* Calculate Speed and Direction from u,v */
404 switch (gridnav->proj.map_proj)
405 {
406 case GRID_MERCATOR:
407 *u_grid = u_earth;
408 *v_grid = v_earth;
409 break;
410 case GRID_POLSTR: case GRID_LAMCON:
411 speed = sqrt(u_earth * u_earth + v_earth * v_earth);
412 dir = RAD_TO_DEG * atan2(-u_earth,-v_earth);
413 while (dir >= 360) dir -= 360;
414 while (dir < 0) dir += 360;
415
416 dir_grid = dir - (lon - gridnav->proj.central_lon) *
417 gridnav->proj_transform.parm2;
418 while (dir_grid >= 360) dir_grid -= 360;
419 while (dir_grid < 0) dir_grid += 360;
420
421 *u_grid = -1. * speed * sin(dir_grid / RAD_TO_DEG);
422 *v_grid = -1. * speed * cos(dir_grid / RAD_TO_DEG);
423 break;
424 default:
425 fprintf(stderr,
426 "GRID_rotate_from_earth_coords: Invalid Projection\n");
427 return 0;
428 } /* End of switch projection */
429
430 return 1;
431 }
432
433 /*****************************************************************************
434 * Rotates u and v components of wind, from being relative to earth, to
435 * relative to grid.
436 *
437 * input:
438 * *gridnav - a pointer to a filled gridnav structure. The gridnav
439 * structure should have been filled using one of the init
440 * functions.
441 * lon - longitude for location
442 * u_grid - the u component of the wind in grid coordinates
443 * v_grid - the v component of the wind in grid coordinates
444 *
445 * output:
446 * *u_earth - pointer to value for u in earth-relative coordinates
447 * *v_grid - pointer to value for v in earth-relative coordinates
448 *
449 *****************************************************************************/
450
451 int GRID_rotate_to_earth_coords(GridNav *gridnav, float lon, float u_grid,
452 float v_grid, float *u_earth, float *v_earth)
453 {
454 float speed, dir_grid;
455 float dir_earth;
456
457 /* Calculate Speed and Direction from u,v */
458 switch (gridnav->proj.map_proj)
459 {
460 case GRID_MERCATOR:
461 *u_earth = u_grid;
462 *v_earth = v_grid;
463 break;
464 case GRID_POLSTR: case GRID_LAMCON:
465 speed = sqrt(u_grid * u_grid + v_grid * v_grid);
466 dir_grid = RAD_TO_DEG * atan2(-u_grid, -v_grid);
467 while (dir_grid >= 360) dir_grid -= 360;
468 while (dir_grid < 0) dir_grid += 360;
469
470 dir_earth = dir_grid + (lon - gridnav->proj.central_lon) *
471 gridnav->proj_transform.parm2;
472
473 while (dir_earth >= 360) dir_earth -= 360;
474 while (dir_earth < 0) dir_earth += 360;
475
476 *u_earth = -1. * speed * sin(dir_earth / RAD_TO_DEG);
477 *v_earth = -1. * speed * cos(dir_earth / RAD_TO_DEG);
478 break;
479 default:
480 fprintf(stderr,
481 "GRID_rotate_to_earth_coords: Invalid Projection\n");
482 return 0;
483 } /* End of switch projection */
484 return 1;
485 }