get_region_center.c

References to this file elsewhere.
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       gridnav_projection = GRID_LATLON;
136       break;
137     case WRF_MERCATOR:
138       gridnav_projection = GRID_MERCATOR;
139       break;
140     case WRF_LAMBERT:
141       gridnav_projection = GRID_LAMCON;
142       break;
143     case WRF_POLAR_STEREO:
144       gridnav_projection = GRID_POLSTR;
145       break;
146     default:
147       fprintf(stderr,"Error, invalid projection: %d\n",wrf_projection);
148       gridnav_projection = -1;
149     }
150   
151   return gridnav_projection;
152 }
153 
154 int GET_LL_LATLON(float *central_lat, float *central_lon, int *projection,
155 		  float *latin1, float *latin2, int *nx, int *ny, 
156 		  float *dx, float *dy, float *center_lat, float *center_lon,
157 		  float *LLLa, float *LLLo, float *URLa, float *URLo, int *ierr)
158 {
159 
160   int grid_projection;
161   float x_center;
162   float y_center;
163   GridNav gridnav;
164   int status;
165 
166   grid_projection = get_gridnav_projection(*projection);
167 
168   /* Get coords of center of grid */
169   x_center = (*nx + 1)/2.;
170   y_center = (*ny + 1)/2.;
171 
172  /* Initialize grid structure */
173   status = GRID_init(*central_lat, *central_lon, grid_projection,
174 		     *latin1, *latin2, *nx, *ny, *dx, *dy,
175 		     *center_lat, *center_lon, x_center, y_center,
176 		     &gridnav);
177   if (!status)
178     {
179       fprintf(stderr,"write_grib: error from GRID_init\n");
180       *ierr = 1;
181       return;
182     }
183 
184   /* get lat/lon of lower left corner */
185   status = GRID_to_latlon(&gridnav, 1, 1, LLLa, LLLo);
186   if (!status)
187     {
188       fprintf(stderr,
189 	      "write_grib: error from GRID_to_latlon for first lat/lon\n");
190       *ierr = 1;
191       return;
192     }
193 
194   /* get lat/lon of upper right corner */
195   status = GRID_to_latlon(&gridnav, *nx, *ny, URLa, URLo);
196   if (!status)
197     {
198       fprintf(stderr,
199 	      "write_grib: error from GRID_to_latlon for first lat/lon\n");
200       *ierr = 1;
201       return;
202     }
203 
204   *ierr = 0;
205   return;
206 
207 }