HCastFile.cxx 66.2 KB
Newer Older
Scott Mindock's avatar
Scott Mindock committed
1
2
3
4
5
6
7
8
9
10
11
#include <cassert>
#include <iostream> // cerr, clog, endl
#include <cmath> // pow, log
#include <fcntl.h>  // open
#include <unistd.h> // close, dup, dup2, execlp, fork
#include <sys/stat.h> // fstat
#include <sys/mman.h> // mmap, munmap
#include <errno.h> // EINTR
#include <algorithm> // pair
#include <vector> // vector
#include <cstdio> // tmpfile
RH6B's avatar
RH6B committed
12
#include <sys/wait.h> // initalized=0waitpid
Scott Mindock's avatar
Scott Mindock committed
13
#include <sstream> // ostringstream
R.K.Garcia's avatar
R.K.Garcia committed
14
#include <strings.h> // memcmp
Scott Mindock's avatar
Scott Mindock committed
15
#include <cstring> // bzero
16
17
#include <stdint.h>
#include <stdlib.h> // size_t, offset_t
Scott Mindock's avatar
Scott Mindock committed
18
#include <string.h>
19

Scott Mindock's avatar
Scott Mindock committed
20
21
22
23
24
// permit external log stream
#ifndef LOG
#define LOG std::clog
#endif

Scott Mindock's avatar
Scott Mindock committed
25
// by default, cache files we've unpacked
R.K.Garcia's avatar
R.K.Garcia committed
26
#ifndef NO_LRU_CACHE
Scott Mindock's avatar
Scott Mindock committed
27
28
29
30
31
#  include "lru/lru_cache.h"
#endif

#include "aHimawariFile.hxx"

32
const double NAND = nan("");
33
// const float NANF = nanf("");
34

Scott Mindock's avatar
Scott Mindock committed
35
36
37
38
39
40
41
using namespace std;

// allow external definition of a log stream
#ifdef LOG_DECL
LOG_DECL
#endif

R.K.Garcia's avatar
R.K.Garcia committed
42
/*
Scott Mindock's avatar
Scott Mindock committed
43
44
45
46
47
48
49
#if __BYTE_ORDER__==__ORDER_LITTLE_ENDIAN__
const int LOCAL_ENDIAN = 0;
#elif __BYTE_ORDER__==__ORDER_BIG_ENDIAN__
const int LOCAL_ENDIAN=1;
#else
#  error Could Not Determine Endianness
#endif
R.K.Garcia's avatar
R.K.Garcia committed
50
*/
Scott Mindock's avatar
Scott Mindock committed
51
52
53
54
55
56
57

const unsigned MAX_IMAGE_LINES = 30000;
const unsigned MAX_IMAGE_COLUMNS = 30000;

// generated code defines the data structures
#include "hcast_file_struct.gxx"

Scott Mindock's avatar
Scott Mindock committed
58
struct count_conversion_point {
R.K.Garcia's avatar
R.K.Garcia committed
59
60
    int count;
    double bt_or_albeto;
Scott Mindock's avatar
Scott Mindock committed
61
62
63
};

struct img_compenstation_data {
R.K.Garcia's avatar
R.K.Garcia committed
64
65
66
    int line;
    double coff;
    double loff;
Scott Mindock's avatar
Scott Mindock committed
67
68
69
};

#define MAX_LINES 25000
70
#define MAX_CONVERSION_POINTS 17000  // this previously caused a bus error/segfault on HCastFile.cxx:117 for IMG_DK01B07_201507241500 when it was 1025; 14-bit B07 requires 16384
71

RH6B's avatar
RH6B committed
72
73
typedef struct conversions {
    int used;
Scott Mindock's avatar
Scott Mindock committed
74

R.K.Garcia's avatar
R.K.Garcia committed
75
76
77
78
79
80
81
82
    char scene_id[25];
    int last_line;
    double conversion_table[MAX_CONVERSION_POINTS];
    double rad_conversion_table[MAX_CONVERSION_POINTS];
    count_conversion_point count_conversion_points[MAX_CONVERSION_POINTS];
    img_compenstation_data image_compenstion_coeffs[MAX_LINES];
    ///double mjd_line_times[MAX_LINES];
    HimawariTime mjd_line_times_struct[MAX_LINES];
RH6B's avatar
RH6B committed
83

RH6B's avatar
RH6B committed
84
} conversions_type;
Scott Mindock's avatar
Scott Mindock committed
85

RH6B's avatar
RH6B committed
86
87
88
89
#define MAX_SCENE 25
static conversions_type scene_conversions[MAX_SCENE];

static int sinitalized = 0;
RH6B's avatar
RH6B committed
90
91
92
93
94
95
96
97
98
99
100
101
const HimawariTime NULL_TIME = {
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    0,
    NAN,
    NAN
};
RH6B's avatar
RH6B committed
102
void init_scene_conversions(conversions_type *p_scene_conversion) {
R.K.Garcia's avatar
R.K.Garcia committed
103
    p_scene_conversion->last_line = 0;
RH6B's avatar
RH6B committed
104
    p_scene_conversion->used=0;
R.K.Garcia's avatar
R.K.Garcia committed
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
    memset(p_scene_conversion->scene_id, 0, 25);
    //for (int i = 0; i < MAX_LINES; i++) {
    //  p_scene_conversion->mjd_line_times[i] = 0.0;
    //}
    for (int i = 0; i < MAX_LINES; i++) {
        p_scene_conversion->mjd_line_times_struct[i] =  NULL_TIME;
    }
    for (int i = 0; i < MAX_LINES; i++) {
        p_scene_conversion->image_compenstion_coeffs[i].line =
                -1;
        p_scene_conversion->image_compenstion_coeffs[i].coff = 0.0;
        p_scene_conversion->image_compenstion_coeffs[i].loff = 0.0;
    }

    //LOG<< "conversion_initialized " << std::endl;
    for (int i = 0; i < MAX_CONVERSION_POINTS; i++) {
        p_scene_conversion->count_conversion_points[i].count = -1;
        p_scene_conversion->count_conversion_points[i].bt_or_albeto = 0.0;


    }
    for (int i = 0; i < MAX_CONVERSION_POINTS; i++) {
        p_scene_conversion->conversion_table[i] = 0.0;
        p_scene_conversion->rad_conversion_table[i] = 0.0;
    }
RH6B's avatar
RH6B committed
130
131
132
133
134
135
136
137

}

// return the scene_conversion or create one
// first cut use a pool of conversion structures
conversions_type * find_scene(char * name) {
    char target[25];

R.K.Garcia's avatar
R.K.Garcia committed
138
139
140
141
142
143
144
145
146
147
148
    int scene_no = 0;
    int scene_found=-1;
    if (sinitalized == 0) {
        for (scene_no = 0; scene_no < MAX_SCENE; scene_no++) {
            //LOG<< "init_scene_conversions " << scene_no << std::endl;
            init_scene_conversions(&scene_conversions[scene_no]);
        }
        sinitalized = 1;
    }
    memset(target,0,25);
    strncpy(target,name,24);
RH6B's avatar
RH6B committed
149

RH6B's avatar
RH6B committed
150
    //LOG << "Find scene target " <<  target << std::endl;
RH6B's avatar
RH6B committed
151

R.K.Garcia's avatar
R.K.Garcia committed
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
    for (scene_no = 0; scene_no < MAX_SCENE; scene_no++) {

        if (scene_conversions[scene_no].used == 0 ) {
            //LOG<< "New scene  " << scene_no << " " << target <<std::endl;
            strncpy(scene_conversions[scene_no].scene_id, target, 24);
            scene_conversions[scene_no].used = 1;
        }

        if (strncmp(target, scene_conversions[scene_no].scene_id,24) == 0) {
            scene_found = scene_no;

            //LOG<< "Scene found " << scene_found <<  " " << target <<std::endl;
            break;
        }
    }
    if (scene_found < 0 ) {
        LOG<< "Warning No scene found " <<   target <<std::endl;
        return NULL;
    }
    return &(scene_conversions[scene_found]);
RH6B's avatar
RH6B committed
172
}
Scott Mindock's avatar
Scott Mindock committed
173

RH6B's avatar
RH6B committed
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
/***

pro planck_ahi, T, chan

chan_idx = chan - 7


;AHI Planck constants in (cm-1) space
 fk1 = [203539.00, 49635.30, 35718.60, 30051.90, 18805.40, 13324.60, 10581.50, 8417.47, 6311.93, 5092.73]   ;16

 fk2 =[3705.95, 2315.34, 2074.83, 1958.73, 1675.37, 1493.60, 1383.14, 1281.58, 1164.32, 1083.93]

 bc1 =[0.45925, 1.62641, 0.30465, 0.05552, 0.13189, 0.09194, 0.08686, 0.17399, 0.23460, 0.05990]

 bc2 =[0.99936, 0.99649, 0.99927, 0.99986, 0.99963, 0.99971, 0.99971, 0.99938, 0.99908, 0.99975]


;Wavenumber
 nu =[2575.74, 1609.03, 1442.11, 1361.42, 1164.48, 1038.14, 961.39, 890.83, 809.40, 753.42]



;Calculate Radiance in NOAA (cm-1) units
adjusted_temp = bc1(chan_idx) + bc2(chan_idx) * T

Rad = fk1(chan_idx) / (exp(fk2(chan_idx) / adjusted_temp) - 1.)

print, "Radiance, NOAA UNITS (cm-1) = ", Rad


;convert into JMA/NASA (um) units
rad_noaa= Rad / (((10000.0 / nu(chan_idx) )^2) / 10.0)

print, "Radiance, JMA/NASA UNITS (um) = ", rad_noaa

stop
end

***/

//;AHI Planck constants in (cm-1) space
215
216
217
218
// static const double fk1[10] = {203539.00, 49635.30, 35718.60, 30051.90, 18805.40, 13324.60, 10581.50, 8417.47, 6311.93, 5092.73};
// static const double fk2[10] = {3705.95, 2315.34, 2074.83, 1958.73, 1675.37, 1493.60, 1383.14, 1281.58, 1164.32, 1083.93};
// static const double bc1[10] = {0.45925, 1.62641, 0.30465, 0.05552, 0.13189, 0.09194, 0.08686, 0.17399, 0.23460, 0.05990};
// static const double bc2[10] = {0.99936, 0.99649, 0.99927, 0.99986, 0.99963, 0.99971, 0.99971, 0.99938, 0.99908, 0.99975};
RH6B's avatar
RH6B committed
219
220

//;Wavenumber
221
// static const double nu[10] ={2575.74, 1609.03, 1442.11, 1361.42, 1164.48, 1038.14, 961.39, 890.83, 809.40, 753.42};
222

223
// convert raw counts to radiance
224
// extracted from table
225
// SHA256(Himawari-8_HRIT_table.xlsx)= 1c95ee65dae21f3d8bff956e7643704a01a29a25951cc49b7c6664e9687794d7
226
227
228
229
230
231
// copy column to column.txt and run this python snippet:
//   from numpy import array, arange, polyfit
//   ys = open("column.txt", 'r').read()
//   y = array([float(n) for n in ys.split() if n])
//   polyfit(arange(len(y),dtype=float), y, 2)
// x**2 (confirmed near 0); then x**1, x**0 terms are then copied into this table for Himawari-8:
232
static const double rad_mb[16][2] = {{ NAND, 0.0 },  // B01
233
                                     { NAND, 0.0 },  // B02
234
235
236
237
                                     { 0.50817698, -0.51934652 },  // B03
                                     { 0.30270543, -0.30935355 },  // B04
                                     { 0.07548224, -0.07713647 },  // B05
                                     { 0.02340196, -0.02391235 },  // B06
238
239
240
241
242
243
244
                                     { -1.24342320e-03, 1.27202222 },  // B07 IR4
                                     { -5.64755417e-03, 5.77772475 },  // B08 IR3
                                     { -1.35314435e-02, 1.38435591e+01 },  // B09
                                     { -1.44305111e-02, 1.47640133e+01 },  // B10
                                     { -1.56364563e-02, 1.60025993e+01 },  // B11
                                     { -1.53364798e-02, 1.57039922e+01 },  // B12
                                     { -1.46494251e-02, 1.50098803e+01 },  // B13 IR1
245
                                     { -1.36611465e-02, 1.40105473e+01 },  // B14
246
                                     { -1.21365498e-02, 1.24694285e+01 },  // B15 IR2
247
                                     { -1.09186756e-02, 1.12390993e+01 }   // B16
248
249
                                    };  // NOTE: H-8 spreadsheet

250
// high bit-rate B07 with 14 bits of resolution
251
static const double rad_b07_hbr_mb[] = { -7.87614900e-04,   1.28775036e+01 };
252

253
// static const double alb_mb[2] = {0.00097849, -0.001};
RH6B's avatar
RH6B committed
254

255
256
257
258
259
// static double bt_to_rad( double T, int band ) {
//     int chan_idx = band-7;
//     double adjusted_temp;
//     double Rad;
//     double rad_noaa;
RH6B's avatar
RH6B committed
260

261
262
//     //;Calculate Radiance in NOAA (cm-1) units
//     adjusted_temp = bc1[chan_idx] + bc2[chan_idx] * T;
RH6B's avatar
RH6B committed
263

264
//     Rad = fk1[chan_idx] / (exp(fk2[chan_idx] / adjusted_temp) - 1.);
RH6B's avatar
RH6B committed
265

266
//     ///LOG << "Radiance, NOAA UNITS (cm-1) = " << Rad << std::endl;
RH6B's avatar
RH6B committed
267
268


269
270
//     //;convert into JMA/NASA (um) units
//     rad_noaa = Rad / (  pow((10000.0 / nu[chan_idx] ),2.0) / 10.0);
RH6B's avatar
RH6B committed
271

272
//     ///LOG << "Radiance,JMA/NASA UNITS (um) = " << rad_noaa << std::endl;
RH6B's avatar
RH6B committed
273

274
275
//     return rad_noaa;
// }
276

277
// static double alb_to_rad( double A, int band ) {
R.K.Garcia's avatar
R.K.Garcia committed
278
//  // reverse transform
RH6B's avatar
RH6B committed
279

280
// }
281

282
static double counts_to_rad(double x, int band, bool high_bit_rate=false) {
R.K.Garcia's avatar
R.K.Garcia committed
283
284
285
286
287
288
289
290
291
292
293
294
    assert(band <= 16);
    assert(band >= 1);

    if (high_bit_rate==true && band==7) {
        const double m(rad_b07_hbr_mb[0]),
                               b(rad_b07_hbr_mb[1]);
        return m * x + b;
    }

    const double m(rad_mb[band-1][0]),
                           b(rad_mb[band-1][1]);
    return m * x + b;
295
296
297
}


298
299
// this used data from message #3  to create a table for conversion
// of counts to bt or allbeto
RH6B's avatar
tweek    
RH6B committed
300
static int fill_conversion_table(double * in_conv_tbl,
R.K.Garcia's avatar
R.K.Garcia committed
301
302
303
                                 count_conversion_point* in_conv_pnts ) {
    int table_size = 0;
    int table_idx;
Scott Mindock's avatar
Scott Mindock committed
304

R.K.Garcia's avatar
R.K.Garcia committed
305
//  LOG<< "fill_conversion_table" << std::endl;
Scott Mindock's avatar
Scott Mindock committed
306

R.K.Garcia's avatar
R.K.Garcia committed
307
308
309
310
    for (table_idx = 0; table_idx < MAX_CONVERSION_POINTS; table_idx++) {
        if (in_conv_pnts[table_idx].count != -1) {
            table_size++;
///         LOG<< "table_idx " << table_idx << " "<< in_conv_pnts[table_idx].bt_or_albeto << std::endl;
Scott Mindock's avatar
Scott Mindock committed
311

R.K.Garcia's avatar
R.K.Garcia committed
312
313
314
        }
    }
    //LOG<< "fill_conversion_table with " << table_size << " Points"<< MAX_CONVERSION_POINTS<< std::endl;
Scott Mindock's avatar
Scott Mindock committed
315

R.K.Garcia's avatar
R.K.Garcia committed
316
317
318
319
    for (table_idx = 0; table_idx < MAX_CONVERSION_POINTS; table_idx++) {
        in_conv_tbl[table_idx] = in_conv_pnts[table_size
                - 1].bt_or_albeto;
    }
Scott Mindock's avatar
Scott Mindock committed
320

R.K.Garcia's avatar
R.K.Garcia committed
321
322
323
324
    double maxa = 0.0;
    double mina = 500.0;
    for (int range = 0; range < table_size - 2; range++) {
        // stepwise interpolate through the table
cent6's avatar
cent6 committed
325

R.K.Garcia's avatar
R.K.Garcia committed
326
327
        int r_begin = in_conv_pnts[range].count;
        int r_end = in_conv_pnts[range + 1].count;
Scott Mindock's avatar
Scott Mindock committed
328

R.K.Garcia's avatar
R.K.Garcia committed
329
330
        double v_begin = in_conv_pnts[range].bt_or_albeto;
        double v_end = in_conv_pnts[range + 1].bt_or_albeto;
Scott Mindock's avatar
Scott Mindock committed
331

cent6's avatar
cent6 committed
332
        // the number of indexes
R.K.Garcia's avatar
R.K.Garcia committed
333
        int points = r_end - r_begin;
Scott Mindock's avatar
Scott Mindock committed
334

cent6's avatar
cent6 committed
335
        // The change in value
R.K.Garcia's avatar
R.K.Garcia committed
336
337
        double span = v_end - v_begin;
        double step = (double) span / (double) points;
Scott Mindock's avatar
Scott Mindock committed
338

R.K.Garcia's avatar
R.K.Garcia committed
339
        //LOG<<in_conv_pnts[range].bt_or_albeto<<" "<< range<< " r_begin "<< r_begin << " " << r_end << " " << points << " " << step<< " " << span<< std::endl;
Scott Mindock's avatar
Scott Mindock committed
340

341
        // int use=r_begin-1;
R.K.Garcia's avatar
R.K.Garcia committed
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
        for (int count = 0; count < points; count++) {
            in_conv_tbl[r_begin + count] =
                    in_conv_pnts[range].bt_or_albeto
                            + step * ((double) count);


            if (in_conv_tbl[r_begin + count] > maxa) {
                maxa = in_conv_tbl[r_begin + count];
                //LOG<<"maxa "<< maxa<<" "<< r_begin+count<< std::endl;
            }
            if (in_conv_tbl[r_begin + count] < mina) {
                mina = in_conv_tbl[r_begin + count];
                //LOG<<"mina "<< mina<<" "<< r_begin+count<< std::endl;
            }
        }

        //LOG<< "radoralb "<<range<< "=" << in_conv_pnts[range].bt_or_albeto << " val " << in_conv_tbl[r_begin]<< " min "<< mina << " max "<< maxa<< std::endl;
        //LOG<< "radoralb "<<range+1<< "=" << in_conv_pnts[range+1].bt_or_albeto << std::endl;
    }

//  LOG<< "fill_conversion_table " << table_size << std::endl;
    return 0;
Scott Mindock's avatar
Scott Mindock committed
364
365
}

RH6B's avatar
RH6B committed
366
struct HCastFileInternal {
R.K.Garcia's avatar
R.K.Garcia committed
367
368
369
370
    int fd; // numerical file descriptor
    void *mm; // memory map start
    off_t size; // how big it is
    struct HimawariSentinels sentinel;
371
    int bpp;  // discovered bits-per-pixel based on translation table content
R.K.Garcia's avatar
R.K.Garcia committed
372
    class HCastNavigationCorrectionTransform *correction; // if NULL, return direct storage values without correction
RH6B's avatar
RH6B committed
373
374
375
376
};



Scott Mindock's avatar
Scott Mindock committed
377
struct HCastFile {
R.K.Garcia's avatar
R.K.Garcia committed
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
    struct HCastFileInternal *impl; // bookkeeping
    primary_header *basic;
    image_structure *data;
    image_navigation *proj;
    image_data_function *idata;
    image_annotation *iano;
    timestamp *itimestamp;
    image_segment *segment;
    image_compensation *compensation;
    image_observation * observation;
    timestamp *obs_time;
    image_quality *error;
    uint16_t *images;

    /* information parsed from text data in file
     * FUTURE: move these into impl, or move contents of impl into this private structure
     *
     */
    std::string hcastname;
    int band_number;
cent6's avatar
cent6 committed
398
399

    conversions_type *scene_conversion_meta;
Scott Mindock's avatar
Scott Mindock committed
400

Scott Mindock's avatar
Scott Mindock committed
401
402
403
404
405
};

#include "hcast_file_iostream.gxx"

inline std::ostream &operator<<(std::ostream &out, const HCastFile &hf) {
R.K.Garcia's avatar
R.K.Garcia committed
406
407
408
409
410
411
412
413
414
    out << "<HCastFile>\n";
    out << *hf.basic;
    out << *hf.data;
    out << *hf.proj;
    out << *hf.segment;
    out << *hf.obs_time;
    out << *hf.error;
    out << "</HCastFile>\n";
    return out;
Scott Mindock's avatar
Scott Mindock committed
415
416
}

R.K.Garcia's avatar
R.K.Garcia committed
417
/*
Scott Mindock's avatar
Scott Mindock committed
418
419
// FUTURE: boot these log functions elsewhere
static void error(const string &txt) {
R.K.Garcia's avatar
R.K.Garcia committed
420
    clog << "ERROR: " << txt << endl;
Scott Mindock's avatar
Scott Mindock committed
421
422
423
}

static void debug(const string &txt) {
R.K.Garcia's avatar
R.K.Garcia committed
424
    clog << "DEBUG: " << txt << endl;
Scott Mindock's avatar
Scott Mindock committed
425
}
R.K.Garcia's avatar
R.K.Garcia committed
426
*/
Scott Mindock's avatar
Scott Mindock committed
427
static void error(const char *txt) {
R.K.Garcia's avatar
R.K.Garcia committed
428
    clog << "ERROR: " << txt << endl;
Scott Mindock's avatar
Scott Mindock committed
429
430
431
432
}

//
// endian_swap for all the himawari blocks
Scott Mindock's avatar
Scott Mindock committed
433
//
Scott Mindock's avatar
Scott Mindock committed
434
435

inline void endian_swap(uint8_t *p, size_t elem, size_t count) {
R.K.Garcia's avatar
R.K.Garcia committed
436
437
438
439
440
441
442
443
444
445
446
447
448
449
    if (elem <= 1)
        return;
    assert((elem & 1) == 0);
    assert(elem <= 16);
    uint8_t t;
    uint8_t *x = p;
    for (; count > 0; --count, x += elem) {
        int d = elem - 1;
        for (unsigned b = 0; b < elem / 2; ++b, --d) {
            t = x[d];
            x[d] = x[b];
            x[b] = t;
        }
    }
Scott Mindock's avatar
Scott Mindock committed
450
451
452
453
}

template<typename BasicType>
void endian_swap(BasicType &q) {
R.K.Garcia's avatar
R.K.Garcia committed
454
    endian_swap((uint8_t *) &q, sizeof(BasicType), 1);
Scott Mindock's avatar
Scott Mindock committed
455
456
}

RH6B's avatar
RH6B committed
457
struct HCastFilscene_conversioneInternal {
R.K.Garcia's avatar
R.K.Garcia committed
458
459
460
461
462
    int fd; // numerical file descriptor
    void *mm; // memory map start
    off_t size; // how big it is
    struct HimawariSentinels sentinel;
    class HCastNavigationCorrectionTransform *correction; // if NULL, return direct storage values without correction
Scott Mindock's avatar
Scott Mindock committed
463
464
465
466
467
};

// bring in the autogenerated endian_swap routines
#include "hcast_file_endian.gxx"

R.K.Garcia's avatar
R.K.Garcia committed
468
// we count zeroes as use-defaults
Scott Mindock's avatar
Scott Mindock committed
469
470
471
// counts values that are still zero are then replaced with corresponding file fields
// FUTURE: make this easier to follow since it's happening in three places (here, attach_basic, map_contents_with_native_order)
static HimawariSentinels merge_sentinels(const HimawariSentinels &a,
R.K.Garcia's avatar
R.K.Garcia committed
472
473
        const HimawariSentinels &b) {
    //derived_invalid=65535;
Scott Mindock's avatar
Scott Mindock committed
474

R.K.Garcia's avatar
R.K.Garcia committed
475
    HimawariSentinels zult;
Scott Mindock's avatar
Scott Mindock committed
476
#   define MERGE(field) zult.field = a.field!=0? a.field : b.field
R.K.Garcia's avatar
R.K.Garcia committed
477
478
479
480
    MERGE(count_invalid);
    MERGE(count_outside_scan);
    MERGE(derived_invalid);  // this causes problem
    MERGE(derived_outside_scan);
Scott Mindock's avatar
Scott Mindock committed
481
#   undef MERGE
R.K.Garcia's avatar
R.K.Garcia committed
482
    return zult;
Scott Mindock's avatar
Scott Mindock committed
483
484
485
486
}

#ifndef NO_LRU_CACHE

Scott Mindock's avatar
Scott Mindock committed
487
const char *ENV_LRU_DEPTH2 = "HIMAWARI_DECOMPRESS_CACHE_SIZE";
Scott Mindock's avatar
Scott Mindock committed
488
489
490
const int DEFAULT_LRU_DEPTH = 40;

// least-recently-used file descriptor cache value:
R.K.Garcia's avatar
R.K.Garcia committed
491
// retain a duplicate file descriptor for an anonymous (tmp) file.
Scott Mindock's avatar
Scott Mindock committed
492
// permit duplicates to be created, replicating open() behavior.
R.K.Garcia's avatar
R.K.Garcia committed
493
494
// when file descriptor falls off LRU, open_fd id destructed
// and private file descriptor is closed.
Scott Mindock's avatar
Scott Mindock committed
495
class open_fd {
R.K.Garcia's avatar
R.K.Garcia committed
496
    const int the_fd;
Scott Mindock's avatar
Scott Mindock committed
497
public:
R.K.Garcia's avatar
R.K.Garcia committed
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
    explicit open_fd(int fd) :
            the_fd(::dup(fd)) {
        assert(fd >= 0);
        assert(the_fd >= 0);
    }

    // default constructor
    open_fd() :
            the_fd(-1) {
    }

    // copy constructor
    open_fd(const open_fd &it) :
            the_fd(it.the_fd >= 0 ? ::dup(it.the_fd) : -1) {
        assert(it.the_fd < 0 || the_fd >= 0);
    }

    int dup() const {
        return ::dup(the_fd);
    }

    ~open_fd() {
    /**
    }
        if (the_fd >= 0)

                if (scene_conversions[scene_no].used == 0 ) {
            //LOG<< "New scene  " << scene_no << " " << target <<std::endl;
            strncpy(scene_conversions[scene_no].scene_id, target, 24);
            scene_conversions[scene_no].used = 1;
        }
cent6's avatar
cent6 committed
529
**/
530
            // LOG<< "close  "  <<std::endl;
cent6's avatar
cent6 committed
531
532


R.K.Garcia's avatar
R.K.Garcia committed
533
534
            ::close(the_fd);
    }
Scott Mindock's avatar
Scott Mindock committed
535
536
537
538
539
};

#endif

static int envint(const char *key, int default_value) {
R.K.Garcia's avatar
R.K.Garcia committed
540
    const char *env = getenv(key);
RH6B's avatar
RH6B committed
541

R.K.Garcia's avatar
R.K.Garcia committed
542
543
544
545
546
    if (env == 0)
        return default_value;
    int value = default_value;
    int rc = sscanf(env, "%d", &value);
    return (rc == 1) ? value : default_value;
Scott Mindock's avatar
Scott Mindock committed
547
548
549
}

// open files that are bzipped or gzipped
Scott Mindock's avatar
Scott Mindock committed
550
int opensesame2(const char *path, int oflag, int mode) {
R.K.Garcia's avatar
R.K.Garcia committed
551
552
    string p(path);
    const char *util = 0;
Scott Mindock's avatar
Scott Mindock committed
553

R.K.Garcia's avatar
R.K.Garcia committed
554
555
556
557
558
559
560
    if (p.rfind(".bz2") == p.length() - 4) {
        util = "bunzip2";
    } else if (p.rfind(".gz") == p.length() - 3) {
        util = "gunzip";
    } else {
        return open(path, oflag, mode);
    }
Scott Mindock's avatar
Scott Mindock committed
561
562

#ifndef NO_LRU_CACHE
R.K.Garcia's avatar
R.K.Garcia committed
563
564
565
566
567
    // include internal least-recently-used cache for files that we had to decompress
    static LRUCache<string, open_fd> cache(
            envint(ENV_LRU_DEPTH2, DEFAULT_LRU_DEPTH));

    if (cache.exists(p)) {
568
        // clog << "decompression cache hit: " << p << endl;
R.K.Garcia's avatar
R.K.Garcia committed
569
570
571
        cache.touch(p);
        return cache.fetch(p).dup();
    }
Scott Mindock's avatar
Scott Mindock committed
572
573
#endif

R.K.Garcia's avatar
R.K.Garcia committed
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
    FILE *fp = 0;
    int fdcomp = -1;
    fp = tmpfile();

    assert(fp != 0);
    fdcomp = dup(fileno(fp));
    assert(fdcomp >= 0);
    fclose(fp);

    /* spin off a process to unpack it into anonymous file */
    int pid = fork();
    if (pid == 0) { // then I am the child
        int fd = open(path, O_RDONLY, 0);
        if (fd < 0) {
            exit(1);
        }
        int rc = dup2(fd, 0); // stdin
        assert(rc >= 0);
        rc = dup2(fdcomp, 1); //stdout
        assert(rc >= 0);
        close(fd);
        close(fdcomp);
        execlp(util, util, 0, NULL);
        exit(1); // never occurs unless execlp fails
    } else if (pid > 0) {
        int rc = 0;
        waitpid(pid, &rc, 0);
        if (rc == 0) {
            goto success;
        } else {
            goto fail;
        }
    }
    // continue if fork failed, which should never happen unless system is leaking process ids
    cerr
            << "ERROR: decompression process fork failed (very unusual), unable to open "
            << p << endl;

    fail: close(fdcomp);
    return -1;

    success:
Scott Mindock's avatar
Scott Mindock committed
616
#ifndef NO_LRU_CACHE
R.K.Garcia's avatar
R.K.Garcia committed
617
    cache.insert(p, open_fd(fdcomp)); // retains its own duplicated file descriptor
Scott Mindock's avatar
Scott Mindock committed
618
#endif
R.K.Garcia's avatar
R.K.Garcia committed
619
    return fdcomp;
Scott Mindock's avatar
Scott Mindock committed
620
621
622
623
624
625
626

}

//
// Attaching the file and populating the public data structure
//

627
static HCastFile *attach_basic(const char *pathname) {
Scott Mindock's avatar
Scott Mindock committed
628

R.K.Garcia's avatar
R.K.Garcia committed
629
    //LOG<< "Attach " << pathname << std::endl;
Scott Mindock's avatar
Scott Mindock committed
630

RH6B's avatar
RH6B committed
631
    conversions_type *scene_conversion = NULL;
R.K.Garcia's avatar
R.K.Garcia committed
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
    char *pch;
    char * name=NULL;
    const char *delimit= "/";
    char buffer[500];
    memset(buffer,0,499);
    strcpy( buffer,pathname);

    pch = strtok ((char*)buffer,delimit);

    while (pch != NULL)
    {
        name=pch;
        pch = strtok (NULL, delimit);
    }
    scene_conversion = find_scene(name);

    if (scene_conversion == NULL ) {
        LOG<< "No scene conversion allocated or found " << pathname << std::endl;
        return 0;
    }

    // open the file
    int fd = opensesame2(pathname, O_RDONLY, 0);
    if (fd < 0) {
        LOG<< "Could not attach " << pathname << " as a HCast file" << std::endl;
        return 0;
    }

    // get information about it
    struct stat stats;
    int rc = fstat(fd, &stats);
    if (rc != 0) {
        LOG<< "invalid file status from " << pathname << std::endl;
        close(fd);
        return 0;
    }

    /*    if (rc!=0 || stats.st_size<off_t(sizeof(cal_info))) {
     LOG << "invalid file status from " << pathname << std::endl;
     close(fd);
     return 0;
     }*/

    // memory-map the entire contents
    void *mm = mmap(0 /*anywhere*/, (size_t) stats.st_size,
            PROT_READ | PROT_WRITE, MAP_PRIVATE, fd, (off_t) 0);
    if (mm == MAP_FAILED) {
        LOG<< "unable to memory-map " << pathname << std::endl;
        close(fd);
        return 0;
    }

    HCastFileInternal *impl = new HCastFileInternal;
    HCastFile *it = new HCastFile;

    impl->fd = fd;
    impl->mm = mm;
    impl->size = (off_t) stats.st_size;
    impl->correction = 0;
    it->impl = impl;
    it->scene_conversion_meta = scene_conversion;

    return it;
Scott Mindock's avatar
Scott Mindock committed
695
696
697
698
}

template<typename Struct>
Struct *mmap_struct(void *mmap, off_t &offset, bool flipping,
R.K.Garcia's avatar
R.K.Garcia committed
699
700
        uint8_t block_number) {
    Struct *it = (Struct *) ((int8_t *) mmap + offset);
Scott Mindock's avatar
Scott Mindock committed
701

R.K.Garcia's avatar
R.K.Garcia committed
702
    if (block_number != it->type) {
RH6B's avatar
RH6B committed
703

R.K.Garcia's avatar
R.K.Garcia committed
704
705
706
707
708
709
710
711
712
        if ( block_number != 0)
        {
          //  LOG<< "hcastfile unexpected block number: wanted " << int(block_number) << " and received " << int(it->type) << std::endl;
        }
        return 0;
    }
    if (flipping) {
        endian_swap(*it); // NOTE: this will also flip the trailing arrays automagically for complex structs
    }
Scott Mindock's avatar
Scott Mindock committed
713

R.K.Garcia's avatar
R.K.Garcia committed
714
    // LOG<< "hcastfile  block number: wanted " << int(block_number) << " and received " << int(it->type) << " " << int(it->length) << std::endl;
Scott Mindock's avatar
Scott Mindock committed
715

R.K.Garcia's avatar
R.K.Garcia committed
716
717
    offset += it->length;
    return it;
Scott Mindock's avatar
Scott Mindock committed
718
719
720
}

class HCastNavigationCorrectionTransform {
R.K.Garcia's avatar
R.K.Garcia committed
721
    // implement position correction equivalent to
R.K.Garcia's avatar
R.K.Garcia committed
722
723
    const double colcent, lincent; // columns, rows
            // rotang;  // radians
R.K.Garcia's avatar
R.K.Garcia committed
724
725
726
727
728
729
730
731
732
733
734
735
736
737
    const double cosrotang, sinrotang;

    const int first_line_of_image_segment;
    const int number_of_lines, number_of_columns;

    struct {
        int first_line_to_correct;
        vector<pair<double, double> > lin_col_adjustment;
    } corr;

    void initialize_adjustment_table(const struct HCastFile *it) {

        return;
        // equivalent of hisd_comp_table
R.K.Garcia's avatar
R.K.Garcia committed
738
        // corr.first_line_to_correct = -1;
R.K.Garcia's avatar
R.K.Garcia committed
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
        /*        if (it->nav_corr->number_of_correction_information_data==0)
         return;
         corr.lin_col_adjustment.reserve(550);
         */
        //      nav_corr_entry cto = {0,0.0,0.0};
        // iterate through the correction table, linearly interpolating the line and column corrections
        /*        for(int i=0; i<int(it->nav_corr->number_of_correction_information_data)-1; ++i) {
         const nav_corr_entry &cfrom(it->nav_corr->corrections[i]);
         cto = it->nav_corr->corrections[i+1];
         if (-1==corr.first_line_to_correct)
         corr.first_line_to_correct = cfrom.line_number_after_rotation;

         for(int line=cfrom.line_number_after_rotation; line<cto.line_number_after_rotation; ++line) {
         double factor = (double(line) - double(cfrom.line_number_after_rotation)) /
         (double(cto.line_number_after_rotation) - double(cfrom.line_number_after_rotation));

         double lincor = (double(cto.shift_amount_for_line_direction) - double(cfrom.shift_amount_for_line_direction))
         * factor + double(cfrom.shift_amount_for_line_direction);

         double colcor = (double(cto.shift_amount_for_column_direction) - double(cfrom.shift_amount_for_column_direction))
         * factor + double(cfrom.shift_amount_for_column_direction);

         corr.lin_col_adjustment.push_back(pair<double, double>(lincor, colcor));
         }
         }
         */
        // handle last line
        /*
         pair<double, double> last(cto.shift_amount_for_line_direction, cto.shift_amount_for_column_direction);
         corr.lin_col_adjustment.push_back(last);
         */
    }

    pair<double, double> lin_col_adjustment_lookup(double line) const {
        const pair<double, double> NO_CORRECTION(0., 0.);
        if (corr.first_line_to_correct == -1) // then no table present!
            return NO_CORRECTION;
        int offset = int(line + 0.5) - corr.first_line_to_correct;
        if (offset < 0 || offset >= int(corr.lin_col_adjustment.size()))
            return NO_CORRECTION;
        return corr.lin_col_adjustment[offset];
    }
Scott Mindock's avatar
Scott Mindock committed
781
782
783

public:

R.K.Garcia's avatar
R.K.Garcia committed
784
    explicit HCastNavigationCorrectionTransform(const HCastFile *it) :
R.K.Garcia's avatar
R.K.Garcia committed
785
            colcent(0), lincent(0), cosrotang(0), sinrotang(0), first_line_of_image_segment(
R.K.Garcia's avatar
R.K.Garcia committed
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
                    it->segment->line_no_image_segm), number_of_lines(
                    it->data->number_of_lines), number_of_columns(
                    it->data->number_of_columns) {
        initialize_adjustment_table(it);
    }

    /// Transform requested 1-based scene line and column in the output image,
    // into storage line and column inside the file
    // will return -1 as line, column if we're out of storage bounds
    inline pair<double, double> coord_adjust(double line_s1,
            double column_s1) const {
        // apply correction table shift
        pair<double, double> adjust = lin_col_adjustment_lookup(line_s1);
        double shifted_lin = line_s1 - adjust.first;
        double shifted_col = column_s1 - adjust.second;

        // adjusted_lin, _col are in scene coordinate system (1-based) - convert from 0-based scene
        double adjusted_lin = (shifted_lin - lincent) * cosrotang
                - (shifted_col - colcent) * sinrotang + lincent;
        double adjusted_col = (shifted_col - colcent) * cosrotang
                + (shifted_lin - lincent) * sinrotang + colcent;

        return pair<double, double>(adjusted_lin, adjusted_col);
    }

    pair<int, int> index_adjust(double line_s1, double column_s1) const {
        pair<double, double> adjusted = coord_adjust(line_s1, column_s1);

        // storage_lin, _col are 0-based storage offsets for packed data in the file
        int storage_lin = int(adjusted.first + 0.5)
                - first_line_of_image_segment;
        int storage_col = int(adjusted.second + 0.5) - 1;

        if ((storage_lin < 0 || storage_lin >= number_of_lines)
                || (storage_col < 0 || storage_col >= number_of_columns)) {
            storage_lin = storage_col = -1;
        }
        return pair<int, int>(storage_lin, storage_col);
    }
Scott Mindock's avatar
Scott Mindock committed
825
826
827

};

828
829
// set default missing value counts to maxint-1 according to bits_per_pixel
// typically this is 1023 (10 bits per pixel) or 16383 (14 bits per pixel)
830
static int init_sentinels_bpp(HCastFile *it, const struct HimawariSentinels *sent, int max_count) {
831
832
    //int32_t missing_value = (1L << (int32_t)it->data->bits_per_pixel) - 1;  // always results in 65535 >_<#
    int32_t bpp = int32_t(::ceil(::log(max_count) / ::log(2.0)));
833
834
835
    it->impl->bpp = bpp;
    // int32_t missing_value = (1L << bpp) - 1;
    int32_t missing_value = max_count;
836
837
838
839
840
841
842
    struct HimawariSentinels dflt = DEFAULT_HIMAWARI_SENTINELS;
    dflt.count_invalid = dflt.count_outside_scan = missing_value;
    if (sent==0) {
        it->impl->sentinel = dflt;
    } else {
        it->impl->sentinel = merge_sentinels(dflt, *sent);
    }
843
    // LOG << "missing value: " << missing_value << " for bpp: " << bpp << " max_count: " << max_count << "\n";
844
845
    return missing_value;
}
RH6B's avatar
RH6B committed
846

847
static int map_contents_with_native_order(HCastFile *it, const HimawariSentinels *optional_sentinels) {
R.K.Garcia's avatar
R.K.Garcia committed
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
    // map basic block
    void *mm = it->impl->mm;
    char buff2[10000];
    char *bufp = NULL;
    bool flipping = true;
    off_t offset1 = 0;

    conversions_type *p_scene_conversion = it->scene_conversion_meta;

    mmap_struct<primary_header>(mm, offset1, flipping, 0);

    mmap_struct<primary_header>(mm, offset1, flipping, 0);

    //if (0 == ret)
    //  LOG<< "exit hcastfile  " << ret<<  std::endl;
    //  return ERROR_ATTACH_DATA_BLOCK;
    ///return ERROR_ATTACH_DATA_BLOCK;
    it->basic = (primary_header *) mm;
    off_t offset = it->basic->length;

    // do we have to flip bytes?
    //bool flipping = (LOCAL_ENDIAN != it->basic->byte_order);

    // FUTURE: get a real logging system
    // LOG << "byte order match: " << flipping << std::endl;

    // data_info *data;
    it->data = mmap_struct<image_structure>(mm, offset, flipping, 1);
    if (0 == it->data)
        return ERROR_ATTACH_DATA_BLOCK;

    //LOG<< "image_structure 1 " << it->data->number_of_columns<< " " << it->data->number_of_lines<< std::endl;

    // proj_info *proj;
    it->proj = mmap_struct<image_navigation>(mm, offset, flipping, 2);
    if (0 == it->proj)
        return ERROR_ATTACH_PROJ_BLOCK;

    //LOG<< "image_structure 1 " << it->proj->COFF<< " " << it->proj->LOFF<< std::endl;
    //LOG<< "image_structure 1 " << it->proj->projection_name<< std::endl;

    // image_data_function *idata;
    it->idata = mmap_struct<image_data_function>(mm, offset, flipping, 3);
    if (0 == it->idata)
        return ERROR_ATTACH_PROJ_BLOCK;

    int numchars = it->idata->length - 3;

    char buff[10000];
    memset(buff, 0, 10000);
    strncpy(buff, (char*) &(it->idata->data_definition_block), numchars);

    for (char *t = buff; *t; t++) {
        if (*t == '\r') {
            *t = '\n';
        }
    }
    //LOG<< "Data_definition_block 3" << std::endl;
    //LOG<< "idata2 " << buff<< std::endl;

    int rc = 2;
    int table_line = 0;
    int now = 0;
911
    int max_count = 0;  // maximum declared counts, used to figure out missing value sentinel
R.K.Garcia's avatar
R.K.Garcia committed
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934

    while (now < numchars) {
        int before = now;
        while (buff[now] != '\n')
            now++;
        now++;

        //LOG<< "before " << before << " "  << now << " "<<now-before << " "  << numchars  <<std::endl;

        int count;
        double bt_or_albeto;

        char b2[400];
        memset(b2, 0, 400);
        strncpy(b2, &buff[before], (now - before) - 1);
        rc = sscanf(b2, "%d:=%lf", &count, &bt_or_albeto);

        if (rc == 2) {
            p_scene_conversion->count_conversion_points[table_line].count = count;
            p_scene_conversion->count_conversion_points[table_line].bt_or_albeto =
                    bt_or_albeto;
            //LOG<< "Data Def 3 table_line " << table_line << " " << count << " " << bt_or_albeto <<std::endl;
            table_line++;
935
            if ((count > max_count) && (count != 65535)) {  // ignore final value and get "designated value"
936
937
                max_count = count;
            }
R.K.Garcia's avatar
R.K.Garcia committed
938
939
940
        }
    }

941
942
    // use discovered maximum count to set sentinel value and bits-per-pixel
    init_sentinels_bpp(it, optional_sentinels, max_count);
943

R.K.Garcia's avatar
R.K.Garcia committed
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
    fill_conversion_table(p_scene_conversion->conversion_table,
            p_scene_conversion->count_conversion_points);

    // image_annotation *iano;
    it->iano = mmap_struct<image_annotation>(mm, offset, flipping, 4);
    if (0 == it->iano)
        return ERROR_ATTACH_PROJ_BLOCK;

    numchars = it->iano->length - 3;

    memset(buff, 0, 10000);
    strncpy(buff, (char*) &(it->iano->data_definition_block), numchars);

    for (char *t = buff; *t; t++) {
        if (*t == '\r') {
            *t = '\n';
        }
    }

    it->hcastname = std::string(buff);
    //LOG<< "image_annotation 4 " << numchars << std::endl;
    //LOG<< "idata2 " << buff<< std::endl;

    //LOG<< "iano " << buff<< it->hcastname<< std::endl;

    // timestamp *itimestamp;
    it->itimestamp = mmap_struct<timestamp>(mm, offset, flipping, 5);
    if (0 == it->itimestamp)
        return ERROR_ATTACH_PROJ_BLOCK;

    // image_segment *segment;
    it->segment = mmap_struct<image_segment>(mm, offset, flipping, 128);
    if (0 == it->segment)
        return ERROR_ATTACH_SEGMENT_BLOCK;

    // image_compensation *compensation;
    it->compensation = mmap_struct<image_compensation>(mm, offset, flipping,
            130);
    if (0 == it->compensation)
        return ERROR_ATTACH_SEGMENT_BLOCK;

    //LOG<< "image_compensation 130" << std::endl;

    bufp = buff;
    memset(buff, 0, 10000);
    numchars = it->compensation->length - 3;
    strncpy(buff, (char*) &(it->compensation->image_compensation_information),
            numchars);

    for (char *t = buff; *t; t++) {
        if (*t == '\r') {
            *t = '\n';
        }
    }

//  LOG<< buff<< std::endl;