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
13 /*
14 This is used by make_class.pl
15
16 START_CLASS_DEF
17 CLASS = iterator
18 SUPER = grib_iterator_class_regular
19 IMPLEMENTS = init;next
20 END_CLASS_DEF
21
22 */
23
24 /* START_CLASS_IMP */
25
26 /*
27
28 Don't edit anything between START_CLASS_IMP and END_CLASS_IMP
29 Instead edit values between START_CLASS_DEF and END_CLASS_DEF
30 or edit "iterator.class" and rerun ./make_class.pl
31
32 */
33
34
35 static void init_class(grib_iterator_class*);
36
37 static int init(grib_iterator* i, grib_handle*, grib_arguments*);
38 static int next(grib_iterator* i, double* lat, double* lon, double* val);
39
40
41 typedef struct grib_iterator_latlon
42 {
43 grib_iterator it;
44 /* Members defined in gen */
45 long carg;
46 const char* missingValue;
47 /* Members defined in regular */
48 double* las;
49 double* los;
50 long Ni;
51 long Nj;
52 long iScansNegatively;
53 long isRotated;
54 double angleOfRotation;
55 double southPoleLat;
56 double southPoleLon;
57 long jPointsAreConsecutive;
58 long disableUnrotate;
59 /* Members defined in latlon */
60 } grib_iterator_latlon;
61
62 extern grib_iterator_class* grib_iterator_class_regular;
63
64 static grib_iterator_class _grib_iterator_class_latlon = {
65 &grib_iterator_class_regular, /* super */
66 "latlon", /* name */
67 sizeof(grib_iterator_latlon), /* size of instance */
68 0, /* inited */
69 &init_class, /* init_class */
70 &init, /* constructor */
71 0, /* destructor */
72 &next, /* Next Value */
73 0, /* Previous Value */
74 0, /* Reset the counter */
75 0, /* has next values */
76 };
77
78 grib_iterator_class* grib_iterator_class_latlon = &_grib_iterator_class_latlon;
79
80
init_class(grib_iterator_class * c)81 static void init_class(grib_iterator_class* c)
82 {
83 c->previous = (*(c->super))->previous;
84 c->reset = (*(c->super))->reset;
85 c->has_next = (*(c->super))->has_next;
86 }
87 /* END_CLASS_IMP */
88
next(grib_iterator * iter,double * lat,double * lon,double * val)89 static int next(grib_iterator* iter, double* lat, double* lon, double* val)
90 {
91 /* GRIB-238: Support rotated lat/lon grids */
92 double ret_lat, ret_lon, ret_val;
93 grib_iterator_latlon* self = (grib_iterator_latlon*)iter;
94
95 if ((long)iter->e >= (long)(iter->nv - 1))
96 return 0;
97
98 iter->e++;
99
100 /* Assumptions:
101 * All rows scan in the same direction (alternativeRowScanning==0)
102 */
103 if (!self->jPointsAreConsecutive) {
104 /* Adjacent points in i (x) direction are consecutive */
105 ret_lat = self->las[(long)floor(iter->e / self->Ni)];
106 ret_lon = self->los[(long)iter->e % self->Ni];
107 ret_val = iter->data[iter->e];
108 }
109 else {
110 /* Adjacent points in j (y) direction is consecutive */
111 ret_lon = self->los[(long)iter->e / self->Nj];
112 ret_lat = self->las[(long)floor(iter->e % self->Nj)];
113 ret_val = iter->data[iter->e];
114 }
115
116 /* See ECC-808: Some users want to disable the unrotate */
117 if (self->isRotated && !self->disableUnrotate) {
118 double new_lat = 0, new_lon = 0;
119 unrotate(ret_lat, ret_lon,
120 self->angleOfRotation, self->southPoleLat, self->southPoleLon,
121 &new_lat, &new_lon);
122 ret_lat = new_lat;
123 ret_lon = new_lon;
124 }
125
126 *lat = ret_lat;
127 *lon = ret_lon;
128 *val = ret_val;
129 return 1;
130 }
131
init(grib_iterator * iter,grib_handle * h,grib_arguments * args)132 static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
133 {
134 grib_iterator_latlon* self = (grib_iterator_latlon*)iter;
135 int err = 0;
136 double jdir;
137 double lat1=0, lat2=0, north=0, south=0;
138 long jScansPositively;
139 long lai;
140
141 const char* s_lat1 = grib_arguments_get_name(h, args, self->carg++);
142 const char* s_jdir = grib_arguments_get_name(h, args, self->carg++);
143 const char* s_jScansPos = grib_arguments_get_name(h, args, self->carg++);
144 const char* s_jPtsConsec = grib_arguments_get_name(h, args, self->carg++);
145 const char* s_isRotatedGrid = grib_arguments_get_name(h, args, self->carg++);
146 const char* s_angleOfRotation= grib_arguments_get_name(h, args, self->carg++);
147 const char* s_latSouthernPole= grib_arguments_get_name(h, args, self->carg++);
148 const char* s_lonSouthernPole= grib_arguments_get_name(h, args, self->carg++);
149
150 self->angleOfRotation = 0;
151 self->isRotated = 0;
152 self->southPoleLat = 0;
153 self->southPoleLon = 0;
154 self->disableUnrotate = 0; /* unrotate enabled by default */
155
156 if ((err = grib_get_long(h, s_isRotatedGrid, &self->isRotated)))
157 return err;
158 if (self->isRotated) {
159 if ((err = grib_get_double_internal(h, s_angleOfRotation, &self->angleOfRotation)))
160 return err;
161 if ((err = grib_get_double_internal(h, s_latSouthernPole, &self->southPoleLat)))
162 return err;
163 if ((err = grib_get_double_internal(h, s_lonSouthernPole, &self->southPoleLon)))
164 return err;
165 }
166
167 if ((err = grib_get_double_internal(h, s_lat1, &lat1)))
168 return err;
169 if ((err = grib_get_double_internal(h, "latitudeLastInDegrees", &lat2)))
170 return err;
171 if ((err = grib_get_double_internal(h, s_jdir, &jdir)))
172 return err;
173 if ((err = grib_get_long_internal(h, s_jScansPos, &jScansPositively)))
174 return err;
175 if ((err = grib_get_long_internal(h, s_jPtsConsec, &self->jPointsAreConsecutive)))
176 return err;
177 if ((err = grib_get_long(h, "iteratorDisableUnrotate", &self->disableUnrotate)))
178 return err;
179
180 /* ECC-984: If jDirectionIncrement is missing, then we cannot use it (See jDirectionIncrementGiven) */
181 /* So try to compute the increment */
182 if (grib_is_missing(h, s_jdir, &err) && err == GRIB_SUCCESS) {
183 const long Nj = self->Nj;
184 Assert(Nj > 1);
185 if (lat1 > lat2) {
186 jdir = (lat1 - lat2) / (Nj - 1);
187 }
188 else {
189 jdir = (lat2 - lat1) / (Nj - 1);
190 }
191 grib_context_log(h->context, GRIB_LOG_INFO,
192 "%s is missing (See jDirectionIncrementGiven). Using value of %.6f obtained from La1, La2 and Nj", s_jdir, jdir);
193 }
194
195 if (jScansPositively) {
196 north = lat2;
197 south = lat1;
198 jdir = -jdir;
199 } else {
200 north = lat1;
201 south = lat2;
202 }
203 if (south > north) {
204 grib_context_log(h->context, GRIB_LOG_ERROR,
205 "First and last latitudes are inconsistent with scanning order: lat1=%g, lat2=%g jScansPositively=%ld",
206 lat1, lat2, jScansPositively);
207 return GRIB_WRONG_GRID;
208 }
209
210 for (lai = 0; lai < self->Nj; lai++) {
211 self->las[lai] = lat1;
212 lat1 -= jdir;
213 }
214
215 iter->e = -1;
216 return err;
217 }
218