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, &centralLongitude)) != 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