external/io_grib_share/get_region_center.c

   1 #include "get_region_center.h"
   2 #include "gridnav.h"
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5
   6 #include "wrf_projection.h"
   7
   8 int get_gridnav_projection(int wrf_projection);
   9
  10 /****************************************************************************
  11  *
  12  * This function calculates the center lat and lon of a region within a larger
  13  *   domain.  It is useful for calculating the center of the boundary regions
  14  *   in a domain.
  15  *
  16  ****************************************************************************/
  17
  18
  19 int GET_REGION_CENTER(char *MemoryOrderIn, int *projection,
  20                       float *domain_center_lat,
  21                       float *domain_center_lon, int *full_xsize,
  22                       int *full_ysize, float *dx, float *dy,
  23                       float *proj_central_lon,
  24                       int *proj_center_flag, float *truelat1,
  25                       float *truelat2, int *region_xsize, int *region_ysize,
  26                       float *region_center_lat, float *region_center_lon,
  27                       int strlen1)
  28 {
  29
  30   char *MemoryOrder;
  31   int grid_projection;
  32   float full_xcenter, full_ycenter;
  33   float region_xcenter, region_ycenter;
  34   int status;
  35   int orig;
  36   int x_pos, y_pos;
  37   GridNav gridnav;
  38
  39   MemoryOrder = (char *)malloc((strlen1+1)*sizeof(char));
  40   memcpy(MemoryOrder,MemoryOrderIn,strlen1);
  41   MemoryOrder[strlen1] = '\0';
  42
  43   grid_projection = get_gridnav_projection(*projection);
  44
  45   full_xcenter = (*full_xsize - 1) / 2.;
  46   full_ycenter = (*full_ysize - 1) / 2.;
  47   region_xcenter = (*region_xsize - 1) / 2.;
  48   region_ycenter = (*region_ysize - 1) / 2.;
  49
  50   orig = 0;
  51
  52   if (strncmp(MemoryOrder,"XS", 2) == 0)
  53     {
  54       x_pos = region_xcenter;
  55       y_pos = full_ycenter;
  56     }
  57   else if (strncmp(MemoryOrder,"XE", 2) == 0)
  58     {
  59       x_pos = (*full_xsize - 1) - region_xcenter;
  60       y_pos = full_ycenter;
  61     }
  62   else if (strncmp(MemoryOrder,"YS", 2) == 0)
  63     {
  64       x_pos = full_xcenter;
  65       y_pos = region_ycenter;
  66     }
  67   else if (strncmp(MemoryOrder,"YE", 2) == 0)
  68     {
  69       x_pos = full_xcenter;
  70       y_pos = (*full_ysize - 1) - region_ycenter;
  71     }
  72   else
  73     {
  74       orig = 1;
  75     }
  76
  77   if (orig == 1)
  78     {
  79       *region_center_lat = *domain_center_lat;
  80       *region_center_lon = *domain_center_lon;
  81       status = 0;
  82     }
  83   else
  84     {
  85       /* Initialize grid structure */
  86       /*
  87       status = GRID_init(grid_info->center_lat, grid_info->central_lon,
  88                          grid_projection,
  89                          grid_info->latin1, grid_info->latin2,
  90                          grid_info->xpoints, grid_info->ypoints,
  91                          grid_info->Di, grid_info->Dj,
  92                          grid_info->center_lat, grid_info->center_lon,
  93                          x_center, y_center,
  94                          &gridnav);
  95       */
  96       status = GRID_init(*domain_center_lat, *proj_central_lon,
  97                          grid_projection,
  98                          *truelat1, *truelat2,
  99                          *full_xsize, *full_ysize, *dx, *dy,
 100                          *domain_center_lat, *domain_center_lon,
 101                          full_xcenter, full_ycenter,
 102                          &gridnav);
 103       if (!status)
 104         {
 105           fprintf(stderr,"get_region_center: error from GRID_init\n");
 106         }
 107
 108       /* get lat/lon of center of region */
 109       status = GRID_to_latlon(&gridnav, x_pos, y_pos, region_center_lat,
 110                               region_center_lon);
 111       if (!status)
 112         {
 113           fprintf(stderr,
 114                   "get_region_cneter: error from GRID_to_latlon for first lat/lon\n");
 115         }
 116
 117     }
 118
 119   free(MemoryOrder);
 120   return status;
 121
 122 }
 123 /******************************************************************************
 124  * translates the grid projection identifier from the WRF id to the grib id.
 125  *****************************************************************************/
 126
 127 int get_gridnav_projection(int wrf_projection)
 128 {
 129   int gridnav_projection;
 130
 131   /* Set the grid projection in the gridnav units */
 132   switch (wrf_projection)
 133     {
 134     case WRF_LATLON:
 135     case WRF_CASSINI:
 136       gridnav_projection = GRID_LATLON;
 137       break;
 138     case WRF_MERCATOR:
 139       gridnav_projection = GRID_MERCATOR;
 140       break;
 141     case WRF_LAMBERT:
 142       gridnav_projection = GRID_LAMCON;
 143       break;
 144     case WRF_POLAR_STEREO:
 145       gridnav_projection = GRID_POLSTR;
 146       break;
 147     default:
 148       fprintf(stderr,"Error, invalid projection: %d\n",wrf_projection);
 149       gridnav_projection = -1;
 150     }
 151
 152   return gridnav_projection;
 153 }
 154
 155 int GET_LL_LATLON(float *central_lat, float *central_lon, int *projection,
 156                   float *latin1, float *latin2, int *nx, int *ny,
 157                   float *dx, float *dy, float *center_lat, float *center_lon,
 158                   float *LLLa, float *LLLo, float *URLa, float *URLo, int *ierr)
 159 {
 160
 161   int grid_projection;
 162   float x_center;
 163   float y_center;
 164   GridNav gridnav;
 165   int status;
 166
 167   grid_projection = get_gridnav_projection(*projection);
 168
 169   /* Get coords of center of grid */
 170   x_center = (*nx + 1)/2.;
 171   y_center = (*ny + 1)/2.;
 172
 173  /* Initialize grid structure */
 174   status = GRID_init(*central_lat, *central_lon, grid_projection,
 175                      *latin1, *latin2, *nx, *ny, *dx, *dy,
 176                      *center_lat, *center_lon, x_center, y_center,
 177                      &gridnav);
 178   if (!status)
 179     {
 180       fprintf(stderr,"write_grib: error from GRID_init\n");
 181       *ierr = 1;
 182       return(0);
 183     }
 184
 185   /* get lat/lon of lower left corner */
 186   status = GRID_to_latlon(&gridnav, 1, 1, LLLa, LLLo);
 187   if (!status)
 188     {
 189       fprintf(stderr,
 190               "write_grib: error from GRID_to_latlon for first lat/lon\n");
 191       *ierr = 1;
 192       return(0);
 193     }
 194
 195   /* get lat/lon of upper right corner */
 196   status = GRID_to_latlon(&gridnav, *nx, *ny, URLa, URLo);
 197   if (!status)
 198     {
 199       fprintf(stderr,
 200               "write_grib: error from GRID_to_latlon for first lat/lon\n");
 201       *ierr = 1;
 202       return(0);
 203     }
 204
 205   *ierr = 0;
 206   return(0);
 207
 208 }