1 /*
2 * (C) Copyright 2005- ECMWF.
3 *
4 * This software is licensed under the terms of the Apache Licence Version 2.0
5 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6 *
7 * In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
8 * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
9 */
10
11 #include "grib_api_internal.h"
12 #include <math.h>
13
14 /*
15 This is used by make_class.pl
16
17 START_CLASS_DEF
18 CLASS = iterator
19 SUPER = grib_iterator_class_gen
20 IMPLEMENTS = destroy
21 IMPLEMENTS = init;next
22 MEMBERS = double *lats
23 MEMBERS = double *lons
24 MEMBERS = long Nj
25 END_CLASS_DEF
26 */
27
28 /* START_CLASS_IMP */
29
30 /*
31
32 Don't edit anything between START_CLASS_IMP and END_CLASS_IMP
33 Instead edit values between START_CLASS_DEF and END_CLASS_DEF
34 or edit "iterator.class" and rerun ./make_class.pl
35
36 */
37
38
39 static void init_class(grib_iterator_class*);
40
41 static int init(grib_iterator* i, grib_handle*, grib_arguments*);
42 static int next(grib_iterator* i, double* lat, double* lon, double* val);
43 static int destroy(grib_iterator* i);
44
45
46 typedef struct grib_iterator_polar_stereographic
47 {
48 grib_iterator it;
49 /* Members defined in gen */
50 long carg;
51 const char* missingValue;
52 /* Members defined in polar_stereographic */
53 double* lats;
54 double* lons;
55 long Nj;
56 } grib_iterator_polar_stereographic;
57
58 extern grib_iterator_class* grib_iterator_class_gen;
59
60 static grib_iterator_class _grib_iterator_class_polar_stereographic = {
61 &grib_iterator_class_gen, /* super */
62 "polar_stereographic", /* name */
63 sizeof(grib_iterator_polar_stereographic), /* size of instance */
64 0, /* inited */
65 &init_class, /* init_class */
66 &init, /* constructor */
67 &destroy, /* destructor */
68 &next, /* Next Value */
69 0, /* Previous Value */
70 0, /* Reset the counter */
71 0, /* has next values */
72 };
73
74 grib_iterator_class* grib_iterator_class_polar_stereographic = &_grib_iterator_class_polar_stereographic;
75
76
init_class(grib_iterator_class * c)77 static void init_class(grib_iterator_class* c)
78 {
79 c->previous = (*(c->super))->previous;
80 c->reset = (*(c->super))->reset;
81 c->has_next = (*(c->super))->has_next;
82 }
83 /* END_CLASS_IMP */
84
next(grib_iterator * i,double * lat,double * lon,double * val)85 static int next(grib_iterator* i, double* lat, double* lon, double* val)
86 {
87 grib_iterator_polar_stereographic* self = (grib_iterator_polar_stereographic*)i;
88
89 if ((long)i->e >= (long)(i->nv - 1))
90 return 0;
91 i->e++;
92
93 *lat = self->lats[i->e];
94 *lon = self->lons[i->e];
95 *val = i->data[i->e];
96
97 return 1;
98 }
99
100 /* Data struct for Forward and Inverse Projections */
101 typedef struct proj_data_t
102 {
103 double centre_lon; /* central longitude */
104 double centre_lat; /* central latitude */
105 double sign; /* sign variable */
106 double ind; /* flag variable */
107 double mcs; /* small m */
108 double tcs; /* small t */
109 double false_northing; /* y offset in meters */
110 double false_easting; /* x offset in meters */
111 } proj_data_t;
112
113 #define RAD2DEG 57.29577951308232087684 /* 180 over pi */
114 #define DEG2RAD 0.01745329251994329576 /* pi over 180 */
115 #define PI_OVER_2 1.5707963267948966 /* half pi */
116 #define EPSILON 1.0e-10
117
init(grib_iterator * iter,grib_handle * h,grib_arguments * args)118 static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
119 {
120 int ret = 0;
121 double *lats, *lons; /* arrays for latitudes and longitudes */
122 double lonFirstInDegrees, latFirstInDegrees, radius;
123 double x, y, Dx, Dy;
124 long nx, ny, centralLongitudeInDegrees, centralLatitudeInDegrees;
125 long alternativeRowScanning, iScansNegatively, i, j;
126 long jScansPositively, jPointsAreConsecutive, southPoleOnPlane;
127 double centralLongitude, centralLatitude; /* in radians */
128 double con1; /* temporary angle */
129 double ts; /* value of small t */
130 double height; /* height above ellipsoid */
131 double x0, y0, lonFirst, latFirst;
132 proj_data_t fwd_proj_data = {0,};
133 proj_data_t inv_proj_data = {0,};
134
135 grib_iterator_polar_stereographic* self = (grib_iterator_polar_stereographic*)iter;
136
137 const char* sradius = grib_arguments_get_name(h, args, self->carg++);
138 const char* snx = grib_arguments_get_name(h, args, self->carg++);
139 const char* sny = grib_arguments_get_name(h, args, self->carg++);
140 const char* slatFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
141 const char* slonFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
142 const char* ssouthPoleOnPlane = grib_arguments_get_name(h, args, self->carg++);
143 const char* scentralLongitude = grib_arguments_get_name(h, args, self->carg++);
144 const char* scentralLatitude = grib_arguments_get_name(h, args, self->carg++);
145 const char* sDx = grib_arguments_get_name(h, args, self->carg++);
146 const char* sDy = grib_arguments_get_name(h, args, self->carg++);
147 const char* siScansNegatively = grib_arguments_get_name(h, args, self->carg++);
148 const char* sjScansPositively = grib_arguments_get_name(h, args, self->carg++);
149 const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, self->carg++);
150 const char* salternativeRowScanning = grib_arguments_get_name(h, args, self->carg++);
151
152 if (grib_is_earth_oblate(h)) {
153 grib_context_log(h->context, GRIB_LOG_ERROR, "Polar stereographic only supported for spherical earth.");
154 return GRIB_GEOCALCULUS_PROBLEM;
155 }
156
157 if ((ret = grib_get_double_internal(h, sradius, &radius)) != GRIB_SUCCESS)
158 return ret;
159 if ((ret = grib_get_long_internal(h, snx, &nx)) != GRIB_SUCCESS)
160 return ret;
161 if ((ret = grib_get_long_internal(h, sny, &ny)) != GRIB_SUCCESS)
162 return ret;
163
164 if (iter->nv != nx * ny) {
165 grib_context_log(h->context, GRIB_LOG_ERROR, "Wrong number of points (%ld!=%ldx%ld)", iter->nv, nx, ny);
166 return GRIB_WRONG_GRID;
167 }
168 if ((ret = grib_get_double_internal(h, slatFirstInDegrees, &latFirstInDegrees)) != GRIB_SUCCESS)
169 return ret;
170 if ((ret = grib_get_double_internal(h, slonFirstInDegrees, &lonFirstInDegrees)) != GRIB_SUCCESS)
171 return ret;
172 if ((ret = grib_get_long_internal(h, ssouthPoleOnPlane, &southPoleOnPlane)) != GRIB_SUCCESS)
173 return ret;
174 if ((ret = grib_get_long_internal(h, scentralLongitude, ¢ralLongitudeInDegrees)) != GRIB_SUCCESS)
175 return ret;
176 if ((ret = grib_get_long_internal(h, scentralLatitude, ¢ralLatitudeInDegrees)) != GRIB_SUCCESS)
177 return ret;
178 if ((ret = grib_get_double_internal(h, sDx, &Dx)) != GRIB_SUCCESS)
179 return ret;
180 if ((ret = grib_get_double_internal(h, sDy, &Dy)) != GRIB_SUCCESS)
181 return ret;
182 if ((ret = grib_get_long_internal(h, sjPointsAreConsecutive, &jPointsAreConsecutive)) != GRIB_SUCCESS)
183 return ret;
184 if ((ret = grib_get_long_internal(h, sjScansPositively, &jScansPositively)) != GRIB_SUCCESS)
185 return ret;
186 if ((ret = grib_get_long_internal(h, siScansNegatively, &iScansNegatively)) != GRIB_SUCCESS)
187 return ret;
188 if ((ret = grib_get_long_internal(h, salternativeRowScanning, &alternativeRowScanning)) != GRIB_SUCCESS)
189 return ret;
190
191 centralLongitude = centralLongitudeInDegrees * DEG2RAD;
192 centralLatitude = centralLatitudeInDegrees * DEG2RAD;
193 lonFirst = lonFirstInDegrees * DEG2RAD;
194 latFirst = latFirstInDegrees * DEG2RAD;
195
196 /* Forward projection initialisation */
197 fwd_proj_data.false_northing = 0;
198 fwd_proj_data.false_easting = 0;
199 fwd_proj_data.centre_lon = centralLongitude;
200 fwd_proj_data.centre_lat = centralLatitude;
201 if (centralLatitude < 0)
202 fwd_proj_data.sign = -1.0;
203 else
204 fwd_proj_data.sign = +1.0;
205 fwd_proj_data.ind = 0;
206 if (fabs(fabs(centralLatitude) - PI_OVER_2) > EPSILON) {
207 /* central latitude different from 90 i.e. not north/south polar */
208 fwd_proj_data.ind = 1;
209 con1 = fwd_proj_data.sign * centralLatitude;
210 fwd_proj_data.mcs = cos(con1);
211 fwd_proj_data.tcs = tan(0.5 * (PI_OVER_2 - con1));
212 }
213
214 /* Forward projection from initial lat,lon to initial x,y */
215 con1 = fwd_proj_data.sign * (lonFirst - fwd_proj_data.centre_lon);
216 ts = tan(0.5 * (PI_OVER_2 - fwd_proj_data.sign * latFirst));
217 if (fwd_proj_data.ind)
218 height = radius * fwd_proj_data.mcs * ts / fwd_proj_data.tcs;
219 else
220 height = 2.0 * radius * ts;
221 x0 = fwd_proj_data.sign * height * sin(con1) + fwd_proj_data.false_easting;
222 y0 = -fwd_proj_data.sign * height * cos(con1) + fwd_proj_data.false_northing;
223
224 x0 = -x0;
225 y0 = -y0;
226
227 /* Inverse projection initialisation */
228 inv_proj_data.false_easting = x0;
229 inv_proj_data.false_northing = y0;
230 inv_proj_data.centre_lon = centralLongitude;
231 inv_proj_data.centre_lat = centralLatitude;
232 if (centralLatitude < 0)
233 inv_proj_data.sign = -1.0;
234 else
235 inv_proj_data.sign = +1.0;
236 inv_proj_data.ind = 0;
237 if (fabs(fabs(centralLatitude) - PI_OVER_2) > EPSILON) {
238 inv_proj_data.ind = 1;
239 con1 = inv_proj_data.sign * inv_proj_data.centre_lat;
240 inv_proj_data.mcs = cos(con1);
241 inv_proj_data.tcs = tan(0.5 * (PI_OVER_2 - con1));
242 }
243 self->lats = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
244 if (!self->lats) {
245 grib_context_log(h->context, GRIB_LOG_ERROR, "Error allocating %ld bytes", iter->nv * sizeof(double));
246 return GRIB_OUT_OF_MEMORY;
247 }
248 self->lons = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
249 if (!self->lats) {
250 grib_context_log(h->context, GRIB_LOG_ERROR, "Error allocating %ld bytes", iter->nv * sizeof(double));
251 return GRIB_OUT_OF_MEMORY;
252 }
253 lats = self->lats;
254 lons = self->lons;
255 Dx = iScansNegatively == 0 ? Dx : -Dx;
256 Dy = jScansPositively == 1 ? Dy : -Dy;
257
258 y = 0;
259 for (j = 0; j < ny; j++) {
260 x = 0;
261 for (i = 0; i < nx; i++) {
262 /* Inverse projection from x,y to lat,lon */
263 /* int index =i+j*nx; */
264 double _x = (x - inv_proj_data.false_easting) * inv_proj_data.sign;
265 double _y = (y - inv_proj_data.false_northing) * inv_proj_data.sign;
266 double rh = sqrt(_x * _x + _y * _y);
267 if (inv_proj_data.ind)
268 ts = rh * inv_proj_data.tcs / (radius * inv_proj_data.mcs);
269 else
270 ts = rh / (radius * 2.0);
271 *lats = inv_proj_data.sign * (PI_OVER_2 - 2 * atan(ts));
272 if (rh == 0) {
273 *lons = inv_proj_data.sign * inv_proj_data.centre_lon;
274 }
275 else {
276 double temp = atan2(_x, -_y);
277 *lons = inv_proj_data.sign * temp + inv_proj_data.centre_lon;
278 }
279 *lats = *lats * RAD2DEG;
280 *lons = *lons * RAD2DEG;
281 while (*lons < 0)
282 *lons += 360;
283 while (*lons > 360)
284 *lons -= 360;
285 lons++;
286 lats++;
287
288 x += Dx;
289 }
290 y += Dy;
291 }
292 #if 0
293 /*standardParallel = (southPoleOnPlane == 1) ? -90 : +90;*/
294 if (jPointsAreConsecutive)
295 {
296 x=xFirst;
297 for (i=0;i<nx;i++) {
298 y=yFirst;
299 for (j=0;j<ny;j++) {
300 rho=sqrt(x*x+y*y);
301 if (rho == 0) {
302 /* indeterminate case */
303 *lats = standardParallel;
304 *lons = centralLongitude;
305 }
306 else {
307 c=2*atan2(rho,(2.0*radius));
308 cosc=cos(c);
309 sinc=sin(c);
310 *lats = asin( cosc*sinphi1 + y*sinc*cosphi1/rho ) * RAD2DEG;
311 *lons = (lambda0+atan2(x*sinc, rho*cosphi1*cosc - y*sinphi1*sinc)) * RAD2DEG;
312 }
313 while (*lons<0) *lons += 360;
314 while (*lons>360) *lons -= 360;
315 lons++;
316 lats++;
317
318 y+=Dy;
319 }
320 x+=Dx;
321 }
322 }
323 else
324 {
325 y=yFirst;
326 for (j=0;j<ny;j++) {
327 x=xFirst;
328 for (i=0;i<nx;i++) {
329 /* int index =i+j*nx; */
330 rho=sqrt(x*x+y*y);
331 if (rho == 0) {
332 /* indeterminate case */
333 *lats = standardParallel;
334 *lons = centralLongitude;
335 }
336 else {
337 c=2*atan2(rho,(2.0*radius));
338 cosc=cos(c);
339 sinc=sin(c);
340 *lats = asin( cosc*sinphi1 + y*sinc*cosphi1/rho ) * RAD2DEG;
341 *lons = (lambda0+atan2(x*sinc, rho*cosphi1*cosc - y*sinphi1*sinc)) * RAD2DEG;
342 }
343 while (*lons<0) *lons += 360;
344 while (*lons>360) *lons -= 360;
345 lons++;
346 lats++;
347
348 x+=Dx;
349 }
350 y+=Dy;
351 }
352 }
353 #endif
354 iter->e = -1;
355
356 return ret;
357 }
358
destroy(grib_iterator * i)359 static int destroy(grib_iterator* i)
360 {
361 grib_iterator_polar_stereographic* self = (grib_iterator_polar_stereographic*)i;
362 const grib_context* c = i->h->context;
363
364 grib_context_free(c, self->lats);
365 grib_context_free(c, self->lons);
366 return 1;
367 }
368