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 #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 = destroy
21 IMPLEMENTS = init;next
22 MEMBERS = double *lats
23 MEMBERS = double *lons
24 MEMBERS = long nam
25 END_CLASS_DEF
26 */
27
28 /* START_CLASS_IMP */
29
30 /*
31
32 Don't edit anything between START_CLASS_IMP and END_CLASS_IMP
33 Instead edit values between START_CLASS_DEF and END_CLASS_DEF
34 or edit "iterator.class" and rerun ./make_class.pl
35
36 */
37
38
39 static void init_class (grib_iterator_class*);
40
41 static int init (grib_iterator* i,grib_handle*,grib_arguments*);
42 static int next (grib_iterator* i, double *lat, double *lon, double *val);
43 static int destroy (grib_iterator* i);
44
45
46 typedef struct grib_iterator_lambert_conformal{
47 grib_iterator it;
48 /* Members defined in gen */
49 long carg;
50 const char* missingValue;
51 /* Members defined in lambert_conformal */
52 double *lats;
53 double *lons;
54 long nam;
55 } grib_iterator_lambert_conformal;
56
57 extern grib_iterator_class* grib_iterator_class_gen;
58
59 static grib_iterator_class _grib_iterator_class_lambert_conformal = {
60 &grib_iterator_class_gen, /* super */
61 "lambert_conformal", /* name */
62 sizeof(grib_iterator_lambert_conformal),/* size of instance */
63 0, /* inited */
64 &init_class, /* init_class */
65 &init, /* constructor */
66 &destroy, /* destructor */
67 &next, /* Next Value */
68 0, /* Previous Value */
69 0, /* Reset the counter */
70 0, /* has next values */
71 };
72
73 grib_iterator_class* grib_iterator_class_lambert_conformal = &_grib_iterator_class_lambert_conformal;
74
75
init_class(grib_iterator_class * c)76 static void init_class(grib_iterator_class* c)
77 {
78 c->previous = (*(c->super))->previous;
79 c->reset = (*(c->super))->reset;
80 c->has_next = (*(c->super))->has_next;
81 }
82 /* END_CLASS_IMP */
83
next(grib_iterator * i,double * lat,double * lon,double * val)84 static int next(grib_iterator* i, double *lat, double *lon, double *val)
85 {
86 grib_iterator_lambert_conformal* self = (grib_iterator_lambert_conformal*)i;
87
88 if((long)i->e >= (long)(i->nv-1))
89 return 0;
90 i->e++;
91
92 *lat = self->lats[i->e];
93 *lon = self->lons[i->e];
94 *val = i->data[i->e];
95
96 return 1;
97 }
98
99 #ifndef M_PI
100 #define M_PI 3.14159265358979323846 /* Whole pie */
101 #endif
102
103 #ifndef M_PI_2
104 #define M_PI_2 1.57079632679489661923 /* Half a pie */
105 #endif
106
107 #ifndef M_PI_4
108 #define M_PI_4 0.78539816339744830962 /* Quarter of a pie */
109 #endif
110
111 #define RAD2DEG 57.29577951308232087684 /* 180 over pi */
112 #define DEG2RAD 0.01745329251994329576 /* pi over 180 */
113
init(grib_iterator * iter,grib_handle * h,grib_arguments * args)114 static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
115 {
116 int i, j, err=0;
117 double *lats, *lons; /* the lat/lon arrays to be populated */
118 long nx,ny,iScansNegatively,jScansPositively,jPointsAreConsecutive,alternativeRowScanning;
119 double LoVInDegrees,LaDInDegrees,Latin1InDegrees,Latin2InDegrees,latFirstInDegrees,
120 lonFirstInDegrees, Dx, Dy, radius=0;
121 double latFirstInRadians, lonFirstInRadians, LoVInRadians, Latin1InRadians, Latin2InRadians,
122 LaDInRadians, lonDiff, lonDeg, latDeg;
123 double f, n, rho, rho0, angle, x0, y0, x, y, tmp, tmp2;
124
125 grib_iterator_lambert_conformal* self = (grib_iterator_lambert_conformal*)iter;
126
127 const char* sradius = grib_arguments_get_name(h,args,self->carg++);
128 const char* snx = grib_arguments_get_name(h,args,self->carg++);
129 const char* sny = grib_arguments_get_name(h,args,self->carg++);
130 const char* sLoVInDegrees = grib_arguments_get_name(h,args,self->carg++);
131 const char* sLaDInDegrees = grib_arguments_get_name(h,args,self->carg++);
132 const char* sLatin1InDegrees = grib_arguments_get_name(h,args,self->carg++);
133 const char* sLatin2InDegrees = grib_arguments_get_name(h,args,self->carg++);
134 const char* slatFirstInDegrees = grib_arguments_get_name(h,args,self->carg++);
135 const char* slonFirstInDegrees = grib_arguments_get_name(h,args,self->carg++);
136 /* Dx and Dy are in Metres */
137 const char* sDx = grib_arguments_get_name(h,args,self->carg++);
138 const char* sDy = grib_arguments_get_name(h,args,self->carg++);
139 const char* siScansNegatively = grib_arguments_get_name(h,args,self->carg++);
140 const char* sjScansPositively = grib_arguments_get_name(h,args,self->carg++);
141 const char* sjPointsAreConsecutive = grib_arguments_get_name(h,args,self->carg++);
142 const char* salternativeRowScanning = grib_arguments_get_name(h,args,self->carg++);
143
144 if((err = grib_get_long_internal(h, snx,&nx)) != GRIB_SUCCESS)
145 return err;
146 if((err = grib_get_long_internal(h, sny,&ny)) != GRIB_SUCCESS)
147 return err;
148
149 if (iter->nv!=nx*ny) {
150 grib_context_log(h->context,GRIB_LOG_ERROR,"Wrong number of points (%ld!=%ldx%ld)",iter->nv,nx,ny);
151 return GRIB_WRONG_GRID;
152 }
153 if((err = grib_get_double_internal(h, sradius,&radius)) != GRIB_SUCCESS)
154 return err;
155 if((err = grib_get_double_internal(h, sLoVInDegrees, &LoVInDegrees))!=GRIB_SUCCESS)
156 return err;
157 if((err = grib_get_double_internal(h, sLaDInDegrees, &LaDInDegrees))!=GRIB_SUCCESS)
158 return err;
159 if((err = grib_get_double_internal(h, sLatin1InDegrees, &Latin1InDegrees))!=GRIB_SUCCESS)
160 return err;
161 if((err = grib_get_double_internal(h, sLatin2InDegrees, &Latin2InDegrees))!=GRIB_SUCCESS)
162 return err;
163 if((err = grib_get_double_internal(h, slatFirstInDegrees,&latFirstInDegrees)) !=GRIB_SUCCESS)
164 return err;
165 if((err = grib_get_double_internal(h, slonFirstInDegrees,&lonFirstInDegrees)) !=GRIB_SUCCESS)
166 return err;
167 if((err = grib_get_double_internal(h, sDx,&Dx)) !=GRIB_SUCCESS)
168 return err;
169 if((err = grib_get_double_internal(h, sDy,&Dy)) !=GRIB_SUCCESS)
170 return err;
171 if((err = grib_get_long_internal(h, sjPointsAreConsecutive,&jPointsAreConsecutive)) !=GRIB_SUCCESS)
172 return err;
173 if((err = grib_get_long_internal(h, sjScansPositively,&jScansPositively)) !=GRIB_SUCCESS)
174 return err;
175 if((err = grib_get_long_internal(h, siScansNegatively,&iScansNegatively)) !=GRIB_SUCCESS)
176 return err;
177 if((err = grib_get_long_internal(h, salternativeRowScanning,&alternativeRowScanning)) !=GRIB_SUCCESS)
178 return err;
179
180 /* See Wolfram MathWorld: http://mathworld.wolfram.com/LambertConformalConicProjection.html */
181 latFirstInRadians = latFirstInDegrees * DEG2RAD;
182 lonFirstInRadians = lonFirstInDegrees * DEG2RAD;
183 Latin1InRadians = Latin1InDegrees * DEG2RAD;
184 Latin2InRadians = Latin2InDegrees * DEG2RAD;
185 LaDInRadians = LaDInDegrees * DEG2RAD;
186 LoVInRadians = LoVInDegrees * DEG2RAD;
187
188 if ( fabs(Latin1InRadians - Latin2InRadians) < 1E-09 ) {
189 n = sin(Latin1InRadians);
190 }
191 else {
192 n = log(cos(Latin1InRadians)/cos(Latin2InRadians)) /
193 log(tan(M_PI_4 + Latin2InRadians/2.0) / tan(M_PI_4 + Latin1InRadians/2.0));
194 }
195
196 f = (cos(Latin1InRadians) * pow(tan(M_PI_4 + Latin1InRadians/2.0), n)) / n;
197 rho = radius * f * pow(tan(M_PI_4 + latFirstInRadians/2.0), -n);
198 rho0 = radius * f * pow(tan(M_PI_4 + LaDInRadians/2.0), -n);
199 if ( n < 0 ) /* adjustment for southern hemisphere */
200 rho0 = -rho0;
201 lonDiff = lonFirstInRadians - LoVInRadians;
202
203 /* Adjust longitude to range -180 to 180 */
204 if (lonDiff > M_PI) lonDiff -= 2*M_PI;
205 if (lonDiff < -M_PI) lonDiff += 2*M_PI;
206 angle = n * lonDiff;
207 x0 = rho * sin(angle);
208 y0 = rho0 - rho * cos(angle);
209 Dx = iScansNegatively == 0 ? Dx : -Dx;
210 /* GRIB-405: Don't change sign of Dy. Latitudes ALWAYS increase from latitudeOfFirstGridPoint */
211 /*Dy = jScansPositively == 1 ? Dy : -Dy;*/
212
213 /* No support (yet) for jPointsAreConsecutive */
214 if (jPointsAreConsecutive) {
215 grib_context_log(h->context,GRIB_LOG_ERROR,"No support for: 'Adjacent points in j (y) direction being consecutive'");
216 Assert(0);
217 }
218
219 /* Allocate latitude and longitude arrays */
220 self->lats = (double*)grib_context_malloc(h->context,iter->nv*sizeof(double));
221 if (!self->lats) {
222 grib_context_log(h->context,GRIB_LOG_ERROR, "Unable to allocate %ld bytes",iter->nv*sizeof(double));
223 return GRIB_OUT_OF_MEMORY;
224 }
225 self->lons = (double*)grib_context_malloc(h->context,iter->nv*sizeof(double));
226 if (!self->lats) {
227 grib_context_log(h->context,GRIB_LOG_ERROR, "Unable to allocate %ld bytes",iter->nv*sizeof(double));
228 return GRIB_OUT_OF_MEMORY;
229 }
230 lats=self->lats;
231 lons=self->lons;
232
233 /* Populate our arrays */
234 for (j = 0; j < ny; j++) {
235 y = y0 + j*Dy;
236 if ( n < 0 ) { /* adjustment for southern hemisphere */
237 y = -y;
238 }
239 tmp = rho0 - y;
240 tmp2 = tmp*tmp;
241 for (i = 0; i < nx; i++) {
242 int index =i+j*nx;
243 x = x0 + i*Dx;
244 if ( n < 0 ) { /* adjustment for southern hemisphere */
245 x = -x;
246 }
247
248 angle = atan2(x, tmp);
249 rho = sqrt(x*x + tmp2);
250 if (n <= 0) rho = -rho;
251 lonDeg = LoVInDegrees + (angle/n) * RAD2DEG;
252 latDeg = (2.0 * atan(pow(radius * f/rho, 1.0/n)) - M_PI_2) * RAD2DEG;
253 while ( lonDeg >= 360.0) lonDeg -= 360.0;
254 while ( lonDeg < 0.0) lonDeg += 360.0;
255 lons[index] = lonDeg;
256 lats[index] = latDeg;
257 }
258 }
259
260 iter->e = -1;
261
262 return err;
263 }
264
destroy(grib_iterator * i)265 static int destroy(grib_iterator* i)
266 {
267 grib_iterator_lambert_conformal* self = (grib_iterator_lambert_conformal*)i;
268 const grib_context *c = i->h->context;
269
270 grib_context_free(c,self->lats);
271 grib_context_free(c,self->lons);
272 return 1;
273 }
274