external/io_grib_share/gridnav.c

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