IASI Code Collection
libiasi.c
Go to the documentation of this file.
1/*
2 This file is part of the IASI Code Collection.
3
4 the IASI Code Collections is free software: you can redistribute it
5 and/or modify it under the terms of the GNU General Public License
6 as published by the Free Software Foundation, either version 3 of
7 the License, or (at your option) any later version.
8
9 The IASI Code Collection is distributed in the hope that it will be
10 useful, but WITHOUT ANY WARRANTY; without even the implied warranty
11 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with the IASI Code Collection. If not, see
16 <http://www.gnu.org/licenses/>.
17
18 Copyright (C) 2019-2026 Forschungszentrum Juelich GmbH
19*/
20
26#include "libiasi.h"
27
28/*****************************************************************************/
29
30void add_var(
31 int ncid,
32 const char *varname,
33 const char *unit,
34 const char *longname,
35 int type,
36 int dimid[],
37 int *varid,
38 int ndims) {
39
40 /* Check if variable exists... */
41 if (nc_inq_varid(ncid, varname, varid) != NC_NOERR) {
42
43 /* Define variable... */
44 NC(nc_def_var(ncid, varname, type, ndims, dimid, varid));
45
46 /* Set long name... */
47 NC(nc_put_att_text
48 (ncid, *varid, "long_name", strlen(longname), longname));
49
50 /* Set units... */
51 NC(nc_put_att_text(ncid, *varid, "units", strlen(unit), unit));
52 }
53}
54
55/*****************************************************************************/
56
57void background_poly_help(
58 double *xx,
59 double *yy,
60 int n,
61 int dim) {
62
63 gsl_multifit_linear_workspace *work;
64 gsl_matrix *cov, *X;
65 gsl_vector *c, *x, *y;
66
67 double chisq, xx2[WX > WY ? WX : WY], yy2[WX > WY ? WX : WY];
68
69 size_t i, i2, n2 = 0;
70
71 /* Check for nan... */
72 for (i = 0; i < (size_t) n; i++)
73 if (gsl_finite(yy[i])) {
74 xx2[n2] = xx[i];
75 yy2[n2] = yy[i];
76 n2++;
77 }
78 if ((int) n2 < dim || n2 < 0.9 * n) {
79 for (i = 0; i < (size_t) n; i++)
80 yy[i] = GSL_NAN;
81 return;
82 }
83
84 /* Allocate... */
85 work = gsl_multifit_linear_alloc((size_t) n2, (size_t) dim);
86 cov = gsl_matrix_alloc((size_t) dim, (size_t) dim);
87 X = gsl_matrix_alloc((size_t) n2, (size_t) dim);
88 c = gsl_vector_alloc((size_t) dim);
89 x = gsl_vector_alloc((size_t) n2);
90 y = gsl_vector_alloc((size_t) n2);
91
92 /* Compute polynomial fit... */
93 for (i = 0; i < (size_t) n2; i++) {
94 gsl_vector_set(x, i, xx2[i]);
95 gsl_vector_set(y, i, yy2[i]);
96 for (i2 = 0; i2 < (size_t) dim; i2++)
97 gsl_matrix_set(X, i, i2, pow(gsl_vector_get(x, i), (double) i2));
98 }
99 gsl_multifit_linear(X, y, c, cov, &chisq, work);
100 for (i = 0; i < (size_t) n; i++)
101 yy[i] = gsl_poly_eval(c->data, (int) dim, xx[i]);
102
103 /* Free... */
104 gsl_multifit_linear_free(work);
105 gsl_matrix_free(cov);
106 gsl_matrix_free(X);
107 gsl_vector_free(c);
108 gsl_vector_free(x);
109 gsl_vector_free(y);
110}
111
112/*****************************************************************************/
113
114void background_poly(
115 wave_t *wave,
116 int dim_x,
117 int dim_y) {
118
119 double help[WX], x[WX], x2[WY], y[WX], y2[WY];
120
121 int ix, iy;
122
123 /* Copy temperatures to background... */
124 for (ix = 0; ix < wave->nx; ix++)
125 for (iy = 0; iy < wave->ny; iy++) {
126 wave->bg[ix][iy] = wave->temp[ix][iy];
127 wave->pt[ix][iy] = 0;
128 }
129
130 /* Check parameters... */
131 if (dim_x <= 0 && dim_y <= 0)
132 return;
133
134 /* Compute fit in x-direction... */
135 if (dim_x > 0)
136 for (iy = 0; iy < wave->ny; iy++) {
137 for (ix = 0; ix <= 53; ix++) {
138 x[ix] = (double) ix;
139 y[ix] = wave->bg[ix][iy];
140 }
141 background_poly_help(x, y, 54, dim_x);
142 for (ix = 0; ix <= 29; ix++)
143 help[ix] = y[ix];
144
145 for (ix = 6; ix <= 59; ix++) {
146 x[ix - 6] = (double) ix;
147 y[ix - 6] = wave->bg[ix][iy];
148 }
149 background_poly_help(x, y, 54, dim_x);
150 for (ix = 30; ix <= 59; ix++)
151 help[ix] = y[ix - 6];
152
153 for (ix = 0; ix < wave->nx; ix++)
154 wave->bg[ix][iy] = help[ix];
155 }
156
157 /* Compute fit in y-direction... */
158 if (dim_y > 0)
159 for (ix = 0; ix < wave->nx; ix++) {
160 for (iy = 0; iy < wave->ny; iy++) {
161 x2[iy] = (int) iy;
162 y2[iy] = wave->bg[ix][iy];
163 }
164 background_poly_help(x2, y2, wave->ny, dim_y);
165 for (iy = 0; iy < wave->ny; iy++)
166 wave->bg[ix][iy] = y2[iy];
167 }
168
169 /* Recompute perturbations... */
170 for (ix = 0; ix < wave->nx; ix++)
171 for (iy = 0; iy < wave->ny; iy++)
172 wave->pt[ix][iy] = wave->temp[ix][iy] - wave->bg[ix][iy];
173}
174
175/*****************************************************************************/
176
177void background_smooth(
178 wave_t *wave,
179 int npts_x,
180 int npts_y) {
181
182 const double dmax = 2500.0;
183
184 static double help[WX][WY];
185
186 /* Check parameters... */
187 if (npts_x <= 0 && npts_y <= 0)
188 return;
189
190 /* Smooth background... */
191 for (int ix = 0; ix < wave->nx; ix++)
192 for (int iy = 0; iy < wave->ny; iy++) {
193
194 /* Init... */
195 int n = 0;
196 help[ix][iy] = 0;
197
198 /* Set maximum range... */
199 const int dx = GSL_MIN(GSL_MIN(npts_x, ix), wave->nx - 1 - ix);
200 const int dy = GSL_MIN(GSL_MIN(npts_y, iy), wave->ny - 1 - iy);
201
202 /* Average... */
203 for (int i = ix - dx; i <= ix + dx; i++)
204 for (int j = iy - dy; j <= iy + dy; j++)
205 if (fabs(wave->x[ix] - wave->x[i]) < dmax &&
206 fabs(wave->y[iy] - wave->y[j]) < dmax) {
207 help[ix][iy] += wave->bg[i][j];
208 n++;
209 }
210
211 /* Normalize... */
212 if (n > 0)
213 help[ix][iy] /= n;
214 else
215 help[ix][iy] = GSL_NAN;
216 }
217
218 /* Recalculate perturbations... */
219 for (int ix = 0; ix < wave->nx; ix++)
220 for (int iy = 0; iy < wave->ny; iy++) {
221 wave->bg[ix][iy] = help[ix][iy];
222 wave->pt[ix][iy] = wave->temp[ix][iy] - wave->bg[ix][iy];
223 }
224}
225
226/*****************************************************************************/
227
228void gauss(
229 wave_t *wave,
230 double fwhm) {
231
232 static double help[WX][WY];
233
234 /* Check parameters... */
235 if (fwhm <= 0)
236 return;
237
238 /* Compute sigma^2... */
239 const double sigma2 = gsl_pow_2(fwhm / 2.3548);
240
241 /* Loop over data points... */
242 for (int ix = 0; ix < wave->nx; ix++)
243 for (int iy = 0; iy < wave->ny; iy++) {
244
245 /* Init... */
246 double wsum = 0;
247 help[ix][iy] = 0;
248
249 /* Average... */
250 for (int ix2 = 0; ix2 < wave->nx; ix2++)
251 for (int iy2 = 0; iy2 < wave->ny; iy2++) {
252 const double d2 = gsl_pow_2(wave->x[ix] - wave->x[ix2])
253 + gsl_pow_2(wave->y[iy] - wave->y[iy2]);
254 if (d2 <= 9 * sigma2) {
255 const double w = exp(-d2 / (2 * sigma2));
256 wsum += w;
257 help[ix][iy] += w * wave->pt[ix2][iy2];
258 }
259 }
260
261 /* Normalize... */
262 wave->pt[ix][iy] = help[ix][iy] / wsum;
263 }
264}
265
266/*****************************************************************************/
267
268void hamming(
269 wave_t *wave,
270 int niter) {
271
272 static double help[WX][WY];
273
274 /* Iterations... */
275 for (int iter = 0; iter < niter; iter++) {
276
277 /* Filter in x direction... */
278 for (int ix = 0; ix < wave->nx; ix++)
279 for (int iy = 0; iy < wave->ny; iy++)
280 help[ix][iy]
281 = 0.23 * wave->pt[ix > 0 ? ix - 1 : ix][iy]
282 + 0.54 * wave->pt[ix][iy]
283 + 0.23 * wave->pt[ix < wave->nx - 1 ? ix + 1 : ix][iy];
284
285 /* Filter in y direction... */
286 for (int ix = 0; ix < wave->nx; ix++)
287 for (int iy = 0; iy < wave->ny; iy++)
288 wave->pt[ix][iy]
289 = 0.23 * help[ix][iy > 0 ? iy - 1 : iy]
290 + 0.54 * help[ix][iy]
291 + 0.23 * help[ix][iy < wave->ny - 1 ? iy + 1 : iy];
292 }
293}
294
295/*****************************************************************************/
296
297void iasi_read(
298 int format,
299 char *filename,
300 iasi_rad_t *iasi_rad) {
301
302 /* Read native file... */
303 if (format == 1)
304 iasi_read_native(filename, iasi_rad);
305
306 /* Read netCDF file... */
307 else if (format == 2)
308 iasi_read_netcdf(filename, iasi_rad);
309
310 /* Error... */
311 else
312 ERRMSG("Unknown IASI Level-1 data format!");
313}
314
315/*****************************************************************************/
316
317void iasi_read_native(
318 char *filename,
319 iasi_rad_t *iasi_rad) {
320
321 const char *product_class;
322
323 coda_product *pf;
324
325 coda_cursor cursor;
326
327 iasi_raw_t *iasi_raw;
328
329 int i, j, w, tr1, tr2, tr1_lpm, tr1_rpm, tr2_lpm, tr2_rpm,
330 ichan, mdr_i, num_dims = 1;
331
332 long dim[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9 };
333
334 short int IDefScaleSondNbScale, IDefScaleSondNsfirst[10],
335 IDefScaleSondNslast[10], IDefScaleSondScaleFactor[10];
336
337 float sc, scaling[IASI_L1_NCHAN];
338
339 /* Initialize CODA... */
340 coda_init();
341
342 /* Allocate... */
343 ALLOC(iasi_raw, iasi_raw_t, 1);
344
345 /* Open IASI file... */
346 CODA(coda_open(filename, &pf));
347 CODA(coda_get_product_class(pf, &product_class));
348 CODA(coda_cursor_set_product(&cursor, pf));
349
350 /* Get scaling parameters... */
351 CODA(coda_cursor_goto_record_field_by_name(&cursor, "GIADR_ScaleFactors"));
352
353 CODA(coda_cursor_goto_record_field_by_name
354 (&cursor, "IDefScaleSondNbScale"));
355 CODA(coda_cursor_read_int16(&cursor, &IDefScaleSondNbScale));
356 CODA(coda_cursor_goto_parent(&cursor));
357
358 CODA(coda_cursor_goto_record_field_by_name
359 (&cursor, "IDefScaleSondNsfirst"));
360 CODA(coda_cursor_read_int16_array
361 (&cursor, IDefScaleSondNsfirst, coda_array_ordering_c));
362 CODA(coda_cursor_goto_parent(&cursor));
363
364 CODA(coda_cursor_goto_record_field_by_name(&cursor, "IDefScaleSondNslast"));
365 CODA(coda_cursor_read_int16_array
366 (&cursor, IDefScaleSondNslast, coda_array_ordering_c));
367 CODA(coda_cursor_goto_parent(&cursor));
368
369 CODA(coda_cursor_goto_record_field_by_name
370 (&cursor, "IDefScaleSondScaleFactor"));
371 CODA(coda_cursor_read_int16_array
372 (&cursor, IDefScaleSondScaleFactor, coda_array_ordering_c));
373
374 /* Compute scaling factors... */
375 for (ichan = 0; ichan < IASI_L1_NCHAN; ichan++)
376 scaling[ichan] = GSL_NAN;
377 for (i = 0; i < IDefScaleSondNbScale; i++) {
378 sc = (float) pow(10.0, -IDefScaleSondScaleFactor[i]);
379 for (ichan = IDefScaleSondNsfirst[i] - 1;
380 ichan < IDefScaleSondNslast[i]; ichan++) {
381 w = ichan - IASI_IDefNsfirst1b + 1;
382 if (w >= 0 && w < IASI_L1_NCHAN)
383 scaling[w] = sc;
384 }
385 }
386
387 /* Get number of tracks in record... */
388 CODA(coda_cursor_goto_root(&cursor));
389 CODA(coda_cursor_goto_record_field_by_name(&cursor, "MDR"));
390 CODA(coda_cursor_get_array_dim(&cursor, &num_dims, dim));
391 iasi_raw->ntrack = dim[0];
392
393 /* Read tracks one by one... */
394 for (mdr_i = 0; mdr_i < iasi_raw->ntrack; mdr_i++) {
395
396 /* Reset cursor position... */
397 CODA(coda_cursor_goto_root(&cursor));
398
399 /* Move cursor to radiation data... */
400 CODA(coda_cursor_goto_record_field_by_name(&cursor, "MDR"));
401 CODA(coda_cursor_goto_array_element_by_index(&cursor, mdr_i));
402 CODA(coda_cursor_goto_record_field_by_name(&cursor, "MDR"));
403 CODA(coda_cursor_goto_record_field_by_name(&cursor, "GS1cSpect"));
404 CODA(coda_cursor_read_int16_array
405 (&cursor, &iasi_raw->Radiation[mdr_i][0][0][0],
406 coda_array_ordering_c));
407
408 /* Read time... */
409 CODA(coda_cursor_goto_parent(&cursor));
410 CODA(coda_cursor_goto_record_field_by_name(&cursor, "OnboardUTC"));
411 CODA(coda_cursor_read_double_array
412 (&cursor, &iasi_raw->Time[mdr_i][0], coda_array_ordering_c));
413
414 /* Read coordinates... */
415 CODA(coda_cursor_goto_parent(&cursor));
416 CODA(coda_cursor_goto_record_field_by_name(&cursor, "GGeoSondLoc"));
417 CODA(coda_cursor_read_double_array
418 (&cursor, &iasi_raw->Loc[mdr_i][0][0][0], coda_array_ordering_c));
419
420 /* Read satellite altitude... */
421 CODA(coda_cursor_goto_parent(&cursor));
422 CODA(coda_cursor_goto_record_field_by_name(&cursor,
423 "EARTH_SATELLITE_DISTANCE"));
424 CODA(coda_cursor_read_uint32(&cursor, &iasi_raw->Sat_z[mdr_i]));
425
426 /* Read spectral range... */
427 iasi_raw->IDefSpectDWn1b[mdr_i] = IASI_IDefSpectDWn1b / 100.0;
428
429 CODA(coda_cursor_goto_parent(&cursor));
430 CODA(coda_cursor_goto_record_field_by_name(&cursor, "IDefNsfirst1b"));
431 CODA(coda_cursor_read_int32(&cursor, &iasi_raw->IDefNsfirst1b[mdr_i]));
432 if (iasi_raw->IDefNsfirst1b[mdr_i] != IASI_IDefNsfirst1b)
433 ERRMSG("Unexpected value for IDefNsfirst1b!");
434
435 CODA(coda_cursor_goto_parent(&cursor));
436 CODA(coda_cursor_goto_record_field_by_name(&cursor, "IDefNslast1b"));
437 CODA(coda_cursor_read_int32(&cursor, &iasi_raw->IDefNslast1b[mdr_i]));
438 if (iasi_raw->IDefNslast1b[mdr_i] != IASI_IDefNslast1b)
439 ERRMSG("Unexpected value for IDefNslast1b!");
440
441 /* Compute wavenumber... */
442 if (mdr_i == 0)
443 for (i = 0; i < IASI_L1_NCHAN; i++)
444 iasi_raw->Wavenumber[i] =
445 iasi_raw->IDefSpectDWn1b[mdr_i] *
446 (float) (iasi_raw->IDefNsfirst1b[mdr_i] + i - 1);
447 }
448
449 /* Close file... */
450 CODA(coda_close(pf));
451
452 /* Finalize CODA... */
453 coda_done();
454
455 /* Set number of tracks... */
456 iasi_rad->ntrack = (int) (iasi_raw->ntrack * 2);
457
458 /* Copy wavenumbers... */
459 for (ichan = 0; ichan < IASI_L1_NCHAN; ichan++)
460 iasi_rad->freq[ichan] = iasi_raw->Wavenumber[ichan];
461
462 /* Copy footprint data... */
463 for (mdr_i = 0; mdr_i < iasi_raw->ntrack; mdr_i++) {
464 tr1 = mdr_i * 2;
465 tr2 = mdr_i * 2 + 1;
466 tr1_lpm = 3;
467 tr1_rpm = 0;
468 tr2_lpm = 2;
469 tr2_rpm = 1;
470
471 /* Copy time (2x2 matrix has same measurement time)... */
472 for (i = 0; i < IASI_NXTRACK; i++) {
473 iasi_rad->Time[tr1][i * 2] = iasi_raw->Time[mdr_i][i];
474 iasi_rad->Time[tr1][i * 2 + 1] = iasi_raw->Time[mdr_i][i];
475 iasi_rad->Time[tr2][i * 2] = iasi_raw->Time[mdr_i][i];
476 iasi_rad->Time[tr2][i * 2 + 1] = iasi_raw->Time[mdr_i][i];
477 }
478
479 /* Copy location... */
480 for (i = 0; i < IASI_NXTRACK; i++) {
481 iasi_rad->Longitude[tr1][i * 2] = iasi_raw->Loc[mdr_i][i][tr1_lpm][0];
482 iasi_rad->Longitude[tr1][i * 2 + 1] =
483 iasi_raw->Loc[mdr_i][i][tr1_rpm][0];
484 iasi_rad->Latitude[tr1][i * 2] = iasi_raw->Loc[mdr_i][i][tr1_lpm][1];
485 iasi_rad->Latitude[tr1][i * 2 + 1] =
486 iasi_raw->Loc[mdr_i][i][tr1_rpm][1];
487
488 iasi_rad->Longitude[tr2][i * 2] = iasi_raw->Loc[mdr_i][i][tr2_lpm][0];
489 iasi_rad->Longitude[tr2][i * 2 + 1] =
490 iasi_raw->Loc[mdr_i][i][tr2_rpm][0];
491 iasi_rad->Latitude[tr2][i * 2] = iasi_raw->Loc[mdr_i][i][tr2_lpm][1];
492 iasi_rad->Latitude[tr2][i * 2 + 1] =
493 iasi_raw->Loc[mdr_i][i][tr2_rpm][1];
494 }
495
496 /* Copy satellite location (we only have one height value)... */
497 iasi_rad->Sat_lon[tr1] = iasi_rad->Longitude[tr1][28];
498 iasi_rad->Sat_lat[tr1] = iasi_rad->Latitude[tr1][28];
499 iasi_rad->Sat_lon[tr2] = iasi_rad->Longitude[tr2][28];
500 iasi_rad->Sat_lat[tr2] = iasi_rad->Latitude[tr2][28];
501 iasi_rad->Sat_z[tr1] =
502 iasi_raw->Sat_z[mdr_i] / 1000.0 - wgs84(iasi_rad->Sat_lat[tr1]);
503 iasi_rad->Sat_z[tr2] =
504 iasi_raw->Sat_z[mdr_i] / 1000.0 - wgs84(iasi_rad->Sat_lat[tr2]);
505
506 /* Copy radiation data... */
507 for (i = 0; i < IASI_NXTRACK; i++) {
508 for (ichan = 0; ichan < IASI_L1_NCHAN; ichan++) {
509 sc = scaling[ichan] * 100.0f;
510 iasi_rad->Rad[tr1][i * 2][ichan] =
511 iasi_raw->Radiation[mdr_i][i][tr1_lpm][ichan] * sc;
512 iasi_rad->Rad[tr1][i * 2 + 1][ichan] =
513 iasi_raw->Radiation[mdr_i][i][tr1_rpm][ichan] * sc;
514 iasi_rad->Rad[tr2][i * 2][ichan] =
515 iasi_raw->Radiation[mdr_i][i][tr2_lpm][ichan] * sc;
516 iasi_rad->Rad[tr2][i * 2 + 1][ichan] =
517 iasi_raw->Radiation[mdr_i][i][tr2_rpm][ichan] * sc;
518 }
519 }
520 }
521
522 /* Check radiance data... */
523 for (i = 0; i < iasi_rad->ntrack; i++)
524 for (j = 0; j < L1_NXTRACK; j++)
525 if (iasi_rad->Rad[i][j][6753] > iasi_rad->Rad[i][j][6757]
526 || iasi_rad->Rad[i][j][6753] < 0)
527 for (ichan = 0; ichan < IASI_L1_NCHAN; ichan++)
528 iasi_rad->Rad[i][j][ichan] = GSL_NAN;
529
530 /* Free... */
531 free(iasi_raw);
532}
533
534/*****************************************************************************/
535
536typedef struct {
537 int orbit;
538 int scan;
539} orbit_scan_t;
540
541static int cmp_orbit_scan(
542 const void *a,
543 const void *b) {
544 const orbit_scan_t *x = (const orbit_scan_t *) a;
545 const orbit_scan_t *y = (const orbit_scan_t *) b;
546 if (x->orbit < y->orbit)
547 return -1;
548 if (x->orbit > y->orbit)
549 return 1;
550 if (x->scan < y->scan)
551 return -1;
552 if (x->scan > y->scan)
553 return 1;
554 return 0;
555}
556
557void iasi_read_netcdf(
558 char *filename,
559 iasi_rad_t *iasi_rad) {
560
561 int ncid = -1, dim_point_id = -1, dim_chan_id = -1, var_lat = -1,
562 var_lon = -1, var_date = -1, var_orbit = -1, var_scan = -1,
563 var_pixel = -1, var_fov = -1, var_channame = -1, var_R = -1,
564 *channel_name = NULL;
565
566 size_t npoint = 0, nchan = 0;
567
568 /* Full arrays... */
569 int *orbit_all = NULL;
570 int *scan_all = NULL;
571 int *pixel_all = NULL;
572 int *fov_all = NULL;
573 double *lat_all = NULL;
574 double *lon_all = NULL;
575 double *date_all = NULL;
576 float *R_all = NULL;
577
578 /* Scanline keys... */
579 orbit_scan_t *keys_all = NULL; /* length npoint */
580 orbit_scan_t *keys_uniq = NULL; /* length <= npoint */
581 size_t nuniq = 0;
582
583 /* Always leave struct safe... */
584 iasi_rad->ntrack = 0;
585
586 /* Standard IASI wavenumber grid... */
587 for (int ichan = 0; ichan < IASI_L1_NCHAN; ichan++)
588 iasi_rad->freq[ichan] = 644.75 + (double) (ichan + 1) * 0.25;
589
590 /* Open file... */
591 if (nc_open(filename, NC_NOWRITE, &ncid) != NC_NOERR) {
592 WARN("Cannot open file %s", filename);
593 return;
594 }
595
596 /* Get dimensions... */
597 NC(nc_inq_dimid(ncid, "point", &dim_point_id));
598 NC(nc_inq_dimlen(ncid, dim_point_id, &npoint));
599 NC(nc_inq_dimid(ncid, "channel", &dim_chan_id));
600 NC(nc_inq_dimlen(ncid, dim_chan_id, &nchan));
601 if (npoint == 0 || nchan == 0)
602 ERRMSG("Empty file, npoint and/or nchan is zero!");
603
604 /* Get variables... */
605 NC(nc_inq_varid(ncid, "lat", &var_lat));
606 NC(nc_inq_varid(ncid, "lon", &var_lon));
607 NC(nc_inq_varid(ncid, "date", &var_date));
608 NC(nc_inq_varid(ncid, "orbit", &var_orbit));
609 NC(nc_inq_varid(ncid, "scan", &var_scan));
610 NC(nc_inq_varid(ncid, "pixel", &var_pixel));
611 NC(nc_inq_varid(ncid, "fov", &var_fov));
612 NC(nc_inq_varid(ncid, "channel_name", &var_channame));
613 NC(nc_inq_varid(ncid, "R", &var_R));
614
615 /* Channel mapping... */
616 ALLOC(channel_name, int,
617 nchan);
618 NC(nc_get_var_int(ncid, var_channame, channel_name));
619
620 /* Allocate full arrays... */
621 ALLOC(orbit_all, int,
622 npoint);
623 ALLOC(scan_all, int,
624 npoint);
625 ALLOC(pixel_all, int,
626 npoint);
627 ALLOC(fov_all, int,
628 npoint);
629 ALLOC(lat_all, double,
630 npoint);
631 ALLOC(lon_all, double,
632 npoint);
633 ALLOC(date_all, double,
634 npoint);
635 ALLOC(R_all, float,
636 npoint * nchan);
637
638 /* Read everything... */
639 NC(nc_get_var_int(ncid, var_orbit, orbit_all));
640 NC(nc_get_var_int(ncid, var_scan, scan_all));
641 NC(nc_get_var_int(ncid, var_pixel, pixel_all));
642 NC(nc_get_var_int(ncid, var_fov, fov_all));
643 NC(nc_get_var_double(ncid, var_lat, lat_all));
644 NC(nc_get_var_double(ncid, var_lon, lon_all));
645 NC(nc_get_var_double(ncid, var_date, date_all));
646 NC(nc_get_var_float(ncid, var_R, R_all));
647
648 /* Build list of (orbit,scan) keys for each point... */
649 ALLOC(keys_all, orbit_scan_t, npoint);
650 for (size_t i = 0; i < npoint; i++) {
651 keys_all[i].orbit = orbit_all[i];
652 keys_all[i].scan = scan_all[i];
653 }
654
655 /* Sort keys_all, then unique into keys_uniq.. */
656 qsort(keys_all, npoint, sizeof(orbit_scan_t), cmp_orbit_scan);
657 ALLOC(keys_uniq, orbit_scan_t, npoint);
658 nuniq = 0;
659 for (size_t i = 0; i < npoint; i++)
660 if (i == 0 || cmp_orbit_scan(&keys_all[i], &keys_all[i - 1]) != 0)
661 keys_uniq[nuniq++] = keys_all[i];
662
663 /* Early check: file must fit completely... */
664 if (nuniq > (size_t) (L1_NTRACK / 2))
665 ERRMSG("Too many scanlines in file: %zu (max %d). Increase L1_NTRACK!",
666 nuniq, L1_NTRACK / 2);
667
668 /* ntrack = 2 * number of unique scanlines */
669 iasi_rad->ntrack = (int) (2 * nuniq);
670
671 /* Initialize radiances with NaN (for missing pixels)... */
672 for (int tr = 0; tr < iasi_rad->ntrack; tr++)
673 for (int ix = 0; ix < L1_NXTRACK; ix++)
674 for (int g = 0; g < IASI_L1_NCHAN; g++)
675 iasi_rad->Rad[tr][ix][g] = GSL_NAN;
676
677 /* No satellite position info... */
678 for (int tr = 0; tr < iasi_rad->ntrack; tr++) {
679 iasi_rad->Sat_z[tr] = GSL_NAN;
680 iasi_rad->Sat_lon[tr] = GSL_NAN;
681 iasi_rad->Sat_lat[tr] = GSL_NAN;
682 }
683
684 /* Fill iasi_rad... */
685 for (size_t i = 0; i < npoint; i++) {
686
687 /* Find scanline index from (orbit,scan)... */
688 orbit_scan_t key;
689 key.orbit = orbit_all[i];
690 key.scan = scan_all[i];
691 orbit_scan_t *found =
692 (orbit_scan_t *) bsearch(&key, keys_uniq, nuniq, sizeof(orbit_scan_t),
693 cmp_orbit_scan);
694 if (!found)
695 continue;
696 size_t scan_idx = (size_t) (found - keys_uniq);
697
698 int f = fov_all[i];
699 if (f < 1 || f > 120)
700 continue;
701
702 /* Get position... */
703 int pos = (f - 1) / 4; /* 0..29 */
704 int sub = (f - 1) % 4; /* 0..3 */
705 if (pos < 0 || pos >= IASI_NXTRACK)
706 continue;
707
708 int track_add, x_add;
709 switch (sub) {
710 case 3:
711 track_add = 0;
712 x_add = 0;
713 break; /* tr1_lpm */
714 case 0:
715 track_add = 0;
716 x_add = 1;
717 break; /* tr1_rpm */
718 case 2:
719 track_add = 1;
720 x_add = 0;
721 break; /* tr2_lpm */
722 case 1:
723 track_add = 1;
724 x_add = 1;
725 break; /* tr2_rpm */
726 default:
727 continue;
728 }
729
730 int tr = 2 * (int) scan_idx + track_add;
731 int ix = 2 * pos + x_add;
732
733 if (tr < 0 || tr >= iasi_rad->ntrack)
734 continue;
735 if (ix < 0 || ix >= L1_NXTRACK)
736 continue;
737
738 /* Save data... */
739 iasi_rad->Time[tr][ix] = date_all[i];
740 iasi_rad->Longitude[tr][ix] = lon_all[i];
741 iasi_rad->Latitude[tr][ix] = lat_all[i];
742 float *row = &R_all[i * nchan];
743 for (size_t k = 0; k < nchan; k++) {
744 int g = channel_name[k] - 1;
745 if (g >= 0 && g < IASI_L1_NCHAN)
746 iasi_rad->Rad[tr][ix][g] = 100.0f * row[k]; /* m^-1 -> cm^-1 */
747 }
748 }
749
750 /* Close file... */
751 NC(nc_close(ncid));
752
753 /* Free... */
754 free(channel_name);
755 free(orbit_all);
756 free(scan_all);
757 free(pixel_all);
758 free(fov_all);
759 free(lat_all);
760 free(lon_all);
761 free(date_all);
762 free(R_all);
763 free(keys_all);
764 free(keys_uniq);
765}
766
767/*****************************************************************************/
768
769void median(
770 wave_t *wave,
771 int dx) {
772
773 static double data[WX * WY], help[WX][WY];
774
775 /* Check parameters... */
776 if (dx <= 0)
777 return;
778
779 /* Loop over data points... */
780 for (int ix = 0; ix < wave->nx; ix++)
781 for (int iy = 0; iy < wave->ny; iy++) {
782
783 /* Init... */
784 size_t n = 0;
785
786 /* Get data... */
787 for (int ix2 = GSL_MAX(ix - dx, 0);
788 ix2 < GSL_MIN(ix + dx, wave->nx - 1); ix2++)
789 for (int iy2 = GSL_MAX(iy - dx, 0);
790 iy2 < GSL_MIN(iy + dx, wave->ny - 1); iy2++) {
791 data[n] = wave->pt[ix2][iy2];
792 n++;
793 }
794
795 /* Normalize... */
796 gsl_sort(data, 1, n);
797 help[ix][iy] = gsl_stats_median_from_sorted_data(data, 1, n);
798 }
799
800 /* Loop over data points... */
801 for (int ix = 0; ix < wave->nx; ix++)
802 for (int iy = 0; iy < wave->ny; iy++)
803 wave->pt[ix][iy] = help[ix][iy];
804}
805
806/*****************************************************************************/
807
808void noise(
809 wave_t *wave,
810 double *mu,
811 double *sig) {
812
813 int ix, ix2, iy, iy2, n = 0, okay;
814
815 /* Init... */
816 *mu = 0;
817 *sig = 0;
818
819 /* Estimate noise (Immerkaer, 1996)... */
820 for (ix = 1; ix < wave->nx - 1; ix++)
821 for (iy = 1; iy < wave->ny - 1; iy++) {
822
823 /* Check data... */
824 okay = 1;
825 for (ix2 = ix - 1; ix2 <= ix + 1; ix2++)
826 for (iy2 = iy - 1; iy2 <= iy + 1; iy2++)
827 if (!gsl_finite(wave->temp[ix2][iy2]))
828 okay = 0;
829 if (!okay)
830 continue;
831
832 /* Get mean noise... */
833 n++;
834 *mu += wave->temp[ix][iy];
835 *sig += gsl_pow_2(+4. / 6. * wave->temp[ix][iy]
836 - 2. / 6. * (wave->temp[ix - 1][iy]
837 + wave->temp[ix + 1][iy]
838 + wave->temp[ix][iy - 1]
839 + wave->temp[ix][iy + 1])
840 + 1. / 6. * (wave->temp[ix - 1][iy - 1]
841 + wave->temp[ix + 1][iy - 1]
842 + wave->temp[ix - 1][iy + 1]
843 + wave->temp[ix + 1][iy + 1]));
844 }
845
846 /* Normalize... */
847 *mu /= (double) n;
848 *sig = sqrt(*sig / (double) n);
849}
850
851/*****************************************************************************/
852
853void pert2wave(
854 pert_t *pert,
855 wave_t *wave,
856 int track0,
857 int track1,
858 int xtrack0,
859 int xtrack1) {
860
861 double x0[3], x1[3];
862
863 int itrack, ixtrack;
864
865 /* Check ranges... */
866 track0 = GSL_MIN(GSL_MAX(track0, 0), pert->ntrack - 1);
867 track1 = GSL_MIN(GSL_MAX(track1, 0), pert->ntrack - 1);
868 xtrack0 = GSL_MIN(GSL_MAX(xtrack0, 0), pert->nxtrack - 1);
869 xtrack1 = GSL_MIN(GSL_MAX(xtrack1, 0), pert->nxtrack - 1);
870
871 /* Set size... */
872 wave->nx = xtrack1 - xtrack0 + 1;
873 if (wave->nx > WX)
874 ERRMSG("Too many across-track values!");
875 wave->ny = track1 - track0 + 1;
876 if (wave->ny > WY)
877 ERRMSG("Too many along-track values!");
878
879 /* Loop over footprints... */
880 for (itrack = track0; itrack <= track1; itrack++)
881 for (ixtrack = xtrack0; ixtrack <= xtrack1; ixtrack++) {
882
883 /* Get distances... */
884 if (itrack == track0) {
885 wave->x[0] = 0;
886 if (ixtrack > xtrack0) {
887 geo2cart(0, pert->lon[itrack][ixtrack - 1],
888 pert->lat[itrack][ixtrack - 1], x0);
889 geo2cart(0, pert->lon[itrack][ixtrack],
890 pert->lat[itrack][ixtrack], x1);
891 wave->x[ixtrack - xtrack0] =
892 wave->x[ixtrack - xtrack0 - 1] + DIST(x0, x1);
893 }
894 }
895 if (ixtrack == xtrack0) {
896 wave->y[0] = 0;
897 if (itrack > track0) {
898 geo2cart(0, pert->lon[itrack - 1][ixtrack],
899 pert->lat[itrack - 1][ixtrack], x0);
900 geo2cart(0, pert->lon[itrack][ixtrack],
901 pert->lat[itrack][ixtrack], x1);
902 wave->y[itrack - track0] =
903 wave->y[itrack - track0 - 1] + DIST(x0, x1);
904 }
905 }
906
907 /* Save geolocation... */
908 wave->time = pert->time[(track0 + track1) / 2][(xtrack0 + xtrack1) / 2];
909 wave->z = 0;
910 wave->lon[ixtrack - xtrack0][itrack - track0] =
911 pert->lon[itrack][ixtrack];
912 wave->lat[ixtrack - xtrack0][itrack - track0] =
913 pert->lat[itrack][ixtrack];
914
915 /* Save temperature data... */
916 wave->temp[ixtrack - xtrack0][itrack - track0]
917 = pert->bt[itrack][ixtrack];
918 wave->bg[ixtrack - xtrack0][itrack - track0]
919 = pert->bt[itrack][ixtrack] - pert->pt[itrack][ixtrack];
920 wave->pt[ixtrack - xtrack0][itrack - track0]
921 = pert->pt[itrack][ixtrack];
922 wave->var[ixtrack - xtrack0][itrack - track0]
923 = pert->var[itrack][ixtrack];
924 }
925}
926
927/*****************************************************************************/
928
929void read_pert(
930 char *filename,
931 char *pertname,
932 pert_t *pert) {
933
934 static char varname[LEN];
935
936 static int dimid[2], ncid, varid;
937
938 static size_t ntrack, nxtrack, start[2] = { 0, 0 }, count[2] = { 1, 1 };
939
940 /* Write info... */
941 printf("Read perturbation data: %s\n", filename);
942
943 /* Open netCDF file... */
944 NC(nc_open(filename, NC_NOWRITE, &ncid));
945
946 /* Get dimensions... */
947 NC(nc_inq_dimid(ncid, "NTRACK", &dimid[0]));
948 NC(nc_inq_dimid(ncid, "NXTRACK", &dimid[1]));
949 NC(nc_inq_dimlen(ncid, dimid[0], &ntrack));
950 NC(nc_inq_dimlen(ncid, dimid[1], &nxtrack));
951 if (nxtrack > PERT_NXTRACK)
952 ERRMSG("Too many tracks!");
953 if (ntrack > PERT_NTRACK)
954 ERRMSG("Too many scans!");
955 pert->ntrack = (int) ntrack;
956 pert->nxtrack = (int) nxtrack;
957 count[1] = nxtrack;
958
959 /* Read data... */
960 NC(nc_inq_varid(ncid, "time", &varid));
961 for (size_t itrack = 0; itrack < ntrack; itrack++) {
962 start[0] = itrack;
963 NC(nc_get_vara_double(ncid, varid, start, count, pert->time[itrack]));
964 }
965
966 NC(nc_inq_varid(ncid, "lon", &varid));
967 for (size_t itrack = 0; itrack < ntrack; itrack++) {
968 start[0] = itrack;
969 NC(nc_get_vara_double(ncid, varid, start, count, pert->lon[itrack]));
970 }
971
972 NC(nc_inq_varid(ncid, "lat", &varid));
973 for (size_t itrack = 0; itrack < ntrack; itrack++) {
974 start[0] = itrack;
975 NC(nc_get_vara_double(ncid, varid, start, count, pert->lat[itrack]));
976 }
977
978 NC(nc_inq_varid(ncid, "bt_8mu", &varid));
979 for (size_t itrack = 0; itrack < ntrack; itrack++) {
980 start[0] = itrack;
981 NC(nc_get_vara_double(ncid, varid, start, count, pert->dc[itrack]));
982 }
983
984 sprintf(varname, "bt_%s", pertname);
985 NC(nc_inq_varid(ncid, varname, &varid));
986 for (size_t itrack = 0; itrack < ntrack; itrack++) {
987 start[0] = itrack;
988 NC(nc_get_vara_double(ncid, varid, start, count, pert->bt[itrack]));
989 }
990
991 sprintf(varname, "bt_%s_pt", pertname);
992 NC(nc_inq_varid(ncid, varname, &varid));
993 for (size_t itrack = 0; itrack < ntrack; itrack++) {
994 start[0] = itrack;
995 NC(nc_get_vara_double(ncid, varid, start, count, pert->pt[itrack]));
996 }
997
998 sprintf(varname, "bt_%s_var", pertname);
999 NC(nc_inq_varid(ncid, varname, &varid));
1000 for (size_t itrack = 0; itrack < ntrack; itrack++) {
1001 start[0] = itrack;
1002 NC(nc_get_vara_double(ncid, varid, start, count, pert->var[itrack]));
1003 }
1004
1005 /* Close file... */
1006 NC(nc_close(ncid));
1007}
1008
1009/*****************************************************************************/
1010
1011void variance(
1012 wave_t *wave,
1013 double dh) {
1014
1015 double dh2, mu, help;
1016
1017 int dx, dy, ix, ix2, iy, iy2, n;
1018
1019 /* Check parameters... */
1020 if (dh <= 0)
1021 return;
1022
1023 /* Compute squared radius... */
1024 dh2 = gsl_pow_2(dh);
1025
1026 /* Get sampling distances... */
1027 dx =
1028 (int) (dh / fabs(wave->x[wave->nx - 1] - wave->x[0]) *
1029 (wave->nx - 1.0) + 1);
1030 dy =
1031 (int) (dh / fabs(wave->y[wave->ny - 1] - wave->y[0]) *
1032 (wave->ny - 1.0) + 1);
1033
1034 /* Loop over data points... */
1035 for (ix = 0; ix < wave->nx; ix++)
1036 for (iy = 0; iy < wave->ny; iy++) {
1037
1038 /* Init... */
1039 mu = help = 0;
1040 n = 0;
1041
1042 /* Get data... */
1043 for (ix2 = GSL_MAX(ix - dx, 0); ix2 <= GSL_MIN(ix + dx, wave->nx - 1);
1044 ix2++)
1045 for (iy2 = GSL_MAX(iy - dy, 0);
1046 iy2 <= GSL_MIN(iy + dy, wave->ny - 1); iy2++)
1047 if ((gsl_pow_2(wave->x[ix] - wave->x[ix2])
1048 + gsl_pow_2(wave->y[iy] - wave->y[iy2])) <= dh2)
1049 if (gsl_finite(wave->pt[ix2][iy2])) {
1050 mu += wave->pt[ix2][iy2];
1051 help += gsl_pow_2(wave->pt[ix2][iy2]);
1052 n++;
1053 }
1054
1055 /* Compute local variance... */
1056 if (n > 1)
1057 wave->var[ix][iy] = help / n - gsl_pow_2(mu / n);
1058 else
1059 wave->var[ix][iy] = GSL_NAN;
1060 }
1061}
1062
1063/*****************************************************************************/
1064
1065double wgs84(
1066 double lat) {
1067
1068 const double a = 6378.1370, b = 6356.7523;
1069
1070 double cphi, sphi;
1071
1072 cphi = cos(lat * M_PI / 180.);
1073 sphi = sin(lat * M_PI / 180.);
1074
1075 return sqrt((gsl_pow_2(a * a * cphi) + gsl_pow_2(b * b * sphi))
1076 / (gsl_pow_2(a * cphi) + gsl_pow_2(b * sphi)));
1077}
1078
1079/*****************************************************************************/
1080
1081void write_l1(
1082 char *filename,
1083 iasi_l1_t *l1) {
1084
1085 int dimid[10], ncid, time_id, lon_id, lat_id,
1086 sat_z_id, sat_lon_id, sat_lat_id, nu_id, rad_id;
1087
1088 /* Open or create netCDF file... */
1089 printf("Write IASI Level-1 file: %s\n", filename);
1090 if (nc_open(filename, NC_WRITE, &ncid) != NC_NOERR) {
1091 NC(nc_create(filename, NC_CLOBBER, &ncid));
1092 } else {
1093 NC(nc_redef(ncid));
1094 }
1095
1096 /* Set dimensions... */
1097 if (nc_inq_dimid(ncid, "L1_NTRACK", &dimid[0]) != NC_NOERR)
1098 NC(nc_def_dim(ncid, "L1_NTRACK", l1->ntrack, &dimid[0]));
1099 if (nc_inq_dimid(ncid, "L1_NXTRACK", &dimid[1]) != NC_NOERR)
1100 NC(nc_def_dim(ncid, "L1_NXTRACK", L1_NXTRACK, &dimid[1]));
1101 if (nc_inq_dimid(ncid, "L1_NCHAN", &dimid[2]) != NC_NOERR)
1102 NC(nc_def_dim(ncid, "L1_NCHAN", L1_NCHAN, &dimid[2]));
1103
1104 /* Add variables... */
1105 add_var(ncid, "l1_time", "s", "time (seconds since 2000-01-01T00:00Z)",
1106 NC_DOUBLE, dimid, &time_id, 2);
1107 add_var(ncid, "l1_lon", "deg", "longitude", NC_DOUBLE, dimid, &lon_id, 2);
1108 add_var(ncid, "l1_lat", "deg", "latitude", NC_DOUBLE, dimid, &lat_id, 2);
1109 add_var(ncid, "l1_sat_z", "km", "satellite altitude",
1110 NC_DOUBLE, dimid, &sat_z_id, 1);
1111 add_var(ncid, "l1_sat_lon", "deg", "(estimated) satellite longitude",
1112 NC_DOUBLE, dimid, &sat_lon_id, 1);
1113 add_var(ncid, "l1_sat_lat", "deg", "(estimated) satellite latitude",
1114 NC_DOUBLE, dimid, &sat_lat_id, 1);
1115 add_var(ncid, "l1_nu", "cm^-1", "channel wavenumber",
1116 NC_DOUBLE, &dimid[2], &nu_id, 1);
1117 add_var(ncid, "l1_rad", "W/(m^2 sr cm^-1)", "channel radiance",
1118 NC_FLOAT, dimid, &rad_id, 3);
1119
1120 /* Leave define mode... */
1121 NC(nc_enddef(ncid));
1122
1123 /* Write data... */
1124 NC(nc_put_var_double(ncid, time_id, l1->time[0]));
1125 NC(nc_put_var_double(ncid, lon_id, l1->lon[0]));
1126 NC(nc_put_var_double(ncid, lat_id, l1->lat[0]));
1127 NC(nc_put_var_double(ncid, sat_z_id, l1->sat_z));
1128 NC(nc_put_var_double(ncid, sat_lon_id, l1->sat_lon));
1129 NC(nc_put_var_double(ncid, sat_lat_id, l1->sat_lat));
1130 NC(nc_put_var_double(ncid, nu_id, l1->nu));
1131 NC(nc_put_var_float(ncid, rad_id, &l1->rad[0][0][0]));
1132
1133 /* Close file... */
1134 NC(nc_close(ncid));
1135}
1136
1137/*****************************************************************************/
1138
1139void write_l2(
1140 char *filename,
1141 iasi_l2_t *l2) {
1142
1143 int dimid[10], ncid, time_id, z_id, lon_id, lat_id, p_id, t_id;
1144
1145 /* Create netCDF file... */
1146 printf("Write IASI Level-2 file: %s\n", filename);
1147 if (nc_open(filename, NC_WRITE, &ncid) != NC_NOERR) {
1148 NC(nc_create(filename, NC_CLOBBER, &ncid));
1149 } else {
1150 NC(nc_redef(ncid));
1151 }
1152
1153 /* Set dimensions... */
1154 if (nc_inq_dimid(ncid, "L2_NTRACK", &dimid[0]) != NC_NOERR)
1155 NC(nc_def_dim(ncid, "L2_NTRACK", l2->ntrack, &dimid[0]));
1156 if (nc_inq_dimid(ncid, "L2_NXTRACK", &dimid[1]) != NC_NOERR)
1157 NC(nc_def_dim(ncid, "L2_NXTRACK", L2_NXTRACK, &dimid[1]));
1158 if (nc_inq_dimid(ncid, "L2_NLAY", &dimid[2]) != NC_NOERR)
1159 NC(nc_def_dim(ncid, "L2_NLAY", L2_NLAY, &dimid[2]));
1160
1161 /* Add variables... */
1162 add_var(ncid, "l2_time", "s", "time (seconds since 2000-01-01T00:00Z)",
1163 NC_DOUBLE, dimid, &time_id, 2);
1164 add_var(ncid, "l2_z", "km", "altitude", NC_DOUBLE, dimid, &z_id, 3);
1165 add_var(ncid, "l2_lon", "deg", "longitude", NC_DOUBLE, dimid, &lon_id, 2);
1166 add_var(ncid, "l2_lat", "deg", "latitude", NC_DOUBLE, dimid, &lat_id, 2);
1167 add_var(ncid, "l2_press", "hPa", "pressure",
1168 NC_DOUBLE, &dimid[2], &p_id, 1);
1169 add_var(ncid, "l2_temp", "K", "temperature", NC_DOUBLE, dimid, &t_id, 3);
1170
1171 /* Leave define mode... */
1172 NC(nc_enddef(ncid));
1173
1174 /* Write data... */
1175 NC(nc_put_var_double(ncid, time_id, l2->time[0]));
1176 NC(nc_put_var_double(ncid, z_id, l2->z[0][0]));
1177 NC(nc_put_var_double(ncid, lon_id, l2->lon[0]));
1178 NC(nc_put_var_double(ncid, lat_id, l2->lat[0]));
1179 NC(nc_put_var_double(ncid, p_id, l2->p));
1180 NC(nc_put_var_double(ncid, t_id, l2->t[0][0]));
1181
1182 /* Close file... */
1183 NC(nc_close(ncid));
1184}
iasi_read
void iasi_read(int format, char *filename, iasi_rad_t *iasi_rad)
Read IASI Level-1 data.
Definition: libiasi.c:297
background_poly_help
void background_poly_help(double *xx, double *yy, int n, int dim)
Get background based on polynomial fits.
Definition: libiasi.c:57
background_smooth
void background_smooth(wave_t *wave, int npts_x, int npts_y)
Smooth background.
Definition: libiasi.c:177
write_l2
void write_l2(char *filename, iasi_l2_t *l2)
Write IASI Level-2 data.
Definition: libiasi.c:1139
noise
void noise(wave_t *wave, double *mu, double *sig)
Estimate noise.
Definition: libiasi.c:808
write_l1
void write_l1(char *filename, iasi_l1_t *l1)
Write IASI Level-1 data.
Definition: libiasi.c:1081
iasi_read_netcdf
void iasi_read_netcdf(char *filename, iasi_rad_t *iasi_rad)
Read IASI Level-1 data from netCDF file.
Definition: libiasi.c:557
hamming
void hamming(wave_t *wave, int niter)
Apply Hamming filter to perturbations...
Definition: libiasi.c:268
background_poly
void background_poly(wave_t *wave, int dim_x, int dim_y)
Get background based on polynomial fits.
Definition: libiasi.c:114
variance
void variance(wave_t *wave, double dh)
Compute local variance.
Definition: libiasi.c:1011
add_var
void add_var(int ncid, const char *varname, const char *unit, const char *longname, int type, int dimid[], int *varid, int ndims)
Add variable to netCDF file.
Definition: libiasi.c:30
pert2wave
void pert2wave(pert_t *pert, wave_t *wave, int track0, int track1, int xtrack0, int xtrack1)
Convert radiance perturbation data to wave analysis struct.
Definition: libiasi.c:853
read_pert
void read_pert(char *filename, char *pertname, pert_t *pert)
Read radiance perturbation data.
Definition: libiasi.c:929
gauss
void gauss(wave_t *wave, double fwhm)
Apply Gaussian filter to perturbations...
Definition: libiasi.c:228
iasi_read_native
void iasi_read_native(char *filename, iasi_rad_t *iasi_rad)
Read IASI Level-1 data from native file.
Definition: libiasi.c:317
median
void median(wave_t *wave, int dx)
Apply median filter to perturbations...
Definition: libiasi.c:769
wgs84
double wgs84(double lat)
Calculate Earth radius according to WGS-84 reference ellipsoid.
Definition: libiasi.c:1065
libiasi.h
IASI Code Collection library declarations.
IASI_NXTRACK
#define IASI_NXTRACK
Raw data across-track size of IASI radiance granule.
Definition: libiasi.h:117
IASI_IDefNslast1b
#define IASI_IDefNslast1b
Expected value for the computation of the last wavenumber.
Definition: libiasi.h:126
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: libiasi.h:155
PERT_NXTRACK
#define PERT_NXTRACK
Across-track size of perturbation data.
Definition: libiasi.h:135
L1_NXTRACK
#define L1_NXTRACK
Across-track size of IASI radiance granule (don't change).
Definition: libiasi.h:102
L1_NTRACK
#define L1_NTRACK
Maximum along-track size of IASI radiance granule (don't change).
Definition: libiasi.h:99
IASI_L1_NCHAN
#define IASI_L1_NCHAN
Number of channels of IASI radiance granule.
Definition: libiasi.h:114
WX
#define WX
Across-track size of wave analysis data.
Definition: libiasi.h:138
L2_NXTRACK
#define L2_NXTRACK
Across-track size of IASI retrieval granule (don't change).
Definition: libiasi.h:111
PERT_NTRACK
#define PERT_NTRACK
Along-track size of perturbation data.
Definition: libiasi.h:132
IASI_IDefSpectDWn1b
#define IASI_IDefSpectDWn1b
Expected value for the interval of the IASI wavenumbers [m^-1].
Definition: libiasi.h:129
L2_NLAY
#define L2_NLAY
Number of IASI pressure layers (don't change).
Definition: libiasi.h:105
IASI_IDefNsfirst1b
#define IASI_IDefNsfirst1b
Expected value for the computation of the first wavenumber.
Definition: libiasi.h:123
L1_NCHAN
#define L1_NCHAN
Number of IASI radiance channels (don't change).
Definition: libiasi.h:96
CODA
#define CODA(cmd)
Execute CODA library command and check result.
Definition: libiasi.h:148
WY
#define WY
Along-track size of wave analysis data.
Definition: libiasi.h:141
iasi_l1_t
IASI Level-1 data.
Definition: libiasi.h:166
iasi_l1_t::sat_z
double sat_z[L1_NTRACK]
Satellite altitude [km].
Definition: libiasi.h:181
iasi_l1_t::nu
double nu[L1_NCHAN]
Channel frequencies [cm^-1].
Definition: libiasi.h:190
iasi_l1_t::sat_lon
double sat_lon[L1_NTRACK]
Satellite longitude [deg].
Definition: libiasi.h:184
iasi_l1_t::sat_lat
double sat_lat[L1_NTRACK]
Satellite latitude [deg].
Definition: libiasi.h:187
iasi_l1_t::lon
double lon[L1_NTRACK][L1_NXTRACK]
Footprint longitude [deg].
Definition: libiasi.h:175
iasi_l1_t::time
double time[L1_NTRACK][L1_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libiasi.h:172
iasi_l1_t::lat
double lat[L1_NTRACK][L1_NXTRACK]
Footprint latitude [deg].
Definition: libiasi.h:178
iasi_l1_t::rad
float rad[L1_NTRACK][L1_NXTRACK][L1_NCHAN]
Radiance [W/(m^2 sr cm^-1)].
Definition: libiasi.h:193
iasi_l1_t::ntrack
size_t ntrack
Number of along-track values.
Definition: libiasi.h:169
iasi_l2_t
IASI Level-2 data.
Definition: libiasi.h:198
iasi_l2_t::time
double time[L2_NTRACK][L2_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libiasi.h:204
iasi_l2_t::lon
double lon[L2_NTRACK][L2_NXTRACK]
Longitude [deg].
Definition: libiasi.h:210
iasi_l2_t::p
double p[L2_NLAY]
Pressure [hPa].
Definition: libiasi.h:216
iasi_l2_t::lat
double lat[L2_NTRACK][L2_NXTRACK]
Latitude [deg].
Definition: libiasi.h:213
iasi_l2_t::t
double t[L2_NTRACK][L2_NXTRACK][L2_NLAY]
Temperature [K].
Definition: libiasi.h:219
iasi_l2_t::z
double z[L2_NTRACK][L2_NXTRACK][L2_NLAY]
Geopotential height [km].
Definition: libiasi.h:207
iasi_l2_t::ntrack
size_t ntrack
Number of along-track values.
Definition: libiasi.h:201
iasi_rad_t
IASI converted Level-1 radiation data.
Definition: libiasi.h:288
iasi_rad_t::Longitude
double Longitude[L1_NTRACK][L1_NXTRACK]
Longitude of the sounder pixel.
Definition: libiasi.h:300
iasi_rad_t::Latitude
double Latitude[L1_NTRACK][L1_NXTRACK]
Latitude of the sounder pixel.
Definition: libiasi.h:303
iasi_rad_t::Time
double Time[L1_NTRACK][L1_NXTRACK]
Seconds since 2000-01-01 for each sounder pixel.
Definition: libiasi.h:297
iasi_rad_t::freq
double freq[IASI_L1_NCHAN]
channel wavenumber [cm^-1]
Definition: libiasi.h:294
iasi_rad_t::ntrack
int ntrack
Number of along-track samples.
Definition: libiasi.h:291
iasi_rad_t::Sat_lon
double Sat_lon[L1_NTRACK]
Estimated longitude of the satellite.
Definition: libiasi.h:312
iasi_rad_t::Sat_z
double Sat_z[L1_NTRACK]
Altitude of the satellite.
Definition: libiasi.h:309
iasi_rad_t::Rad
float Rad[L1_NTRACK][L1_NXTRACK][IASI_L1_NCHAN]
Radiance [W/(m^2 sr cm^-1)].
Definition: libiasi.h:306
iasi_rad_t::Sat_lat
double Sat_lat[L1_NTRACK]
Estimated latitude of the satellite.
Definition: libiasi.h:315
iasi_raw_t
IASI raw Level-1 data.
Definition: libiasi.h:256
iasi_raw_t::Wavenumber
float Wavenumber[IASI_L1_NCHAN]
Wavenumbers are computed with the expected values.
Definition: libiasi.h:277
iasi_raw_t::IDefSpectDWn1b
float IDefSpectDWn1b[L1_NTRACK]
Constants for radiation spectrum (must be equal to the expected constant).
Definition: libiasi.h:262
iasi_raw_t::IDefNslast1b
int32_t IDefNslast1b[L1_NTRACK]
Constants for radiation spectrum (must be equal to the expected constant).
Definition: libiasi.h:268
iasi_raw_t::Time
double Time[L1_NTRACK][IASI_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libiasi.h:271
iasi_raw_t::ntrack
long ntrack
Number of along-track samples.
Definition: libiasi.h:259
iasi_raw_t::Radiation
short int Radiation[L1_NTRACK][IASI_NXTRACK][IASI_PM][IASI_L1_NCHAN]
Radiance [W/(m^2 sr m^-1)].
Definition: libiasi.h:280
iasi_raw_t::Sat_z
unsigned int Sat_z[L1_NTRACK]
Satellite altitude [m].
Definition: libiasi.h:283
iasi_raw_t::IDefNsfirst1b
int32_t IDefNsfirst1b[L1_NTRACK]
Constants for radiation spectrum (must be equal to the expected constant).
Definition: libiasi.h:265
iasi_raw_t::Loc
double Loc[L1_NTRACK][IASI_NXTRACK][IASI_PM][2]
Location of the sounder pixel (long,lat).
Definition: libiasi.h:274
orbit_scan_t
Definition: libiasi.c:536
orbit_scan_t::scan
int scan
Definition: libiasi.c:538
orbit_scan_t::orbit
int orbit
Definition: libiasi.c:537
pert_t
Perturbation data.
Definition: libiasi.h:224
pert_t::var
double var[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature variance (4 or 15 micron) [K].
Definition: libiasi.h:251
pert_t::pt
double pt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature perturbation (4 or 15 micron) [K].
Definition: libiasi.h:248
pert_t::time
double time[PERT_NTRACK][PERT_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libiasi.h:233
pert_t::dc
double dc[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (8 micron) [K].
Definition: libiasi.h:242
pert_t::ntrack
int ntrack
Number of along-track values.
Definition: libiasi.h:227
pert_t::nxtrack
int nxtrack
Number of across-track values.
Definition: libiasi.h:230
pert_t::lon
double lon[PERT_NTRACK][PERT_NXTRACK]
Longitude [deg].
Definition: libiasi.h:236
pert_t::bt
double bt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (4 or 15 micron) [K].
Definition: libiasi.h:245
pert_t::lat
double lat[PERT_NTRACK][PERT_NXTRACK]
Latitude [deg].
Definition: libiasi.h:239
wave_t
Wave analysis data.
Definition: libiasi.h:320
wave_t::nx
int nx
Number of across-track values.
Definition: libiasi.h:323
wave_t::ny
int ny
Number of along-track values.
Definition: libiasi.h:326
wave_t::var
double var[WX][WY]
Variance [K].
Definition: libiasi.h:356
wave_t::temp
double temp[WX][WY]
Temperature [K].
Definition: libiasi.h:347
wave_t::x
double x[WX]
Across-track distance [km].
Definition: libiasi.h:341
wave_t::lon
double lon[WX][WY]
Longitude [deg].
Definition: libiasi.h:335
wave_t::y
double y[WY]
Along-track distance [km].
Definition: libiasi.h:344
wave_t::lat
double lat[WX][WY]
Latitude [deg].
Definition: libiasi.h:338
wave_t::z
double z
Altitude [km].
Definition: libiasi.h:332
wave_t::bg
double bg[WX][WY]
Background [K].
Definition: libiasi.h:350
wave_t::time
double time
Time (seconds since 2000-01-01T00:00Z).
Definition: libiasi.h:329
wave_t::pt
double pt[WX][WY]
Perturbation [K].
Definition: libiasi.h:353
Generated by doxygen 1.9.4