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 = previous;next
21    IMPLEMENTS = init;destroy
22    MEMBERS    = double   *las
23    MEMBERS    = double   *los
24    MEMBERS    = long      Ni
25    MEMBERS    = long      Nj
26    MEMBERS    = long iScansNegatively
27    MEMBERS    = long isRotated
28    MEMBERS    = double angleOfRotation
29    MEMBERS    = double southPoleLat
30    MEMBERS    = double southPoleLon
31    MEMBERS    = long jPointsAreConsecutive
32    MEMBERS    = long disableUnrotate
33    END_CLASS_DEF
34 
35  */
36 
37 /* START_CLASS_IMP */
38 
39 /*
40 
41 Don't edit anything between START_CLASS_IMP and END_CLASS_IMP
42 Instead edit values between START_CLASS_DEF and END_CLASS_DEF
43 or edit "iterator.class" and rerun ./make_class.pl
44 
45 */
46 
47 
48 static void init_class(grib_iterator_class*);
49 
50 static int init(grib_iterator* i, grib_handle*, grib_arguments*);
51 static int next(grib_iterator* i, double* lat, double* lon, double* val);
52 static int previous(grib_iterator* ei, double* lat, double* lon, double* val);
53 static int destroy(grib_iterator* i);
54 
55 
56 typedef struct grib_iterator_regular
57 {
58     grib_iterator it;
59     /* Members defined in gen */
60     long carg;
61     const char* missingValue;
62     /* Members defined in regular */
63     double* las;
64     double* los;
65     long Ni;
66     long Nj;
67     long iScansNegatively;
68     long isRotated;
69     double angleOfRotation;
70     double southPoleLat;
71     double southPoleLon;
72     long jPointsAreConsecutive;
73     long disableUnrotate;
74 } grib_iterator_regular;
75 
76 extern grib_iterator_class* grib_iterator_class_gen;
77 
78 static grib_iterator_class _grib_iterator_class_regular = {
79     &grib_iterator_class_gen,      /* super                     */
80     "regular",                     /* name                      */
81     sizeof(grib_iterator_regular), /* size of instance          */
82     0,                             /* inited */
83     &init_class,                   /* init_class */
84     &init,                         /* constructor               */
85     &destroy,                      /* destructor                */
86     &next,                         /* Next Value                */
87     &previous,                     /*  Previous Value           */
88     0,                             /* Reset the counter         */
89     0,                             /* has next values           */
90 };
91 
92 grib_iterator_class* grib_iterator_class_regular = &_grib_iterator_class_regular;
93 
94 
init_class(grib_iterator_class * c)95 static void init_class(grib_iterator_class* c)
96 {
97     c->reset    = (*(c->super))->reset;
98     c->has_next = (*(c->super))->has_next;
99 }
100 /* END_CLASS_IMP */
101 
102 
next(grib_iterator * i,double * lat,double * lon,double * val)103 static int next(grib_iterator* i, double* lat, double* lon, double* val)
104 {
105     grib_iterator_regular* self = (grib_iterator_regular*)i;
106 
107     if ((long)i->e >= (long)(i->nv - 1))
108         return 0;
109 
110     i->e++;
111 
112     *lat = self->las[(long)floor(i->e / self->Ni)];
113     *lon = self->los[(long)i->e % self->Ni];
114     *val = i->data[i->e];
115 
116     return 1;
117 }
118 
previous(grib_iterator * i,double * lat,double * lon,double * val)119 static int previous(grib_iterator* i, double* lat, double* lon, double* val)
120 {
121     grib_iterator_regular* self = (grib_iterator_regular*)i;
122 
123     if (i->e < 0)
124         return 0;
125     *lat = self->las[(long)floor(i->e / self->Ni)];
126     *lon = self->los[i->e % self->Ni];
127     *val = i->data[i->e];
128     i->e--;
129 
130     return 1;
131 }
132 
destroy(grib_iterator * i)133 static int destroy(grib_iterator* i)
134 {
135     grib_iterator_regular* self = (grib_iterator_regular*)i;
136     const grib_context* c       = i->h->context;
137     grib_context_free(c, self->las);
138     grib_context_free(c, self->los);
139     return GRIB_SUCCESS;
140 }
141 
init(grib_iterator * i,grib_handle * h,grib_arguments * args)142 static int init(grib_iterator* i, grib_handle* h, grib_arguments* args)
143 {
144     grib_iterator_regular* self = (grib_iterator_regular*)i;
145     int ret                     = GRIB_SUCCESS;
146 
147     long Ni; /* Number of points along a parallel = Nx */
148     long Nj; /* Number of points along a meridian = Ny */
149     double idir, lon1, lon2;
150     long loi;
151 
152     const char* s_lon1      = grib_arguments_get_name(h, args, self->carg++);
153     const char* s_idir      = grib_arguments_get_name(h, args, self->carg++);
154     const char* s_Ni        = grib_arguments_get_name(h, args, self->carg++);
155     const char* s_Nj        = grib_arguments_get_name(h, args, self->carg++);
156     const char* s_iScansNeg = grib_arguments_get_name(h, args, self->carg++);
157 
158     if ((ret = grib_get_double_internal(h, s_lon1, &lon1)))
159         return ret;
160     if ((ret = grib_get_double_internal(h, "longitudeOfLastGridPointInDegrees", &lon2)))
161         return ret;
162     if ((ret = grib_get_double_internal(h, s_idir, &idir)))
163         return ret;
164     if ((ret = grib_get_long_internal(h, s_Ni, &Ni)))
165         return ret;
166     if (grib_is_missing(h, s_Ni, &ret) && ret == GRIB_SUCCESS) {
167         grib_context_log(h->context, GRIB_LOG_ERROR, "Key %s cannot be 'missing' for a regular grid!", s_Ni);
168         return GRIB_WRONG_GRID;
169     }
170 
171     if ((ret = grib_get_long_internal(h, s_Nj, &Nj)))
172         return ret;
173     if (grib_is_missing(h, s_Nj, &ret) && ret == GRIB_SUCCESS) {
174         grib_context_log(h->context, GRIB_LOG_ERROR, "Key %s cannot be 'missing' for a regular grid!", s_Nj);
175         return GRIB_WRONG_GRID;
176     }
177 
178     if ((ret = grib_get_long_internal(h, s_iScansNeg, &self->iScansNegatively)))
179         return ret;
180 
181     /* GRIB-801: Careful of case with a single point! Ni==1 */
182     if (Ni > 1) {
183         /* Note: If first and last longitudes are equal I assume you wanna go round the globe */
184         if (self->iScansNegatively) {
185             if (lon1 > lon2) {
186                 idir = (lon1 - lon2) / (Ni - 1);
187             }
188             else {
189                 idir = (lon1 + 360.0 - lon2) / (Ni - 1);
190             }
191         }
192         else {
193             if (lon2 > lon1) {
194                 idir = (lon2 - lon1) / (Ni - 1);
195             }
196             else {
197                 idir = (lon2 + 360.0 - lon1) / (Ni - 1);
198             }
199         }
200     }
201     if (self->iScansNegatively) {
202         idir = -idir;
203     }
204     else {
205         if (lon1 + (Ni - 2) * idir > 360)
206             lon1 -= 360;
207         /*See ECC-704, GRIB-396*/
208         /*else if ( (lon1+(Ni-1)*idir)-360 > epsilon ){
209             idir=360.0/(float)Ni;
210         }*/
211     }
212 
213     self->Ni = Ni;
214     self->Nj = Nj;
215 
216     self->las = (double*)grib_context_malloc(h->context, Nj * sizeof(double));
217     self->los = (double*)grib_context_malloc(h->context, Ni * sizeof(double));
218 
219     for (loi = 0; loi < Ni; loi++) {
220         self->los[loi] = lon1;
221         lon1 += idir;
222     }
223 
224     return ret;
225 }
226