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_space_view
47 {
48 grib_iterator it;
49 /* Members defined in gen */
50 long carg;
51 const char* missingValue;
52 /* Members defined in space_view */
53 double* lats;
54 double* lons;
55 long Nj;
56 } grib_iterator_space_view;
57
58 extern grib_iterator_class* grib_iterator_class_gen;
59
60 static grib_iterator_class _grib_iterator_class_space_view = {
61 &grib_iterator_class_gen, /* super */
62 "space_view", /* name */
63 sizeof(grib_iterator_space_view), /* 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_space_view = &_grib_iterator_class_space_view;
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_space_view* self = (grib_iterator_space_view*)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 #define RAD2DEG 57.29577951308232087684 /* 180 over pi */
101 #define DEG2RAD 0.01745329251994329576 /* pi over 180 */
102
init(grib_iterator * iter,grib_handle * h,grib_arguments * args)103 static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
104 {
105 /* REFERENCE:
106 * LRIT/HRIT Global Specification (CGMS 03, Issue 2.6, 12.08.1999)
107 */
108 int ret = GRIB_SUCCESS;
109 double *lats, *lons; /* arrays of latitudes and longitudes */
110 double latOfSubSatellitePointInDegrees, lonOfSubSatellitePointInDegrees;
111 double orientationInDegrees, nrInRadiusOfEarth;
112 double radius = 0, xpInGridLengths = 0, ypInGridLengths = 0;
113 long nx, ny, earthIsOblate = 0;
114 long alternativeRowScanning, iScansNegatively;
115 long Xo, Yo, jScansPositively, jPointsAreConsecutive, i;
116
117 double major = 0, minor = 0, r_eq, r_pol, height;
118 double lap, lop, orient_angle, angular_size;
119 double xp, yp, dx, dy, rx, ry, x, y;
120 double cos_x, cos_y, sin_x, sin_y;
121 double factor_1, factor_2, tmp1, Sd, Sn, Sxy, S1, S2, S3;
122 int x0, y0, ix, iy;
123 double *s_x, *c_x; /* arrays storing sin and cos values */
124 size_t array_size = (iter->nv * sizeof(double));
125
126 grib_iterator_space_view* self = (grib_iterator_space_view*)iter;
127
128 const char* sradius = grib_arguments_get_name(h, args, self->carg++);
129 const char* sEarthIsOblate = grib_arguments_get_name(h, args, self->carg++);
130 const char* sMajorAxisInMetres = grib_arguments_get_name(h, args, self->carg++);
131 const char* sMinorAxisInMetres = grib_arguments_get_name(h, args, self->carg++);
132 const char* snx = grib_arguments_get_name(h, args, self->carg++);
133 const char* sny = grib_arguments_get_name(h, args, self->carg++);
134 const char* sLatOfSubSatellitePointInDegrees = grib_arguments_get_name(h, args, self->carg++);
135 const char* sLonOfSubSatellitePointInDegrees = grib_arguments_get_name(h, args, self->carg++);
136 const char* sDx = grib_arguments_get_name(h, args, self->carg++);
137 const char* sDy = grib_arguments_get_name(h, args, self->carg++);
138 const char* sXpInGridLengths = grib_arguments_get_name(h, args, self->carg++);
139 const char* sYpInGridLengths = grib_arguments_get_name(h, args, self->carg++);
140 const char* sOrientationInDegrees = grib_arguments_get_name(h, args, self->carg++);
141 const char* sNrInRadiusOfEarth = grib_arguments_get_name(h, args, self->carg++);
142 const char* sXo = grib_arguments_get_name(h, args, self->carg++);
143 const char* sYo = grib_arguments_get_name(h, args, self->carg++);
144
145 const char* siScansNegatively = grib_arguments_get_name(h, args, self->carg++);
146 const char* sjScansPositively = grib_arguments_get_name(h, args, self->carg++);
147 const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, self->carg++);
148 const char* sAlternativeRowScanning = grib_arguments_get_name(h, args, self->carg++);
149
150 if ((ret = grib_get_long_internal(h, snx, &nx)) != GRIB_SUCCESS)
151 return ret;
152 if ((ret = grib_get_long_internal(h, sny, &ny)) != GRIB_SUCCESS)
153 return ret;
154 if ((ret = grib_get_long_internal(h, sEarthIsOblate, &earthIsOblate)) != GRIB_SUCCESS)
155 return ret;
156
157 if (earthIsOblate) {
158 if ((ret = grib_get_double_internal(h, sMajorAxisInMetres, &major)) != GRIB_SUCCESS)
159 return ret;
160 if ((ret = grib_get_double_internal(h, sMinorAxisInMetres, &minor)) != GRIB_SUCCESS)
161 return ret;
162 }
163 else {
164 if ((ret = grib_get_double_internal(h, sradius, &radius)) != GRIB_SUCCESS)
165 return ret;
166 }
167
168 if (iter->nv != nx * ny) {
169 grib_context_log(h->context, GRIB_LOG_ERROR, "Wrong number of points (%ld!=%ldx%ld)", iter->nv, nx, ny);
170 return GRIB_WRONG_GRID;
171 }
172 if ((ret = grib_get_double_internal(h, sLatOfSubSatellitePointInDegrees, &latOfSubSatellitePointInDegrees)) != GRIB_SUCCESS)
173 return ret;
174 if ((ret = grib_get_double_internal(h, sLonOfSubSatellitePointInDegrees, &lonOfSubSatellitePointInDegrees)) != GRIB_SUCCESS)
175 return ret;
176 if ((ret = grib_get_double_internal(h, sDx, &dx)) != GRIB_SUCCESS)
177 return ret;
178 if ((ret = grib_get_double_internal(h, sDy, &dy)) != GRIB_SUCCESS)
179 return ret;
180 if ((ret = grib_get_double_internal(h, sXpInGridLengths, &xpInGridLengths)) != GRIB_SUCCESS)
181 return ret;
182 if ((ret = grib_get_double_internal(h, sYpInGridLengths, &ypInGridLengths)) != GRIB_SUCCESS)
183 return ret;
184 if ((ret = grib_get_double_internal(h, sOrientationInDegrees, &orientationInDegrees)) != GRIB_SUCCESS)
185 return ret;
186
187 /* Orthographic not supported. This happens when Nr (camera altitude) is missing */
188 if (grib_is_missing(h, sNrInRadiusOfEarth, &ret)) {
189 grib_context_log(h->context, GRIB_LOG_ERROR, "Orthographic view (Nr missing) not supported");
190 return GRIB_NOT_IMPLEMENTED;
191 }
192 if ((ret = grib_get_double_internal(h, sNrInRadiusOfEarth, &nrInRadiusOfEarth)) != GRIB_SUCCESS)
193 return ret;
194
195 if ((ret = grib_get_long_internal(h, sXo, &Xo)) != GRIB_SUCCESS)
196 return ret;
197 if ((ret = grib_get_long_internal(h, sYo, &Yo)) != GRIB_SUCCESS)
198 return ret;
199 if ((ret = grib_get_long_internal(h, sjPointsAreConsecutive, &jPointsAreConsecutive)) != GRIB_SUCCESS)
200 return ret;
201 if ((ret = grib_get_long_internal(h, sjScansPositively, &jScansPositively)) != GRIB_SUCCESS)
202 return ret;
203 if ((ret = grib_get_long_internal(h, siScansNegatively, &iScansNegatively)) != GRIB_SUCCESS)
204 return ret;
205 if ((ret = grib_get_long_internal(h, sAlternativeRowScanning, &alternativeRowScanning)) != GRIB_SUCCESS)
206 return ret;
207
208 if (earthIsOblate) {
209 r_eq = major; /* In km */
210 r_pol = minor;
211 }
212 else {
213 r_eq = r_pol = radius * 0.001; /*conv to km*/
214 }
215 angular_size = 2.0 * asin(1.0 / nrInRadiusOfEarth);
216 height = nrInRadiusOfEarth * r_eq;
217
218 lap = latOfSubSatellitePointInDegrees;
219 lop = lonOfSubSatellitePointInDegrees;
220 lap *= 1e-6; /* default scaling factor */
221 lop *= 1e-6;
222 if (lap != 0.0)
223 return GRIB_NOT_IMPLEMENTED;
224 /*lap *= DEG2RAD;*/
225 lop *= DEG2RAD;
226
227 orient_angle = orientationInDegrees;
228 if (orient_angle != 0.0)
229 return GRIB_NOT_IMPLEMENTED;
230
231 xp = xpInGridLengths;
232 yp = ypInGridLengths;
233 x0 = Xo;
234 y0 = Yo;
235
236 rx = angular_size / dx;
237 ry = (r_pol / r_eq) * angular_size / dy;
238
239 self->lats = (double*)grib_context_malloc(h->context, array_size);
240 if (!self->lats) {
241 grib_context_log(h->context, GRIB_LOG_ERROR, "Error allocating %ld bytes", array_size);
242 return GRIB_OUT_OF_MEMORY;
243 }
244 self->lons = (double*)grib_context_malloc(h->context, array_size);
245 if (!self->lats) {
246 grib_context_log(h->context, GRIB_LOG_ERROR, "Error allocating %ld bytes", array_size);
247 return GRIB_OUT_OF_MEMORY;
248 }
249 lats = self->lats;
250 lons = self->lons;
251
252 if (!iScansNegatively) {
253 xp = xp - x0;
254 }
255 else {
256 xp = (nx - 1) - (xp - x0);
257 }
258 if (jScansPositively) {
259 yp = yp - y0;
260 }
261 else {
262 yp = (ny - 1) - (yp - y0);
263 }
264 i = 0;
265 factor_2 = (r_eq / r_pol) * (r_eq / r_pol);
266 factor_1 = height * height - r_eq * r_eq;
267
268 /* Store array of sin and cosine values to avoid recalculation */
269 s_x = (double*)grib_context_malloc(h->context, nx * sizeof(double));
270 if (!s_x) {
271 grib_context_log(h->context, GRIB_LOG_ERROR, "Error allocating %ld bytes", nx * sizeof(double));
272 return GRIB_OUT_OF_MEMORY;
273 }
274 c_x = (double*)grib_context_malloc(h->context, nx * sizeof(double));
275 if (!c_x) {
276 grib_context_log(h->context, GRIB_LOG_ERROR, "Error allocating %ld bytes", nx * sizeof(double));
277 return GRIB_OUT_OF_MEMORY;
278 }
279
280 for (ix = 0; ix < nx; ix++) {
281 x = (ix - xp) * rx;
282 s_x[ix] = sin(x);
283 c_x[ix] = sqrt(1.0 - s_x[ix] * s_x[ix]);
284 }
285
286 /*for (iy = 0; iy < ny; iy++) {*/
287 for (iy = ny - 1; iy >= 0; --iy) {
288 y = (iy - yp) * ry;
289 sin_y = sin(y);
290 cos_y = sqrt(1.0 - sin_y * sin_y);
291
292 tmp1 = (1 + (factor_2 - 1.0) * sin_y * sin_y);
293
294 for (ix = 0; ix < nx; ix++, i++) {
295 /*x = (ix - xp) * rx;*/
296 /* Use sin/cos previously computed */
297 sin_x = s_x[ix];
298 cos_x = c_x[ix];
299
300 Sd = height * cos_x * cos_y;
301 Sd = Sd * Sd - tmp1 * factor_1;
302 if (Sd <= 0.0) { /* outside of view */
303 lats[i] = lons[i] = 0; /* TODO: error? */
304 }
305 else {
306 Sd = sqrt(Sd);
307 Sn = (height * cos_x * cos_y - Sd) / tmp1;
308 S1 = height - Sn * cos_x * cos_y;
309 S2 = Sn * sin_x * cos_y;
310 S3 = Sn * sin_y;
311 Sxy = sqrt(S1 * S1 + S2 * S2);
312 lons[i] = atan(S2 / S1) * (RAD2DEG) + lop;
313 lats[i] = atan(factor_2 * S3 / Sxy) * (RAD2DEG);
314 /*fprintf(stderr, "lat=%g lon=%g\n", lats[i], lons[i]);*/
315 }
316 while (lons[i] < 0)
317 lons[i] += 360;
318 while (lons[i] > 360)
319 lons[i] -= 360;
320 }
321 }
322 grib_context_free(h->context, s_x);
323 grib_context_free(h->context, c_x);
324 iter->e = -1;
325
326 return ret;
327 }
328
destroy(grib_iterator * i)329 static int destroy(grib_iterator* i)
330 {
331 grib_iterator_space_view* self = (grib_iterator_space_view*)i;
332 const grib_context* c = i->h->context;
333
334 grib_context_free(c, self->lats);
335 grib_context_free(c, self->lons);
336 return 1;
337 }
338