← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to imagery/i.atcorr/geomcond.cpp

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
imagery/i.atcorr/geomcond.cpp
 
1
extern "C" {
 
2
#include <grass/gis.h>
 
3
#include <grass/glocale.h>
 
4
}
 
5
 
 
6
#include "geomcond.h"
 
7
#include "common.h"
 
8
 
 
9
/* **********************************************************************c */
 
10
/*                                                                      c */
 
11
/*                                                *     sun             c */
 
12
/*                                              \ * /                   c */
 
13
/*                                            * * * * *                 c */
 
14
/*                                   z          / * \                   c */
 
15
/*                                   +           /+                     c */
 
16
/*            satellite    /         +          /                       c */
 
17
/*                       o/          +         /                        c */
 
18
/*                      /.\          +        /.                        c */
 
19
/*                     / . \  _avis-_+_-asol_/ .                        c */
 
20
/*                       .  \-      -+      /  .    north               c */
 
21
/*                       .   \       +     /   .  +                     c */
 
22
/*                       .    \      +    /    .+                       c */
 
23
/*                       .     \     +   /    +.                        c */
 
24
/*                       .      \    +  /   +  .                        c */
 
25
/*                       .       \   + /  +    .                        c */
 
26
/*                       .        \  +/ +      .                        c */
 
27
/*    west + + + + + + + . + + + + +\+ + + + + . + + + + + + + + east   c */
 
28
/*                       .          +..        .                        c */
 
29
/*                       .        + .   .      .                        c */
 
30
/*                       .      +  .      .    .                        c */
 
31
/*                       .    +   .       .'.  .                        c */
 
32
/*                       .  +    .. . , '     ..                        c */
 
33
/*                       .+     .       \       .                       c */
 
34
/*                      +.     .         \        .                     c */
 
35
/*                    +  .    .           \         .                   c */
 
36
/*             south     .   .       (phiv-phi0)                        c */
 
37
/*                                                                      c */
 
38
/*                                                                      c */
 
39
/*                                                                      c */
 
40
/* **********************************************************************c */
 
41
 
 
42
 
 
43
/*      To take into account the variation of the solar constant as a function 
 
44
  of the Julian day. 
 
45
 
 
46
  return dsol           
 
47
  dsol is a multiplicative factor to apply to the mean value of solar constant 
 
48
*/
 
49
float GeomCond::varsol ()
 
50
{
 
51
/* calculation of the variability of the solar constant during the year. 
 
52
   jday is the number of the day in the month   */
 
53
    long int j;
 
54
    if (month <= 2) j = (month - 1) * 31 + jday;
 
55
    else if (month > 8) j = (month - 1) * 31 - (month - 2) / 2 - 2 + jday;
 
56
    else j = (month - 1) * 31 - (month - 1) / 2 - 2 + jday;
 
57
 
 
58
/* Computing 2nd power */
 
59
    double tmp = 1.f - cos ((float) (j - 4) * 0.9856f * M_PI / 180.f) * .01673f;
 
60
    return 1.f / (float)(tmp * tmp);
 
61
}
 
62
 
 
63
 
 
64
/* spot, landsat5 and landsat7 is handled the same way */
 
65
void GeomCond::landsat(float tu)
 
66
{
 
67
/*     warning !!! */
 
68
/*     xlon and xlat are the coordinates of the scene center. */
 
69
    avis = 0.f;
 
70
    phiv = 0.f;
 
71
    possol(tu);
 
72
}
 
73
 
 
74
/*
 
75
  To compute the solar azimuthal and zenithal angles (in degrees) for a point over
 
76
  the globe defined by its longitude and its latitude (in dec. degrees) for a day of the year (fixed by
 
77
  number of the month and number of the day in the month) at any Greenwich Meridian Time (GMT
 
78
  dec. hour).
 
79
*/
 
80
void GeomCond::possol(float tu)
 
81
{
 
82
    long int ia = 0;
 
83
    long int nojour;
 
84
/*     solar position (zenithal angle asol,azimuthal angle phi0 */
 
85
/*                     in degrees) */
 
86
/*     jday is the number of the day in the month */
 
87
    day_number(ia, nojour);
 
88
    pos_fft (nojour, tu);
 
89
    if (asol > 90.f)
 
90
        G_warning(_("The sun is not raised"));
 
91
}
 
92
 
 
93
void GeomCond::day_number(long int ia, long int& j)
 
94
{
 
95
    if (month <= 2)
 
96
    {
 
97
        j = (month - 1) * 31 + jday;
 
98
        return;
 
99
    }
 
100
 
 
101
    if (month > 8) j = (month - 1) * 31 - (month - 2) / 2 - 2 + jday;
 
102
    else j = (month - 1) * 31 - (month - 1) / 2 - 2 + jday;
 
103
 
 
104
    if (ia != 0 && ia % 4 == 0) ++j;
 
105
}
 
106
 
 
107
/* returns the sign of the element */
 
108
#define SIGN(X) (((X) >= 0) ? 1. : -1.) 
 
109
 
 
110
void GeomCond::pos_fft (long int j, float tu)
 
111
{
 
112
    /* Local variables */
 
113
    double ah, et, az, caz, xla, tet, tsm, tsv, elev, azim, delta, amuzero;
 
114
 
 
115
    /*     solar position (zenithal angle asol,azimuthal angle phi0 */
 
116
    /*                     in degrees) */
 
117
    /*     j is the day number in the year */
 
118
 
 
119
    /* mean solar time (heure decimale) */
 
120
    tsm = tu + xlon / 15.;
 
121
    xla = xlat * M_PI / 180.;
 
122
    tet = (float)(j) * M_PI2 / 365.;
 
123
 
 
124
    /* time equation (in mn.dec) */
 
125
    et = 7.5e-5f + 0.001868f * cos (tet) - 0.032077f * sin (tet) - 
 
126
        0.014615f * cos (tet * 2.f) - 0.040849f * sin (tet * 2.f);
 
127
 
 
128
    et = et * 12.f * 60.f / M_PI;
 
129
 
 
130
    /* true solar time */
 
131
    tsv = tsm + et / 60.f;
 
132
    tsv += -12.f;
 
133
 
 
134
    /* hour angle */
 
135
    ah = tsv * 15.f * M_PI / 180.f;
 
136
 
 
137
    /* solar declination   (in radian) */
 
138
    delta = 0.006918f - 0.399912f * cos (tet) + 0.070257f * sin (tet) - 
 
139
        0.006758f * cos (tet * 2.f) + 9.07e-4f * sin (tet * 2.f) - 
 
140
        0.002697f * cos (tet * 3.f) + 0.00148f * sin (tet * 3.f);
 
141
 
 
142
    /* elevation,azimuth */
 
143
    amuzero = sin (xla) * sin (delta) + cos (xla) * cos (delta) * cos (ah);
 
144
    elev = asin (amuzero);
 
145
    az = cos (delta) * sin (ah) / cos (elev);
 
146
  
 
147
    if (fabs (az) - 1.f > 0.f) az = SIGN(az);
 
148
 
 
149
    caz = (-cos (xla) * sin (delta) + sin (xla) * cos (delta) * cos (ah)) / cos (elev);
 
150
    azim = asin (az);
 
151
    if (caz <= 0.f) azim = M_PI - azim;
 
152
 
 
153
    if (caz > 0.f && az <= 0.f) azim += M_PI2;
 
154
 
 
155
    azim += M_PI;
 
156
    if (azim > M_PI2) azim -= M_PI2;
 
157
        
 
158
    elev = elev * 180. / M_PI;
 
159
        
 
160
    /*     conversion in degrees */
 
161
    asol = (float)(90. - elev);
 
162
    phi0 = (float)(azim * 180. / M_PI);
 
163
}
 
164
 
 
165
/*
 
166
  convert:
 
167
  1 = meteosat observation 
 
168
  2 = goes east observation
 
169
  3 = goes west observation
 
170
*/
 
171
void GeomCond::posobs(float tu, int nc, int nl)
 
172
{
 
173
    double yr, xr, alti;
 
174
 
 
175
    if(igeom == 1) /* meteosat observation */
 
176
    {
 
177
        yr = nl - 1250.5;
 
178
        xr = nc - 2500.5;
 
179
        alti = 42164.0 - 6378.155;
 
180
    } 
 
181
    else if(igeom == 2) /* goes east observation */
 
182
    {
 
183
        yr = nl - 8665.5;
 
184
        xr = nc - 6498.5;
 
185
        alti = 42107.0 - 6378.155;
 
186
    }
 
187
    else /* goes west observation */
 
188
    {
 
189
        yr = nl - 8665.5;
 
190
        xr = nc - 6498.5;
 
191
        alti = 42147.0 - 6378.155;
 
192
    }
 
193
 
 
194
 
 
195
    const double re = 6378.155;
 
196
    const double aaa = 1. / 297.;
 
197
    const double rp = re / (1.f + aaa);
 
198
    const double cdr = M_PI / 180.;
 
199
    const double crd = 180. / M_PI;
 
200
 
 
201
    double deltax;
 
202
    double deltay;
 
203
 
 
204
    if(igeom == 1) 
 
205
    {
 
206
        deltax = 18.0 / 5000.0;
 
207
        deltay = 18.0 / 2500.0;
 
208
    }
 
209
    else
 
210
    {
 
211
        deltax = 18.0 / 12997.0;
 
212
        deltay = 20.0 / 17331.0;
 
213
    }
 
214
 
 
215
    double x = xr * deltax * cdr;
 
216
    double y = yr * deltay * cdr;
 
217
    double rs = re + alti;
 
218
    double tanx = tan(x);
 
219
    double tany = tan(y);
 
220
    double val1 = 1.0 + (tanx * tanx);
 
221
    double val2 = 1.0 + (tany * (1.0 + aaa)) * (tany * (1.0 + aaa));
 
222
    double yk = rs / re;
 
223
    double cosx2 = 1. / (val1 * val2);
 
224
      
 
225
    double sn, zt, xt, yt, teta, ylat, ylon;
 
226
    if((1. / cosx2) > ((yk * yk) / (yk*yk - 1.)))
 
227
    {
 
228
        G_warning(_("No possibility to compute lat. and long."));
 
229
        return;
 
230
    }
 
231
    else
 
232
    {
 
233
        sn = (rs - (re * (sqrt((yk * yk) - (yk*yk - 1.) * (1. / cosx2))))) / (1. / cosx2);
 
234
        zt = rs - sn;
 
235
        xt = -(sn * tanx);
 
236
        yt = sn * tany / cos(x);
 
237
        teta = asin(yt / rp);
 
238
        ylat = (atan(((tan(teta)) * rp) / re));
 
239
        ylon = atan(xt / zt);
 
240
    }
 
241
 
 
242
    xlat = (float)(ylat * crd);
 
243
 
 
244
    if(igeom == 1) xlon = (float)(ylon * crd);
 
245
    else if(igeom == 2) xlon = (float)(ylon * crd - 75.);
 
246
    else xlon = (float)(ylon * crd - 135.);
 
247
 
 
248
    possol(tu);
 
249
 
 
250
    if(igeom == 1) ylon = xlon * M_PI / 180.;
 
251
    else if(igeom == 2) ylon = xlon * M_PI / 180. + 75. * cdr;
 
252
    else ylon = xlon * M_PI / 180. + 135. * cdr;
 
253
 
 
254
    ylat = xlat * M_PI / 180.;
 
255
    double gam = sqrt(((1. / cosx2) - 1.) * cosx2);
 
256
    avis = (float)(asin((1. + alti / re) * (gam)) * 180. / M_PI);
 
257
    phiv = (float)((atan2(tan(ylon),sin(ylat)) + M_PI) * 180. / M_PI);
 
258
}
 
259
 
 
260
void GeomCond::posnoa(float tu, int nc, float xlonan, float campm, float hna)
 
261
{
 
262
/*     noaa 6 definition
 
263
       orbite inclination ai in radians
 
264
       hor mouv in rad/s  an
 
265
       h/r=860/6378
 
266
       campm allows the user to switch to pm platforms */
 
267
 
 
268
    const double r = 860. / 6378.155;
 
269
    const double ai = 98.96 * M_PI / 180.;
 
270
    const double an = 360. * M_PI / (6119. * 180.);
 
271
    double ylonan = xlonan * M_PI / 180.;
 
272
    double t = tu * 3600;
 
273
    double hnam = hna;
 
274
    hnam = hnam * 3600;
 
275
    double u = t - hnam;
 
276
    u = campm * u * an;
 
277
    double delt = ((nc - (2048 + 1) / 2.) * 55.385 / ((2048. - 1) / 2.));
 
278
    delt = campm * delt * M_PI / 180.;
 
279
    avis = (float)asin((1 + r) * sin(delt));
 
280
    double d = avis - delt;
 
281
    double y = cos(d) * cos(ai) * sin(u) - sin(ai) * sin(d);
 
282
    double z = cos(d) * sin(ai) * sin(u) + cos(ai) * sin(d);
 
283
    double ylat = asin(z);
 
284
    double cosy = cos(d) * cos(u) / cos(ylat);
 
285
 
 
286
    double siny = y / cos(ylat);
 
287
    double ylon = asin(siny);
 
288
    if(cosy <= 0.)
 
289
    {
 
290
        if(siny > 0) ylon = M_PI - ylon;
 
291
        if(siny <= 0) ylon = -(M_PI + ylon);
 
292
    }
 
293
    double ylo1 = ylon + ylonan - (t - hnam) * 2. * M_PI / 86400.;
 
294
    xlat = (float)(ylat * 180. / M_PI);
 
295
    xlon = (float)(ylo1 * 180. / M_PI);
 
296
 
 
297
 
 
298
 
 
299
    possol(tu);
 
300
 
 
301
    double zlat = asin(sin(ai) * sin(u));
 
302
    double zlon = atan2(cos(ai) * sin(u),cos(u));
 
303
    if(nc != 1024)
 
304
    {
 
305
        double xnum = sin(zlon - ylon) * cos(zlat) / sin(fabs(d));
 
306
        double xden = (sin(zlat) - sin(ylat) * cos(d)) / cos(ylat) / sin(fabs(d));
 
307
        phiv = (float)atan2(xnum,xden);
 
308
    }
 
309
    else phiv = 0.;
 
310
    phiv = (float)(phiv * 180. / M_PI);
 
311
    avis = (float)(fabs(avis) * 180. / M_PI);
 
312
}
 
313
 
 
314
void GeomCond::parse()
 
315
{
 
316
    cin >> igeom;
 
317
    cin.ignore(numeric_limits<int>::max(),'\n');  /* read the rest of the scraps, like comments */
 
318
 
 
319
    float campm = -1.0f;        /* initialize in case igeom == 5 */
 
320
    float tu, xlonan, hna;
 
321
    int nc, nl;
 
322
 
 
323
    switch(igeom)
 
324
    {
 
325
    case 0: /* internal format */
 
326
    {
 
327
        cin >> asol;
 
328
        cin >> phi0;
 
329
        cin >> avis;
 
330
        cin >> phiv;
 
331
        cin >> month;
 
332
        cin >> jday;
 
333
        cin.ignore(numeric_limits<int>::max(),'\n');  /* read the rest of the scraps, like comments */
 
334
        break;
 
335
    }
 
336
    case 1: /* meteosat observation */
 
337
    case 2: /* goes east observation */
 
338
    case 3: /* goes west observation  */
 
339
    {
 
340
        cin >> month;
 
341
        cin >> jday;
 
342
        cin >> tu;
 
343
        cin >> nc;
 
344
        cin >> nl;
 
345
        cin.ignore(numeric_limits<int>::max(),'\n');
 
346
        posobs(tu, nc, nl);
 
347
        break;
 
348
    }
 
349
    case 4: campm = 1.0f;
 
350
    case 5: 
 
351
    {
 
352
        cin >> month;
 
353
        cin >> jday;
 
354
        cin >> tu;
 
355
        cin >> nc;
 
356
        cin >> xlonan;
 
357
        cin >> hna;
 
358
        cin.ignore(numeric_limits<int>::max(),'\n');
 
359
        posnoa(tu, nc, xlonan, campm, hna);
 
360
        break;
 
361
    }
 
362
    case 6: /* hrv   ( spot )    * enter month,day,hh.ddd,long.,lat. */
 
363
    case 7: /* tm    ( landsat ) * enter month,day,hh.ddd,long.,lat. */
 
364
    case 8: /* etm+  ( landsat7) * enter month,day,hh.ddd,long.,lat. */
 
365
    case 9: /* liss  ( IRS 1C)   * enter month,day,hh.ddd,long.,lat. */
 
366
    case 10: /* aster            * enter month,day,hh.ddd,long.,lat. */
 
367
    case 11: /* avnir            * enter month,day,hh.ddd,long.,lat. */
 
368
    case 12: /* ikonos           * enter month,day,hh.ddd,long.,lat. */
 
369
    case 13: /* rapideye         * enter month,day,hh.ddd,long.,lat. */
 
370
    case 14: /* vgt1_spot4       * enter month,day,hh.ddd,long.,lat. */
 
371
    case 15: /* vgt2_spot5       * enter month,day,hh.ddd,long.,lat. */
 
372
        case 16: /* worldview2       * enter month,day,hh.ddd,long.,lat. */
 
373
        case 17: /* quickbird2       * enter month,day,hh.ddd,long.,lat. */
 
374
    case 18: /* Landsat 8        * enter month,day,hh.ddd,long.,lat. */
 
375
    {
 
376
        cin >> month;
 
377
        cin >> jday;
 
378
        cin >> tu;
 
379
        cin >> xlon;
 
380
        cin >> xlat;
 
381
        cin.ignore(numeric_limits<int>::max(),'\n');  /* read the rest of the scraps, like comments */
 
382
        landsat(tu);
 
383
        break;
 
384
    }
 
385
    default: G_fatal_error(_("Unsupported/unreadable format in control file (found igeom=%ld)"), igeom);
 
386
    }
 
387
 
 
388
 
 
389
    /* ********************************************************************** */
 
390
    /*                                                                        */
 
391
    /*                                 / scattered direction                  */
 
392
    /*                               /                                        */
 
393
    /*                             /                                          */
 
394
    /*                           / adif                                       */
 
395
    /*    incident   + + + + + + + + + + + + + + +                            */
 
396
    /*    direction                                                           */
 
397
    /*                                                                        */
 
398
    /* ********************************************************************** */
 
399
    phi = (float)fabs(phiv - phi0);
 
400
    phirad = (phi0 - phiv) * (float)M_PI / 180.f;
 
401
    if (phirad < 0.f) phirad += (float)M_PI2;
 
402
    if (phirad > M_PI2) phirad -= (float)M_PI2;
 
403
 
 
404
    xmus = (float)cos (asol * M_PI / 180.f);
 
405
    xmuv = (float)cos (avis * M_PI / 180.f);
 
406
    xmup = (float)cos (phirad);
 
407
    xmud = -xmus * xmuv - (float)sqrt (1.f - xmus * xmus) * (float)sqrt (1.f - xmuv * xmuv) * xmup;
 
408
 
 
409
    /* test vermote bug */
 
410
    if (xmud > 1.f)  xmud = 1.f;
 
411
    if (xmud < -1.f) xmud = -1.f;
 
412
    adif = (float)acos (xmud) * 180.f / (float)M_PI;
 
413
 
 
414
    dsol = varsol();
 
415
}
 
416
 
 
417
/* ---- print geometrical conditions ---- */
 
418
void GeomCond::print()
 
419
{
 
420
    static const string etiq1[] = {
 
421
        string(" user defined conditions     "),
 
422
        string(" Meteosat observation        "),
 
423
        string(" GOES east observation       "),
 
424
        string(" GOES west observation       "),
 
425
        string(" AVHRR (AM noaa) observation "),
 
426
        string(" AVHRR (PM noaa) observation "),
 
427
        string(" H.R.V.   observation        "),
 
428
        string(" T.M.     observation        "),
 
429
        string(" ETM+     observation        "),
 
430
        string(" LISS     observation        "),
 
431
        string(" ASTER    observation        "),
 
432
        string(" AVNIR    observation        "),
 
433
        string(" Ikonos   observation        "),
 
434
        string(" Rapideye observation        "),
 
435
        string(" VGT1-SPOT4 observation      "),
 
436
        string(" VGT2-SPOT5 observation      "),
 
437
        string(" Worldview2 observation      "),
 
438
        string(" Quickbird2 observation      "),
 
439
    string(" Landsat 8 observation       ")
 
440
        };
 
441
 
 
442
    static const string head(" geometrical conditions identity  ");
 
443
    static const string line(" -------------------------------  ");
 
444
    Output::Begin(); Output::Repeat(22,' '); Output::Print(head); Output::End();
 
445
    Output::Begin(); Output::Repeat(22,' '); Output::Print(line); Output::End();
 
446
 
 
447
        
 
448
    Output::Begin(); Output::Repeat(22,' '); Output::Print(etiq1[igeom]); Output::End();
 
449
    Output::Begin(); Output::End();
 
450
 
 
451
        
 
452
    Output::Begin(); Output::Repeat(2,' ');
 
453
    ostringstream s1;
 
454
    s1.setf(ios::fixed, ios::floatfield);
 
455
    s1 << " month: " << month << " day: " << jday;
 
456
    s1 << ends;
 
457
    Output::Print(s1.str());
 
458
    Output::End();
 
459
 
 
460
 
 
461
    Output::Begin(); Output::Repeat(2,' ');
 
462
    ostringstream s2;
 
463
    s2.setf(ios::fixed, ios::floatfield);
 
464
    s2 << setprecision(2);
 
465
 
 
466
 
 
467
    s2 << " solar zenith angle:  " << setw(6) << asol << " deg ";
 
468
    s2 << " solar azimuthal angle:      " << setw(6) << phi0 << " deg";
 
469
    s2 << ends;
 
470
    Output::Print(s2.str());
 
471
    Output::End();
 
472
 
 
473
        
 
474
    Output::Begin(); Output::Repeat(2,' ');
 
475
    ostringstream s3;
 
476
    s3.setf(ios::fixed, ios::floatfield);
 
477
    s3 << setprecision(2);
 
478
    s3 << " view zenith angle:   " << setw(6) << avis << " deg ";
 
479
    s3 << " view azimuthal angle:       " << setw(6) << phiv << " deg ";
 
480
    s3 << ends;
 
481
    Output::Print(s3.str());
 
482
    Output::End();
 
483
    Output::Begin(); Output::Repeat(2,' ');
 
484
    ostringstream s4;
 
485
    s4.setf(ios::fixed, ios::floatfield);
 
486
    s4 << setprecision(2);
 
487
    s4 << " scattering angle:    " << setw(6) << adif << " deg ";
 
488
    s4 << " azimuthal angle difference: " << setw(6) << phi << " deg ";
 
489
    s4 << ends;
 
490
    Output::Print(s4.str());
 
491
    Output::End();
 
492
        
 
493
    Output::Begin(); Output::End();
 
494
}
 
495
 
 
496
GeomCond GeomCond::Parse()
 
497
{
 
498
    GeomCond geom;
 
499
    geom.parse();
 
500
    return geom;
 
501
}