source: trunk/WRF.COMMON/WRFV2/external/io_grib_share/gridnav.c @ 3333

Last change on this file since 3333 was 11, checked in by aslmd, 15 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 15.9 KB
Line 
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
11int 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
40int 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
85int 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
177int 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
256int 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
397int 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
451int 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}
Note: See TracBrowser for help on using the repository browser.