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 /**************************************
12 * Enrico Fucile
13 **************************************/
14
15
16 #include "grib_api_internal.h"
17 #include <math.h>
18
19 /*
20 This is used by make_class.pl
21
22 START_CLASS_DEF
23 CLASS = iterator
24 SUPER = grib_iterator_class_gen
25 IMPLEMENTS = destroy
26 IMPLEMENTS = init;next
27 MEMBERS = double *lats
28 MEMBERS = double *lons
29 MEMBERS = long Nj
30 END_CLASS_DEF
31 */
32
33 /* START_CLASS_IMP */
34
35 /*
36
37 Don't edit anything between START_CLASS_IMP and END_CLASS_IMP
38 Instead edit values between START_CLASS_DEF and END_CLASS_DEF
39 or edit "iterator.class" and rerun ./make_class.pl
40
41 */
42
43
44 static void init_class(grib_iterator_class*);
45
46 static int init(grib_iterator* i, grib_handle*, grib_arguments*);
47 static int next(grib_iterator* i, double* lat, double* lon, double* val);
48 static int destroy(grib_iterator* i);
49
50
51 typedef struct grib_iterator_lambert_azimuthal_equal_area
52 {
53 grib_iterator it;
54 /* Members defined in gen */
55 long carg;
56 const char* missingValue;
57 /* Members defined in lambert_azimuthal_equal_area */
58 double* lats;
59 double* lons;
60 long Nj;
61 } grib_iterator_lambert_azimuthal_equal_area;
62
63 extern grib_iterator_class* grib_iterator_class_gen;
64
65 static grib_iterator_class _grib_iterator_class_lambert_azimuthal_equal_area = {
66 &grib_iterator_class_gen, /* super */
67 "lambert_azimuthal_equal_area", /* name */
68 sizeof(grib_iterator_lambert_azimuthal_equal_area), /* size of instance */
69 0, /* inited */
70 &init_class, /* init_class */
71 &init, /* constructor */
72 &destroy, /* destructor */
73 &next, /* Next Value */
74 0, /* Previous Value */
75 0, /* Reset the counter */
76 0, /* has next values */
77 };
78
79 grib_iterator_class* grib_iterator_class_lambert_azimuthal_equal_area = &_grib_iterator_class_lambert_azimuthal_equal_area;
80
81
init_class(grib_iterator_class * c)82 static void init_class(grib_iterator_class* c)
83 {
84 c->previous = (*(c->super))->previous;
85 c->reset = (*(c->super))->reset;
86 c->has_next = (*(c->super))->has_next;
87 }
88 /* END_CLASS_IMP */
89
next(grib_iterator * i,double * lat,double * lon,double * val)90 static int next(grib_iterator* i, double* lat, double* lon, double* val)
91 {
92 grib_iterator_lambert_azimuthal_equal_area* self = (grib_iterator_lambert_azimuthal_equal_area*)i;
93
94 if ((long)i->e >= (long)(i->nv - 1))
95 return 0;
96 i->e++;
97
98 *lat = self->lats[i->e];
99 *lon = self->lons[i->e];
100 *val = i->data[i->e];
101
102 return 1;
103 }
104
init(grib_iterator * iter,grib_handle * h,grib_arguments * args)105 static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
106 {
107 int ret = 0;
108 double *lats, *lons;
109 double lonFirstInDegrees, latFirstInDegrees, lonFirst, latFirst, radius = 0;
110 long nx, ny, standardParallel, centralLongitude;
111 double phi1, lambda0, xFirst, yFirst, x, y, Dx, Dy;
112 double kp, sinphi1, cosphi1;
113 long alternativeRowScanning, iScansNegatively;
114 long jScansPositively, jPointsAreConsecutive;
115 double sinphi, cosphi, cosdlambda, sindlambda;
116 double cosc, sinc;
117 long i, j;
118
119 grib_iterator_lambert_azimuthal_equal_area* self =
120 (grib_iterator_lambert_azimuthal_equal_area*)iter;
121
122 const char* sradius = grib_arguments_get_name(h, args, self->carg++);
123 const char* snx = grib_arguments_get_name(h, args, self->carg++);
124 const char* sny = grib_arguments_get_name(h, args, self->carg++);
125 const char* slatFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
126 const char* slonFirstInDegrees = grib_arguments_get_name(h, args, self->carg++);
127 const char* sstandardParallel = grib_arguments_get_name(h, args, self->carg++);
128 const char* scentralLongitude = grib_arguments_get_name(h, args, self->carg++);
129 const char* sDx = grib_arguments_get_name(h, args, self->carg++);
130 const char* sDy = grib_arguments_get_name(h, args, self->carg++);
131 const char* siScansNegatively = grib_arguments_get_name(h, args, self->carg++);
132 const char* sjScansPositively = grib_arguments_get_name(h, args, self->carg++);
133 const char* sjPointsAreConsecutive = grib_arguments_get_name(h, args, self->carg++);
134 const char* salternativeRowScanning = grib_arguments_get_name(h, args, self->carg++);
135 double c, rho;
136 double epsilon = 1.0e-20;
137 double d2r = acos(0.0) / 90.0;
138
139 if ((ret = grib_get_double_internal(h, sradius, &radius)) != GRIB_SUCCESS) {
140 /* Check if it's an oblate spheroid */
141 if (grib_is_earth_oblate(h))
142 grib_context_log(h->context, GRIB_LOG_ERROR, "Lambert Azimuthal Equal Area only supported for spherical earth.");
143 return ret;
144 }
145
146 if ((ret = grib_get_long_internal(h, snx, &nx)) != GRIB_SUCCESS)
147 return ret;
148 if ((ret = grib_get_long_internal(h, sny, &ny)) != GRIB_SUCCESS)
149 return ret;
150
151 if (iter->nv != nx * ny) {
152 grib_context_log(h->context, GRIB_LOG_ERROR,
153 "Wrong number of points (%ld!=%ldx%ld)",
154 iter->nv, nx, ny);
155 return GRIB_WRONG_GRID;
156 }
157 if ((ret = grib_get_double_internal(h, slatFirstInDegrees, &latFirstInDegrees)) != GRIB_SUCCESS)
158 return ret;
159 if ((ret = grib_get_double_internal(h, slonFirstInDegrees, &lonFirstInDegrees)) != GRIB_SUCCESS)
160 return ret;
161 if ((ret = grib_get_long_internal(h, sstandardParallel, &standardParallel)) != GRIB_SUCCESS)
162 return ret;
163 if ((ret = grib_get_long_internal(h, scentralLongitude, ¢ralLongitude)) != GRIB_SUCCESS)
164 return ret;
165 if ((ret = grib_get_double_internal(h, sDx, &Dx)) != GRIB_SUCCESS)
166 return ret;
167 if ((ret = grib_get_double_internal(h, sDy, &Dy)) != GRIB_SUCCESS)
168 return ret;
169 if ((ret = grib_get_long_internal(h,
170 sjPointsAreConsecutive, &jPointsAreConsecutive)) != GRIB_SUCCESS)
171 return ret;
172 if ((ret = grib_get_long_internal(h, sjScansPositively, &jScansPositively)) != GRIB_SUCCESS)
173 return ret;
174 if ((ret = grib_get_long_internal(h, siScansNegatively, &iScansNegatively)) != GRIB_SUCCESS)
175 return ret;
176 if ((ret = grib_get_long_internal(h,
177 salternativeRowScanning, &alternativeRowScanning)) != GRIB_SUCCESS)
178 return ret;
179
180 lambda0 = d2r * centralLongitude / 1000000;
181 phi1 = d2r * standardParallel / 1000000;
182 latFirst = latFirstInDegrees * d2r;
183 lonFirst = lonFirstInDegrees * d2r;
184
185 cosphi1 = cos(phi1);
186 sinphi1 = sin(phi1);
187
188 Dx = iScansNegatively == 0 ? Dx / 1000 : -Dx / 1000;
189 Dy = jScansPositively == 1 ? Dy / 1000 : -Dy / 1000;
190 self->lats = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
191 if (!self->lats) {
192 grib_context_log(h->context, GRIB_LOG_ERROR,
193 "Error allocating %ld bytes", iter->nv * sizeof(double));
194 return GRIB_OUT_OF_MEMORY;
195 }
196 self->lons = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
197 if (!self->lats) {
198 grib_context_log(h->context, GRIB_LOG_ERROR,
199 "Error allocating %ld bytes", iter->nv * sizeof(double));
200 return GRIB_OUT_OF_MEMORY;
201 }
202 lats = self->lats;
203 lons = self->lons;
204
205 /* compute xFirst,yFirst in metres */
206 sinphi = sin(latFirst);
207 cosphi = cos(latFirst);
208 cosdlambda = cos(lonFirst - lambda0);
209 sindlambda = sin(lonFirst - lambda0);
210 kp = radius * sqrt(2.0 / (1 + sinphi1 * sinphi + cosphi1 * cosphi * cosdlambda));
211 xFirst = kp * cosphi * sindlambda;
212 yFirst = kp * (cosphi1 * sinphi - sinphi1 * cosphi * cosdlambda);
213
214 if (jPointsAreConsecutive) {
215 x = xFirst;
216 for (i = 0; i < nx; i++) {
217 double xsq = x * x;
218 y = yFirst;
219 for (j = 0; j < ny; j++) {
220 rho = sqrt(xsq + y * y);
221 if (rho > epsilon) {
222 c = 2 * asin(rho / (2.0 * radius));
223 cosc = cos(c);
224 sinc = sin(c);
225 *lats = asin(cosc * sinphi1 + y * sinc * cosphi1 / rho) / d2r;
226 *lons = (lambda0 + atan2(x * sinc, rho * cosphi1 * cosc - y * sinphi1 * sinc)) / d2r;
227 }
228 else {
229 *lats = phi1 / d2r;
230 *lons = lambda0 / d2r;
231 }
232 if (*lons < 0)
233 *lons += 360;
234 lons++;
235 lats++;
236
237 y += Dy;
238 }
239 x += Dx;
240 }
241 }
242 else {
243 y = yFirst;
244 for (j = 0; j < ny; j++) {
245 double ysq = y * y;
246 x = xFirst;
247 for (i = 0; i < nx; i++) {
248 rho = sqrt(x * x + ysq);
249 if (rho > epsilon) {
250 c = 2 * asin(rho / (2.0 * radius));
251 cosc = cos(c);
252 sinc = sin(c);
253 *lats = asin(cosc * sinphi1 + y * sinc * cosphi1 / rho) / d2r;
254 *lons = (lambda0 + atan2(x * sinc, rho * cosphi1 * cosc - y * sinphi1 * sinc)) / d2r;
255 }
256 else {
257 *lats = phi1 / d2r;
258 *lons = lambda0 / d2r;
259 }
260 if (*lons < 0)
261 *lons += 360;
262 lons++;
263 lats++;
264
265 x += Dx;
266 }
267 y += Dy;
268 }
269 }
270
271 iter->e = -1;
272
273 return ret;
274 }
275
destroy(grib_iterator * i)276 static int destroy(grib_iterator* i)
277 {
278 grib_iterator_lambert_azimuthal_equal_area* self = (grib_iterator_lambert_azimuthal_equal_area*)i;
279 const grib_context* c = i->h->context;
280
281 grib_context_free(c, self->lats);
282 grib_context_free(c, self->lons);
283 return 1;
284 }
285