1 /*
2 * Copyright 2005-2018 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 grib_iterator it;
43 /* Members defined in gen */
44 long carg;
45 const char* missingValue;
46 /* Members defined in regular */
47 double *las;
48 double *los;
49 long nap;
50 long nam;
51 long iScansNegatively;
52 long isRotated;
53 double angleOfRotation;
54 double southPoleLat;
55 double southPoleLon;
56 /* Members defined in latlon */
57 } grib_iterator_latlon;
58
59 extern grib_iterator_class* grib_iterator_class_regular;
60
61 static grib_iterator_class _grib_iterator_class_latlon = {
62 &grib_iterator_class_regular, /* super */
63 "latlon", /* name */
64 sizeof(grib_iterator_latlon),/* size of instance */
65 0, /* inited */
66 &init_class, /* init_class */
67 &init, /* constructor */
68 0, /* destructor */
69 &next, /* Next Value */
70 0, /* Previous Value */
71 0, /* Reset the counter */
72 0, /* has next values */
73 };
74
75 grib_iterator_class* grib_iterator_class_latlon = &_grib_iterator_class_latlon;
76
77
init_class(grib_iterator_class * c)78 static void init_class(grib_iterator_class* c)
79 {
80 c->previous = (*(c->super))->previous;
81 c->reset = (*(c->super))->reset;
82 c->has_next = (*(c->super))->has_next;
83 }
84 /* END_CLASS_IMP */
85
86 #define RAD2DEG 57.29577951308232087684 /* 180 over pi */
87 #define DEG2RAD 0.01745329251994329576 /* pi over 180 */
88
unrotate(grib_handle * h,const double inlat,const double inlon,const double angleOfRot,const double southPoleLat,const double southPoleLon,double * outlat,double * outlon)89 void unrotate(grib_handle* h,
90 const double inlat, const double inlon,
91 const double angleOfRot, const double southPoleLat, const double southPoleLon,
92 double* outlat, double* outlon)
93 {
94 /* Algorithm taken from ecKit */
95 const double lon_x = inlon;
96 const double lat_y = inlat;
97 /* First convert the data point from spherical lat lon to (x',y',z') */
98 double latr = lat_y * DEG2RAD;
99 double lonr = lon_x * DEG2RAD;
100 double xd = cos(lonr)*cos(latr);
101 double yd = sin(lonr)*cos(latr);
102 double zd = sin(latr);
103
104 double t = -(90.0 + southPoleLat);
105 double o = -southPoleLon;
106
107 double sin_t = sin(DEG2RAD * t);
108 double cos_t = cos(DEG2RAD * t);
109 double sin_o = sin(DEG2RAD * o);
110 double cos_o = cos(DEG2RAD * o);
111
112 double x = cos_t*cos_o*xd + sin_o*yd + sin_t*cos_o*zd;
113 double y = -cos_t*sin_o*xd + cos_o*yd - sin_t*sin_o*zd;
114 double z = -sin_t*xd + cos_t*zd;
115
116 double ret_lat=0, ret_lon=0;
117
118 /* Then convert back to 'normal' (lat,lon)
119 * Uses arcsin, to convert back to degrees, put in range -1 to 1 in case of slight rounding error
120 * avoid error on calculating e.g. asin(1.00000001) */
121 if (z > 1.0) z = 1.0;
122 if (z < -1.0) z = -1.0;
123
124 ret_lat = asin(z) * RAD2DEG;
125 ret_lon = atan2(y, x) * RAD2DEG;
126
127 /* Still get a very small rounding error, round to 6 decimal places */
128 /* TODO: roundf implementation missing on MSVC! */
129 #ifndef GRIB_ON_WINDOWS
130 ret_lat = roundf( ret_lat * 1000000.0 )/1000000.0;
131 ret_lon = roundf( ret_lon * 1000000.0 )/1000000.0;
132 #endif
133
134 ret_lon -= angleOfRot;
135
136 /* Make sure ret_lon is in range*/
137 /*
138 while (ret_lon < lonmin_) ret_lon += 360.0;
139 while (ret_lon >= lonmax_) ret_lon -= 360.0;
140 */
141 *outlat = ret_lat;
142 *outlon = ret_lon;
143 }
144
next(grib_iterator * i,double * lat,double * lon,double * val)145 static int next(grib_iterator* i, double *lat, double *lon, double *val)
146 {
147 /* GRIB-238: Support rotated lat/lon grids */
148
149 double ret_lat, ret_lon, ret_val;
150 grib_iterator_latlon* self = (grib_iterator_latlon*)i;
151
152 if((long)i->e >= (long)(i->nv-1)) return 0;
153
154 i->e++;
155
156 ret_lat = self->las[(long)floor(i->e/self->nap)];
157 ret_lon = self->los[(long)i->e%self->nap];
158 ret_val = i->data[i->e];
159
160 if (self->isRotated)
161 {
162 double new_lat = 0, new_lon = 0;
163 unrotate(i->h, ret_lat, ret_lon,
164 self->angleOfRotation, self->southPoleLat, self->southPoleLon,
165 &new_lat, &new_lon);
166 ret_lat = new_lat;
167 ret_lon = new_lon;
168 }
169
170 *lat = ret_lat;
171 *lon = ret_lon;
172 *val = ret_val;
173 return 1;
174 }
175
init(grib_iterator * i,grib_handle * h,grib_arguments * args)176 static int init(grib_iterator* i,grib_handle* h,grib_arguments* args)
177 {
178 grib_iterator_latlon* self = (grib_iterator_latlon*)i;
179 int ret = GRIB_SUCCESS;
180 double jdir;
181 double laf;
182 long jScansPositively;
183 long lai;
184
185 const char* latofirst = grib_arguments_get_name(h,args,self->carg++);
186 const char* jdirec = grib_arguments_get_name(h,args,self->carg++);
187 const char* s_jScansPositively = grib_arguments_get_name(h,args,self->carg++);
188 self->angleOfRotation = 0;
189 self->isRotated = 0;
190 self->southPoleLat = 0;
191 self->southPoleLon = 0;
192
193 if ((ret = grib_get_long(h, "is_rotated_grid", &self->isRotated))) return ret;
194 if (self->isRotated) {
195 if ((ret = grib_get_double_internal(h,"angleOfRotation", &self->angleOfRotation))) return ret;
196 if ((ret = grib_get_double_internal(h,"latitudeOfSouthernPoleInDegrees", &self->southPoleLat))) return ret;
197 if ((ret = grib_get_double_internal(h,"longitudeOfSouthernPoleInDegrees", &self->southPoleLon))) return ret;
198 }
199
200 #if 0
201 ret = grib_get_double(h,"angleOfRotation", &self->angleOfRotation);
202 if (ret != GRIB_SUCCESS) {
203 if (ret == GRIB_NOT_FOUND) {
204 ret = GRIB_SUCCESS;
205 } else {
206 return ret;
207 }
208 }
209 else {
210 self->isRotated = 1;
211 if ((ret = grib_get_double_internal(h,"latitudeOfSouthernPoleInDegrees", &self->southPoleLat))) return ret;
212 if ((ret = grib_get_double_internal(h,"longitudeOfSouthernPoleInDegrees", &self->southPoleLon))) return ret;
213 }
214 #endif
215 if((ret = grib_get_double_internal(h,latofirst, &laf))) return ret;
216 if((ret = grib_get_double_internal(h,jdirec, &jdir))) return ret;
217 if((ret = grib_get_long_internal(h,s_jScansPositively,&jScansPositively)))
218 return ret;
219
220 if (jScansPositively) jdir=-jdir;
221
222 for( lai = 0; lai < self->nam; lai++ ) {
223 self->las[lai] = laf;
224 laf -= jdir ;
225 }
226
227 i->e = -1;
228 return ret;
229 }
230