wrfv2_fire/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
  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                         fabs(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 }