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 = destroy
21    IMPLEMENTS = init;next
22    MEMBERS     =   double *lats
23    MEMBERS     =   double *lons
24    MEMBERS     =   long Nj
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_polar_stereographic
47 {
48     grib_iterator it;
49     /* Members defined in gen */
50     long carg;
51     const char* missingValue;
52     /* Members defined in polar_stereographic */
53     double* lats;
54     double* lons;
55     long Nj;
56 } grib_iterator_polar_stereographic;
57 
58 extern grib_iterator_class* grib_iterator_class_gen;
59 
60 static grib_iterator_class _grib_iterator_class_polar_stereographic = {
61     &grib_iterator_class_gen,                  /* super                     */
62     "polar_stereographic",                     /* name                      */
63     sizeof(grib_iterator_polar_stereographic), /* size of instance          */
64     0,                                         /* inited */
65     &init_class,                               /* init_class */
66     &init,                                     /* constructor               */
67     &destroy,                                  /* destructor                */
68     &next,                                     /* Next Value                */
69     0,                                         /*  Previous Value           */
70     0,                                         /* Reset the counter         */
71     0,                                         /* has next values           */
72 };
73 
74 grib_iterator_class* grib_iterator_class_polar_stereographic = &_grib_iterator_class_polar_stereographic;
75 
76 
init_class(grib_iterator_class * c)77 static void init_class(grib_iterator_class* c)
78 {
79     c->previous = (*(c->super))->previous;
80     c->reset    = (*(c->super))->reset;
81     c->has_next = (*(c->super))->has_next;
82 }
83 /* END_CLASS_IMP */
84 
next(grib_iterator * i,double * lat,double * lon,double * val)85 static int next(grib_iterator* i, double* lat, double* lon, double* val)
86 {
87     grib_iterator_polar_stereographic* self = (grib_iterator_polar_stereographic*)i;
88 
89     if ((long)i->e >= (long)(i->nv - 1))
90         return 0;
91     i->e++;
92 
93     *lat = self->lats[i->e];
94     *lon = self->lons[i->e];
95     *val = i->data[i->e];
96 
97     return 1;
98 }
99 
100 /* Data struct for Forward and Inverse Projections */
101 typedef struct proj_data_t
102 {
103     double centre_lon;     /* central longitude */
104     double centre_lat;     /* central latitude */
105     double sign;           /* sign variable */
106     double ind;            /* flag variable */
107     double mcs;            /* small m */
108     double tcs;            /* small t */
109     double false_northing; /* y offset in meters */
110     double false_easting;  /* x offset in meters */
111 } proj_data_t;
112 
113 #define RAD2DEG 57.29577951308232087684 /* 180 over pi */
114 #define DEG2RAD 0.01745329251994329576  /* pi over 180 */
115 #define PI_OVER_2 1.5707963267948966    /* half pi */
116 #define EPSILON 1.0e-10
117 
init(grib_iterator * iter,grib_handle * h,grib_arguments * args)118 static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
119 {
120     int ret = 0;
121     double *lats, *lons; /* arrays for latitudes and longitudes */
122     double lonFirstInDegrees, latFirstInDegrees, radius;
123     double x, y, Dx, Dy;
124     long nx, ny, centralLongitudeInDegrees, centralLatitudeInDegrees;
125     long alternativeRowScanning, iScansNegatively, i, j;
126     long jScansPositively, jPointsAreConsecutive, southPoleOnPlane;
127     double centralLongitude, centralLatitude; /* in radians */
128     double con1;                              /* temporary angle */
129     double ts;                                /* value of small t */
130     double height;                            /* height above ellipsoid */
131     double x0, y0, lonFirst, latFirst;
132     proj_data_t fwd_proj_data = {0,};
133     proj_data_t inv_proj_data = {0,};
134 
135     grib_iterator_polar_stereographic* self = (grib_iterator_polar_stereographic*)iter;
136 
137     const char* sradius                 = grib_arguments_get_name(h, args, self->carg++);
138     const char* snx                     = grib_arguments_get_name(h, args, self->carg++);
139     const char* sny                     = grib_arguments_get_name(h, args, self->carg++);
140     const char* slatFirstInDegrees      = grib_arguments_get_name(h, args, self->carg++);
141     const char* slonFirstInDegrees      = grib_arguments_get_name(h, args, self->carg++);
142     const char* ssouthPoleOnPlane       = grib_arguments_get_name(h, args, self->carg++);
143     const char* scentralLongitude       = grib_arguments_get_name(h, args, self->carg++);
144     const char* scentralLatitude        = grib_arguments_get_name(h, args, self->carg++);
145     const char* sDx                     = grib_arguments_get_name(h, args, self->carg++);
146     const char* sDy                     = grib_arguments_get_name(h, args, self->carg++);
147     const char* siScansNegatively       = grib_arguments_get_name(h, args, self->carg++);
148     const char* sjScansPositively       = grib_arguments_get_name(h, args, self->carg++);
149     const char* sjPointsAreConsecutive  = grib_arguments_get_name(h, args, self->carg++);
150     const char* salternativeRowScanning = grib_arguments_get_name(h, args, self->carg++);
151 
152     if (grib_is_earth_oblate(h)) {
153         grib_context_log(h->context, GRIB_LOG_ERROR, "Polar stereographic only supported for spherical earth.");
154         return GRIB_GEOCALCULUS_PROBLEM;
155     }
156 
157     if ((ret = grib_get_double_internal(h, sradius, &radius)) != GRIB_SUCCESS)
158         return ret;
159     if ((ret = grib_get_long_internal(h, snx, &nx)) != GRIB_SUCCESS)
160         return ret;
161     if ((ret = grib_get_long_internal(h, sny, &ny)) != GRIB_SUCCESS)
162         return ret;
163 
164     if (iter->nv != nx * ny) {
165         grib_context_log(h->context, GRIB_LOG_ERROR, "Wrong number of points (%ld!=%ldx%ld)", iter->nv, nx, ny);
166         return GRIB_WRONG_GRID;
167     }
168     if ((ret = grib_get_double_internal(h, slatFirstInDegrees, &latFirstInDegrees)) != GRIB_SUCCESS)
169         return ret;
170     if ((ret = grib_get_double_internal(h, slonFirstInDegrees, &lonFirstInDegrees)) != GRIB_SUCCESS)
171         return ret;
172     if ((ret = grib_get_long_internal(h, ssouthPoleOnPlane, &southPoleOnPlane)) != GRIB_SUCCESS)
173         return ret;
174     if ((ret = grib_get_long_internal(h, scentralLongitude, &centralLongitudeInDegrees)) != GRIB_SUCCESS)
175         return ret;
176     if ((ret = grib_get_long_internal(h, scentralLatitude, &centralLatitudeInDegrees)) != GRIB_SUCCESS)
177         return ret;
178     if ((ret = grib_get_double_internal(h, sDx, &Dx)) != GRIB_SUCCESS)
179         return ret;
180     if ((ret = grib_get_double_internal(h, sDy, &Dy)) != GRIB_SUCCESS)
181         return ret;
182     if ((ret = grib_get_long_internal(h, sjPointsAreConsecutive, &jPointsAreConsecutive)) != GRIB_SUCCESS)
183         return ret;
184     if ((ret = grib_get_long_internal(h, sjScansPositively, &jScansPositively)) != GRIB_SUCCESS)
185         return ret;
186     if ((ret = grib_get_long_internal(h, siScansNegatively, &iScansNegatively)) != GRIB_SUCCESS)
187         return ret;
188     if ((ret = grib_get_long_internal(h, salternativeRowScanning, &alternativeRowScanning)) != GRIB_SUCCESS)
189         return ret;
190 
191     centralLongitude = centralLongitudeInDegrees * DEG2RAD;
192     centralLatitude  = centralLatitudeInDegrees * DEG2RAD;
193     lonFirst         = lonFirstInDegrees * DEG2RAD;
194     latFirst         = latFirstInDegrees * DEG2RAD;
195 
196     /* Forward projection initialisation */
197     fwd_proj_data.false_northing = 0;
198     fwd_proj_data.false_easting  = 0;
199     fwd_proj_data.centre_lon     = centralLongitude;
200     fwd_proj_data.centre_lat     = centralLatitude;
201     if (centralLatitude < 0)
202         fwd_proj_data.sign = -1.0;
203     else
204         fwd_proj_data.sign = +1.0;
205     fwd_proj_data.ind = 0;
206     if (fabs(fabs(centralLatitude) - PI_OVER_2) > EPSILON) {
207         /* central latitude different from 90 i.e. not north/south polar */
208         fwd_proj_data.ind = 1;
209         con1              = fwd_proj_data.sign * centralLatitude;
210         fwd_proj_data.mcs = cos(con1);
211         fwd_proj_data.tcs = tan(0.5 * (PI_OVER_2 - con1));
212     }
213 
214     /* Forward projection from initial lat,lon to initial x,y */
215     con1 = fwd_proj_data.sign * (lonFirst - fwd_proj_data.centre_lon);
216     ts   = tan(0.5 * (PI_OVER_2 - fwd_proj_data.sign * latFirst));
217     if (fwd_proj_data.ind)
218         height = radius * fwd_proj_data.mcs * ts / fwd_proj_data.tcs;
219     else
220         height = 2.0 * radius * ts;
221     x0 = fwd_proj_data.sign * height * sin(con1) + fwd_proj_data.false_easting;
222     y0 = -fwd_proj_data.sign * height * cos(con1) + fwd_proj_data.false_northing;
223 
224     x0 = -x0;
225     y0 = -y0;
226 
227     /* Inverse projection initialisation */
228     inv_proj_data.false_easting  = x0;
229     inv_proj_data.false_northing = y0;
230     inv_proj_data.centre_lon     = centralLongitude;
231     inv_proj_data.centre_lat     = centralLatitude;
232     if (centralLatitude < 0)
233         inv_proj_data.sign = -1.0;
234     else
235         inv_proj_data.sign = +1.0;
236     inv_proj_data.ind = 0;
237     if (fabs(fabs(centralLatitude) - PI_OVER_2) > EPSILON) {
238         inv_proj_data.ind = 1;
239         con1              = inv_proj_data.sign * inv_proj_data.centre_lat;
240         inv_proj_data.mcs = cos(con1);
241         inv_proj_data.tcs = tan(0.5 * (PI_OVER_2 - con1));
242     }
243     self->lats = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
244     if (!self->lats) {
245         grib_context_log(h->context, GRIB_LOG_ERROR, "Error allocating %ld bytes", iter->nv * sizeof(double));
246         return GRIB_OUT_OF_MEMORY;
247     }
248     self->lons = (double*)grib_context_malloc(h->context, iter->nv * sizeof(double));
249     if (!self->lats) {
250         grib_context_log(h->context, GRIB_LOG_ERROR, "Error allocating %ld bytes", iter->nv * sizeof(double));
251         return GRIB_OUT_OF_MEMORY;
252     }
253     lats = self->lats;
254     lons = self->lons;
255     Dx   = iScansNegatively == 0 ? Dx : -Dx;
256     Dy   = jScansPositively == 1 ? Dy : -Dy;
257 
258     y = 0;
259     for (j = 0; j < ny; j++) {
260         x = 0;
261         for (i = 0; i < nx; i++) {
262             /* Inverse projection from x,y to lat,lon */
263             /* int index =i+j*nx; */
264             double _x = (x - inv_proj_data.false_easting) * inv_proj_data.sign;
265             double _y = (y - inv_proj_data.false_northing) * inv_proj_data.sign;
266             double rh = sqrt(_x * _x + _y * _y);
267             if (inv_proj_data.ind)
268                 ts = rh * inv_proj_data.tcs / (radius * inv_proj_data.mcs);
269             else
270                 ts = rh / (radius * 2.0);
271             *lats = inv_proj_data.sign * (PI_OVER_2 - 2 * atan(ts));
272             if (rh == 0) {
273                 *lons = inv_proj_data.sign * inv_proj_data.centre_lon;
274             }
275             else {
276                 double temp = atan2(_x, -_y);
277                 *lons       = inv_proj_data.sign * temp + inv_proj_data.centre_lon;
278             }
279             *lats = *lats * RAD2DEG;
280             *lons = *lons * RAD2DEG;
281             while (*lons < 0)
282                 *lons += 360;
283             while (*lons > 360)
284                 *lons -= 360;
285             lons++;
286             lats++;
287 
288             x += Dx;
289         }
290         y += Dy;
291     }
292 #if 0
293     /*standardParallel = (southPoleOnPlane == 1) ? -90 : +90;*/
294     if (jPointsAreConsecutive)
295     {
296         x=xFirst;
297         for (i=0;i<nx;i++) {
298             y=yFirst;
299             for (j=0;j<ny;j++) {
300                 rho=sqrt(x*x+y*y);
301                 if (rho == 0) {
302                     /* indeterminate case */
303                     *lats = standardParallel;
304                     *lons = centralLongitude;
305                 }
306                 else {
307                     c=2*atan2(rho,(2.0*radius));
308                     cosc=cos(c);
309                     sinc=sin(c);
310                     *lats = asin( cosc*sinphi1 + y*sinc*cosphi1/rho ) * RAD2DEG;
311                     *lons = (lambda0+atan2(x*sinc, rho*cosphi1*cosc - y*sinphi1*sinc)) * RAD2DEG;
312                 }
313                 while (*lons<0)   *lons += 360;
314                 while (*lons>360) *lons -= 360;
315                 lons++;
316                 lats++;
317 
318                 y+=Dy;
319             }
320             x+=Dx;
321         }
322     }
323     else
324     {
325         y=yFirst;
326         for (j=0;j<ny;j++) {
327             x=xFirst;
328             for (i=0;i<nx;i++) {
329                 /* int index =i+j*nx; */
330                 rho=sqrt(x*x+y*y);
331                 if (rho == 0) {
332                     /* indeterminate case */
333                     *lats = standardParallel;
334                     *lons = centralLongitude;
335                 }
336                 else {
337                     c=2*atan2(rho,(2.0*radius));
338                     cosc=cos(c);
339                     sinc=sin(c);
340                     *lats = asin( cosc*sinphi1 + y*sinc*cosphi1/rho ) * RAD2DEG;
341                     *lons = (lambda0+atan2(x*sinc, rho*cosphi1*cosc - y*sinphi1*sinc)) * RAD2DEG;
342                 }
343                 while (*lons<0)   *lons += 360;
344                 while (*lons>360) *lons -= 360;
345                 lons++;
346                 lats++;
347 
348                 x+=Dx;
349             }
350             y+=Dy;
351         }
352     }
353 #endif
354     iter->e = -1;
355 
356     return ret;
357 }
358 
destroy(grib_iterator * i)359 static int destroy(grib_iterator* i)
360 {
361     grib_iterator_polar_stereographic* self = (grib_iterator_polar_stereographic*)i;
362     const grib_context* c                   = i->h->context;
363 
364     grib_context_free(c, self->lats);
365     grib_context_free(c, self->lons);
366     return 1;
367 }
368