IASI Code Collection
Data Structures | Macros | Functions
libiasi.h File Reference

IASI Code Collection library declarations. More...

#include <netcdf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_spline.h>
#include "coda.h"
#include "jurassic.h"

Go to the source code of this file.

Data Structures

struct  iasi_l1_t
 IASI Level-1 data. More...
 
struct  iasi_l2_t
 IASI Level-2 data. More...
 
struct  pert_t
 Perturbation data. More...
 
struct  iasi_raw_t
 IASI raw Level-1 data. More...
 
struct  iasi_rad_t
 IASI converted Level-1 radiation data. More...
 
struct  wave_t
 Wave analysis data. More...
 

Macros

#define L1_NCHAN   33
 Number of IASI radiance channels (don't change). More...
 
#define L1_NTRACK   1800
 Maximum along-track size of IASI radiance granule (don't change). More...
 
#define L1_NXTRACK   60
 Across-track size of IASI radiance granule (don't change). More...
 
#define L2_NLAY   27
 Number of IASI pressure layers (don't change). More...
 
#define L2_NTRACK   1800
 Maximum along-track size of IASI retrieval granule (don't change). More...
 
#define L2_NXTRACK   60
 Across-track size of IASI retrieval granule (don't change). More...
 
#define IASI_L1_NCHAN   8700
 Number of channels of IASI radiance granule. More...
 
#define IASI_NXTRACK   30
 Raw data across-track size of IASI radiance granule. More...
 
#define IASI_PM   4
 Raw data size of measurement matrix (2x2). More...
 
#define IASI_IDefNsfirst1b   2581
 Expected value for the computation of the first wavenumber. More...
 
#define IASI_IDefNslast1b   11041
 Expected value for the computation of the last wavenumber. More...
 
#define IASI_IDefSpectDWn1b   25
 Expected value for the interval of the IASI wavenumbers [m^-1]. More...
 
#define PERT_NTRACK   132000
 Along-track size of perturbation data. More...
 
#define PERT_NXTRACK   360
 Across-track size of perturbation data. More...
 
#define WX   300
 Across-track size of wave analysis data. More...
 
#define WY   33000
 Along-track size of wave analysis data. More...
 
#define CODA(cmd)
 Execute CODA library command and check result. More...
 
#define NC(cmd)
 Execute netCDF library command and check result. More...
 

Functions

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. More...
 
void background_poly (wave_t *wave, int dim_x, int dim_y)
 Get background based on polynomial fits. More...
 
void background_poly_help (double *xx, double *yy, int n, int dim)
 Get background based on polynomial fits. More...
 
void background_smooth (wave_t *wave, int npts_x, int npts_y)
 Smooth background. More...
 
void gauss (wave_t *wave, double fwhm)
 Apply Gaussian filter to perturbations... More...
 
void hamming (wave_t *wave, int nit)
 Apply Hamming filter to perturbations... More...
 
void iasi_read (int format, char *filename, iasi_rad_t *iasi_rad)
 Read IASI Level-1 data. More...
 
void iasi_read_native (char *filename, iasi_rad_t *iasi_rad)
 Read IASI Level-1 data from native file. More...
 
void iasi_read_netcdf (char *filename, iasi_rad_t *iasi_rad)
 Read IASI Level-1 data from netCDF file. More...
 
void median (wave_t *wave, int dx)
 Apply median filter to perturbations... More...
 
void noise (wave_t *wave, double *mu, double *sig)
 Estimate noise. More...
 
void pert2wave (pert_t *pert, wave_t *wave, int track0, int track1, int xtrack0, int xtrack1)
 Convert radiance perturbation data to wave analysis struct. More...
 
void read_pert (char *filename, char *pertname, pert_t *pert)
 Read radiance perturbation data. More...
 
void variance (wave_t *wave, double dh)
 Compute local variance. More...
 
double wgs84 (double lat)
 Calculate Earth radius according to WGS-84 reference ellipsoid. More...
 
void write_l1 (char *filename, iasi_l1_t *l1)
 Write IASI Level-1 data. More...
 
void write_l2 (char *filename, iasi_l2_t *l2)
 Write IASI Level-2 data. More...
 

Detailed Description

IASI Code Collection library declarations.

Definition in file libiasi.h.

Macro Definition Documentation

◆ L1_NCHAN

#define L1_NCHAN   33

Number of IASI radiance channels (don't change).

Definition at line 96 of file libiasi.h.

◆ L1_NTRACK

#define L1_NTRACK   1800

Maximum along-track size of IASI radiance granule (don't change).

Definition at line 99 of file libiasi.h.

◆ L1_NXTRACK

#define L1_NXTRACK   60

Across-track size of IASI radiance granule (don't change).

Definition at line 102 of file libiasi.h.

◆ L2_NLAY

#define L2_NLAY   27

Number of IASI pressure layers (don't change).

Definition at line 105 of file libiasi.h.

◆ L2_NTRACK

#define L2_NTRACK   1800

Maximum along-track size of IASI retrieval granule (don't change).

Definition at line 108 of file libiasi.h.

◆ L2_NXTRACK

#define L2_NXTRACK   60

Across-track size of IASI retrieval granule (don't change).

Definition at line 111 of file libiasi.h.

◆ IASI_L1_NCHAN

#define IASI_L1_NCHAN   8700

Number of channels of IASI radiance granule.

Definition at line 114 of file libiasi.h.

◆ IASI_NXTRACK

#define IASI_NXTRACK   30

Raw data across-track size of IASI radiance granule.

Definition at line 117 of file libiasi.h.

◆ IASI_PM

#define IASI_PM   4

Raw data size of measurement matrix (2x2).

Definition at line 120 of file libiasi.h.

◆ IASI_IDefNsfirst1b

#define IASI_IDefNsfirst1b   2581

Expected value for the computation of the first wavenumber.

Definition at line 123 of file libiasi.h.

◆ IASI_IDefNslast1b

#define IASI_IDefNslast1b   11041

Expected value for the computation of the last wavenumber.


Definition at line 126 of file libiasi.h.

◆ IASI_IDefSpectDWn1b

#define IASI_IDefSpectDWn1b   25

Expected value for the interval of the IASI wavenumbers [m^-1].

Definition at line 129 of file libiasi.h.

◆ PERT_NTRACK

#define PERT_NTRACK   132000

Along-track size of perturbation data.

Definition at line 132 of file libiasi.h.

◆ PERT_NXTRACK

#define PERT_NXTRACK   360

Across-track size of perturbation data.

Definition at line 135 of file libiasi.h.

◆ WX

#define WX   300

Across-track size of wave analysis data.

Definition at line 138 of file libiasi.h.

◆ WY

#define WY   33000

Along-track size of wave analysis data.

Definition at line 141 of file libiasi.h.

◆ CODA

#define CODA (   cmd)
Value:
{ \
int coda_result=(cmd); \
if(coda_result!=0) \
ERRMSG("%s", coda_errno_to_string(coda_errno)); \
}

Execute CODA library command and check result.

Definition at line 148 of file libiasi.h.

◆ NC

#define NC (   cmd)
Value:
{ \
int nc_result=(cmd); \
if(nc_result!=NC_NOERR) \
ERRMSG("%s", nc_strerror(nc_result)); \
}

Execute netCDF library command and check result.

Definition at line 155 of file libiasi.h.

Function Documentation

◆ 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.

Add variable to netCDF file.

Definition at line 30 of file libiasi.c.

38 {
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}
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: libiasi.h:155

◆ background_poly()

void background_poly ( wave_t *  wave,
int  dim_x,
int  dim_y 
)

Get background based on polynomial fits.

Definition at line 114 of file libiasi.c.

117 {
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}
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
WX
#define WX
Across-track size of wave analysis data.
Definition: libiasi.h:138
WY
#define WY
Along-track size of wave analysis data.
Definition: libiasi.h:141
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::temp
double temp[WX][WY]
Temperature [K].
Definition: libiasi.h:347
wave_t::bg
double bg[WX][WY]
Background [K].
Definition: libiasi.h:350
wave_t::pt
double pt[WX][WY]
Perturbation [K].
Definition: libiasi.h:353
Here is the call graph for this function:

◆ background_poly_help()

void background_poly_help ( double *  xx,
double *  yy,
int  n,
int  dim 
)

Get background based on polynomial fits.

Definition at line 57 of file libiasi.c.

61 {
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}

◆ background_smooth()

void background_smooth ( wave_t *  wave,
int  npts_x,
int  npts_y 
)

Smooth background.

Definition at line 177 of file libiasi.c.

180 {
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}
wave_t::x
double x[WX]
Across-track distance [km].
Definition: libiasi.h:341
wave_t::y
double y[WY]
Along-track distance [km].
Definition: libiasi.h:344

◆ gauss()

void gauss ( wave_t *  wave,
double  fwhm 
)

Apply Gaussian filter to perturbations...

Definition at line 228 of file libiasi.c.

230 {
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}

◆ hamming()

void hamming ( wave_t *  wave,
int  nit 
)

Apply Hamming filter to perturbations...

Definition at line 268 of file libiasi.c.

270 {
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}

◆ iasi_read()

void iasi_read ( int  format,
char *  filename,
iasi_rad_t *  iasi_rad 
)

Read IASI Level-1 data.

Definition at line 297 of file libiasi.c.

300 {
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}
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
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
Here is the call graph for this function:

◆ iasi_read_native()

void iasi_read_native ( char *  filename,
iasi_rad_t *  iasi_rad 
)

Read IASI Level-1 data from native file.

Definition at line 317 of file libiasi.c.

319 {
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}
wgs84
double wgs84(double lat)
Calculate Earth radius according to WGS-84 reference ellipsoid.
Definition: libiasi.c:1143
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
L1_NXTRACK
#define L1_NXTRACK
Across-track size of IASI radiance granule (don't change).
Definition: libiasi.h:102
IASI_L1_NCHAN
#define IASI_L1_NCHAN
Number of channels of IASI radiance granule.
Definition: libiasi.h:114
IASI_IDefSpectDWn1b
#define IASI_IDefSpectDWn1b
Expected value for the interval of the IASI wavenumbers [m^-1].
Definition: libiasi.h:129
IASI_IDefNsfirst1b
#define IASI_IDefNsfirst1b
Expected value for the computation of the first wavenumber.
Definition: libiasi.h:123
CODA
#define CODA(cmd)
Execute CODA library command and check result.
Definition: libiasi.h:148
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
Here is the call graph for this function:

◆ iasi_read_netcdf()

void iasi_read_netcdf ( char *  filename,
iasi_rad_t *  iasi_rad 
)

Read IASI Level-1 data from netCDF file.

Definition at line 557 of file libiasi.c.

559 {
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 var_qualflag = -1, *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 int *qualflag_all = NULL;
578
579 /* Scanline keys... */
580 orbit_scan_t *keys_all = NULL; /* length npoint */
581 orbit_scan_t *keys_uniq = NULL; /* length <= npoint */
582 size_t nuniq = 0;
583
584 /* Always leave struct safe... */
585 iasi_rad->ntrack = 0;
586
587 /* Standard IASI wavenumber grid... */
588 for (int ichan = 0; ichan < IASI_L1_NCHAN; ichan++)
589 iasi_rad->freq[ichan] = 644.75 + (double) (ichan + 1) * 0.25;
590
591 /* Open file... */
592 if (nc_open(filename, NC_NOWRITE, &ncid) != NC_NOERR) {
593 WARN("Cannot open file %s", filename);
594 return;
595 }
596
597 /* Get dimensions... */
598 NC(nc_inq_dimid(ncid, "point", &dim_point_id));
599 NC(nc_inq_dimlen(ncid, dim_point_id, &npoint));
600 LOG(2, "number of points: %lu", npoint);
601 NC(nc_inq_dimid(ncid, "channel", &dim_chan_id));
602 NC(nc_inq_dimlen(ncid, dim_chan_id, &nchan));
603 LOG(2, "number of channels: %lu", nchan);
604 if (npoint == 0 || nchan == 0) {
605 WARN("Empty file, npoint and/or nchan is zero!");
606 NC(nc_close(ncid));
607 return;
608 }
609
610 /* Get variables... */
611 NC(nc_inq_varid(ncid, "lat", &var_lat));
612 NC(nc_inq_varid(ncid, "lon", &var_lon));
613 NC(nc_inq_varid(ncid, "date", &var_date));
614 NC(nc_inq_varid(ncid, "orbit", &var_orbit));
615 NC(nc_inq_varid(ncid, "scan", &var_scan));
616 NC(nc_inq_varid(ncid, "pixel", &var_pixel));
617 NC(nc_inq_varid(ncid, "fov", &var_fov));
618 NC(nc_inq_varid(ncid, "channel_name", &var_channame));
619 NC(nc_inq_varid(ncid, "R", &var_R));
620 NC(nc_inq_varid(ncid, "qualflag", &var_qualflag));
621
622 /* Channel mapping... */
623 ALLOC(channel_name, int,
624 nchan);
625 NC(nc_get_var_int(ncid, var_channame, channel_name));
626 LOG(2, "channels: name= %d ... %d | freq= %g ... %g cm^-1", channel_name[0],
627 channel_name[nchan - 1], iasi_rad->freq[channel_name[0] - 1],
628 iasi_rad->freq[channel_name[nchan - 1] - 1]);
629
630 /* Allocate full arrays... */
631 ALLOC(orbit_all, int,
632 npoint);
633 ALLOC(scan_all, int,
634 npoint);
635 ALLOC(pixel_all, int,
636 npoint);
637 ALLOC(fov_all, int,
638 npoint);
639 ALLOC(lat_all, double,
640 npoint);
641 ALLOC(lon_all, double,
642 npoint);
643 ALLOC(date_all, double,
644 npoint);
645 ALLOC(R_all, float,
646 npoint * nchan);
647 ALLOC(qualflag_all, int,
648 npoint);
649
650 /* Read everything... */
651 NC(nc_get_var_int(ncid, var_orbit, orbit_all));
652 NC(nc_get_var_int(ncid, var_scan, scan_all));
653 NC(nc_get_var_int(ncid, var_pixel, pixel_all));
654 NC(nc_get_var_int(ncid, var_fov, fov_all));
655 NC(nc_get_var_double(ncid, var_lat, lat_all));
656 NC(nc_get_var_double(ncid, var_lon, lon_all));
657 NC(nc_get_var_double(ncid, var_date, date_all));
658 NC(nc_get_var_float(ncid, var_R, R_all));
659 NC(nc_get_var_int(ncid, var_qualflag, qualflag_all));
660
661 /* Build list of (orbit,scan) keys for each point... */
662 ALLOC(keys_all, orbit_scan_t, npoint);
663 for (size_t i = 0; i < npoint; i++) {
664 keys_all[i].orbit = orbit_all[i];
665 keys_all[i].scan = scan_all[i];
666 }
667
668 /* Sort keys_all, then unique into keys_uniq.. */
669 qsort(keys_all, npoint, sizeof(orbit_scan_t), cmp_orbit_scan);
670 ALLOC(keys_uniq, orbit_scan_t, npoint);
671 nuniq = 0;
672 for (size_t i = 0; i < npoint; i++)
673 if (i == 0 || cmp_orbit_scan(&keys_all[i], &keys_all[i - 1]) != 0)
674 keys_uniq[nuniq++] = keys_all[i];
675
676 /* Early check: file must fit completely... */
677 if (nuniq > (size_t) (L1_NTRACK / 2))
678 ERRMSG("Too many scanlines in file: %zu (max %d). Increase L1_NTRACK!",
679 nuniq, L1_NTRACK / 2);
680
681 /* ntrack = 2 * number of unique scanlines */
682 iasi_rad->ntrack = (int) (2 * nuniq);
683
684 /* Initialize radiances with NaN (for missing pixels)... */
685 for (int tr = 0; tr < iasi_rad->ntrack; tr++)
686 for (int ix = 0; ix < L1_NXTRACK; ix++)
687 for (int g = 0; g < IASI_L1_NCHAN; g++)
688 iasi_rad->Rad[tr][ix][g] = GSL_NAN;
689
690 /* Initialize footprint time and location with NaN... */
691 for (int track = 0; track < iasi_rad->ntrack; track++)
692 for (int xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
693 iasi_rad->Time[track][xtrack] = GSL_NAN;
694 iasi_rad->Longitude[track][xtrack] = GSL_NAN;
695 iasi_rad->Latitude[track][xtrack] = GSL_NAN;
696 }
697
698 /* No satellite position info, set to NaN... */
699 for (int tr = 0; tr < iasi_rad->ntrack; tr++) {
700 iasi_rad->Sat_z[tr] = GSL_NAN;
701 iasi_rad->Sat_lon[tr] = GSL_NAN;
702 iasi_rad->Sat_lat[tr] = GSL_NAN;
703 }
704
705 /* Fill iasi_rad... */
706 for (size_t i = 0; i < npoint; i++) {
707
708 /* Find scanline index from (orbit,scan)... */
709 orbit_scan_t key;
710 key.orbit = orbit_all[i];
711 key.scan = scan_all[i];
712 orbit_scan_t *found =
713 (orbit_scan_t *) bsearch(&key, keys_uniq, nuniq, sizeof(orbit_scan_t),
714 cmp_orbit_scan);
715 if (!found)
716 continue;
717 size_t scan_idx = (size_t) (found - keys_uniq);
718
719 int f = fov_all[i];
720 if (f < 1 || f > 120)
721 continue;
722
723 /* Get position... */
724 int pos = (f - 1) / 4; /* 0..29 */
725 int sub = (f - 1) % 4; /* 0..3 */
726 if (pos < 0 || pos >= IASI_NXTRACK)
727 continue;
728
729 int track_add, x_add;
730 switch (sub) {
731 case 3:
732 track_add = 0;
733 x_add = 0;
734 break; /* tr1_lpm */
735 case 0:
736 track_add = 0;
737 x_add = 1;
738 break; /* tr1_rpm */
739 case 2:
740 track_add = 1;
741 x_add = 0;
742 break; /* tr2_lpm */
743 case 1:
744 track_add = 1;
745 x_add = 1;
746 break; /* tr2_rpm */
747 default:
748 continue;
749 }
750
751 int tr = 2 * (int) scan_idx + track_add;
752 int ix = 2 * pos + x_add;
753
754 if (tr < 0 || tr >= iasi_rad->ntrack)
755 continue;
756 if (ix < 0 || ix >= L1_NXTRACK)
757 continue;
758
759 /* Save data... */
760 iasi_rad->Time[tr][ix] = date_all[i];
761 iasi_rad->Longitude[tr][ix] = lon_all[i];
762 iasi_rad->Latitude[tr][ix] = lat_all[i];
763 float *row = &R_all[i * nchan];
764 if (qualflag_all[i] != 0)
765 continue; // keep NaNs for this point
766 for (size_t k = 0; k < nchan; k++) {
767 const int g = channel_name[k] - 1;
768 if (g >= 0 && g < IASI_L1_NCHAN)
769 iasi_rad->Rad[tr][ix][g] = 100.0f * row[k]; /* m^-1 -> cm^-1 */
770 }
771 }
772
773 /* Close file... */
774 NC(nc_close(ncid));
775
776 /* Check time... */
777 double tmin = 1e100, tmax = -1e100;
778 for (int track = 0; track < iasi_rad->ntrack; track++)
779 for (int xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
780 if (gsl_finite(iasi_rad->Time[track][xtrack])) {
781 tmin = GSL_MIN(tmin, iasi_rad->Time[track][xtrack]);
782 tmax = GSL_MAX(tmax, iasi_rad->Time[track][xtrack]);
783 }
784 if (track > 0 && gsl_finite(iasi_rad->Time[track][xtrack])
785 && gsl_finite(iasi_rad->Time[track - 1][xtrack])
786 && iasi_rad->Time[track][xtrack] <
787 iasi_rad->Time[track - 1][xtrack])
788 ERRMSG("Time is not ascending!");
789 }
790 LOG(2, "time= %.2f ... %.2f s", tmin, tmax);
791
792 /* Free... */
793 free(channel_name);
794 free(orbit_all);
795 free(scan_all);
796 free(pixel_all);
797 free(fov_all);
798 free(lat_all);
799 free(lon_all);
800 free(date_all);
801 free(R_all);
802 free(qualflag_all);
803 free(keys_all);
804 free(keys_uniq);
805}
L1_NTRACK
#define L1_NTRACK
Maximum along-track size of IASI radiance granule (don't change).
Definition: libiasi.h:99
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

◆ median()

void median ( wave_t *  wave,
int  dx 
)

Apply median filter to perturbations...

Definition at line 809 of file libiasi.c.

811 {
812
813 static double data[WX * WY], help[WX][WY];
814
815 /* Check parameters... */
816 if (dx <= 0)
817 return;
818
819 /* Loop over data points... */
820 for (int ix = 0; ix < wave->nx; ix++)
821 for (int iy = 0; iy < wave->ny; iy++) {
822
823 /* Init... */
824 size_t n = 0;
825
826 /* Get data... */
827 for (int ix2 = GSL_MAX(ix - dx, 0);
828 ix2 < GSL_MIN(ix + dx, wave->nx - 1); ix2++)
829 for (int iy2 = GSL_MAX(iy - dx, 0);
830 iy2 < GSL_MIN(iy + dx, wave->ny - 1); iy2++) {
831 data[n] = wave->pt[ix2][iy2];
832 n++;
833 }
834
835 /* Normalize... */
836 gsl_sort(data, 1, n);
837 help[ix][iy] = gsl_stats_median_from_sorted_data(data, 1, n);
838 }
839
840 /* Loop over data points... */
841 for (int ix = 0; ix < wave->nx; ix++)
842 for (int iy = 0; iy < wave->ny; iy++)
843 wave->pt[ix][iy] = help[ix][iy];
844}

◆ noise()

void noise ( wave_t *  wave,
double *  mu,
double *  sig 
)

Estimate noise.

Definition at line 848 of file libiasi.c.

851 {
852
853 int ix, ix2, iy, iy2, n = 0, okay;
854
855 /* Init... */
856 *mu = 0;
857 *sig = 0;
858
859 /* Estimate noise (Immerkaer, 1996)... */
860 for (ix = 1; ix < wave->nx - 1; ix++)
861 for (iy = 1; iy < wave->ny - 1; iy++) {
862
863 /* Check data... */
864 okay = 1;
865 for (ix2 = ix - 1; ix2 <= ix + 1; ix2++)
866 for (iy2 = iy - 1; iy2 <= iy + 1; iy2++)
867 if (!gsl_finite(wave->temp[ix2][iy2]))
868 okay = 0;
869 if (!okay)
870 continue;
871
872 /* Get mean noise... */
873 n++;
874 *mu += wave->temp[ix][iy];
875 *sig += gsl_pow_2(+4. / 6. * wave->temp[ix][iy]
876 - 2. / 6. * (wave->temp[ix - 1][iy]
877 + wave->temp[ix + 1][iy]
878 + wave->temp[ix][iy - 1]
879 + wave->temp[ix][iy + 1])
880 + 1. / 6. * (wave->temp[ix - 1][iy - 1]
881 + wave->temp[ix + 1][iy - 1]
882 + wave->temp[ix - 1][iy + 1]
883 + wave->temp[ix + 1][iy + 1]));
884 }
885
886 /* Normalize... */
887 *mu /= (double) n;
888 *sig = sqrt(*sig / (double) n);
889}

◆ 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 at line 893 of file libiasi.c.

899 {
900
901 double x0[3], x1[3];
902
903 int itrack, ixtrack;
904
905 /* Check ranges... */
906 track0 = GSL_MIN(GSL_MAX(track0, 0), pert->ntrack - 1);
907 track1 = GSL_MIN(GSL_MAX(track1, 0), pert->ntrack - 1);
908 xtrack0 = GSL_MIN(GSL_MAX(xtrack0, 0), pert->nxtrack - 1);
909 xtrack1 = GSL_MIN(GSL_MAX(xtrack1, 0), pert->nxtrack - 1);
910
911 /* Set size... */
912 wave->nx = xtrack1 - xtrack0 + 1;
913 if (wave->nx > WX)
914 ERRMSG("Too many across-track values!");
915 wave->ny = track1 - track0 + 1;
916 if (wave->ny > WY)
917 ERRMSG("Too many along-track values!");
918
919 /* ------------------------------------------------------------------------- */
920 /* Compute wave->x[] and wave->y[] as mean step sizes over the whole window, */
921 /* ignoring any lon/lat pairs that contain NaNs (or non-finite values). */
922 /* ------------------------------------------------------------------------- */
923
924 wave->x[0] = 0.0;
925 for (ixtrack = xtrack0 + 1; ixtrack <= xtrack1; ixtrack++) {
926
927 double sum = 0.0;
928 int n = 0;
929
930 /* mean distance between columns (ixtrack-1) and (ixtrack), averaged over all rows */
931 for (itrack = track0; itrack <= track1; itrack++) {
932
933 double lon0 = pert->lon[itrack][ixtrack - 1];
934 double lat0 = pert->lat[itrack][ixtrack - 1];
935 double lon1 = pert->lon[itrack][ixtrack];
936 double lat1 = pert->lat[itrack][ixtrack];
937
938 if (!isfinite(lon0) || !isfinite(lat0) || !isfinite(lon1)
939 || !isfinite(lat1))
940 continue;
941
942 geo2cart(0, lon0, lat0, x0);
943 geo2cart(0, lon1, lat1, x1);
944 sum += DIST(x0, x1);
945 n++;
946 }
947
948 /* If no valid pairs exist for this step, add 0 (keeps x finite/monotonic). */
949 wave->x[ixtrack - xtrack0] =
950 wave->x[ixtrack - xtrack0 - 1] + ((n > 0) ? (sum / (double) n) : 0.0);
951 }
952
953 wave->y[0] = 0.0;
954 for (itrack = track0 + 1; itrack <= track1; itrack++) {
955
956 double sum = 0.0;
957 int n = 0;
958
959 /* mean distance between rows (itrack-1) and (itrack), averaged over all columns */
960 for (ixtrack = xtrack0; ixtrack <= xtrack1; ixtrack++) {
961
962 double lon0 = pert->lon[itrack - 1][ixtrack];
963 double lat0 = pert->lat[itrack - 1][ixtrack];
964 double lon1 = pert->lon[itrack][ixtrack];
965 double lat1 = pert->lat[itrack][ixtrack];
966
967 if (!isfinite(lon0) || !isfinite(lat0) || !isfinite(lon1)
968 || !isfinite(lat1))
969 continue;
970
971 geo2cart(0, lon0, lat0, x0);
972 geo2cart(0, lon1, lat1, x1);
973 sum += DIST(x0, x1);
974 n++;
975 }
976
977 wave->y[itrack - track0] =
978 wave->y[itrack - track0 - 1] + ((n > 0) ? (sum / (double) n) : 0.0);
979 }
980
981 /* Loop over footprints... */
982 for (itrack = track0; itrack <= track1; itrack++)
983 for (ixtrack = xtrack0; ixtrack <= xtrack1; ixtrack++) {
984
985 /* Save geolocation... */
986 wave->time = pert->time[(track0 + track1) / 2][(xtrack0 + xtrack1) / 2];
987 wave->z = 0;
988 wave->lon[ixtrack - xtrack0][itrack - track0] =
989 pert->lon[itrack][ixtrack];
990 wave->lat[ixtrack - xtrack0][itrack - track0] =
991 pert->lat[itrack][ixtrack];
992
993 /* Save temperature data... */
994 wave->temp[ixtrack - xtrack0][itrack - track0]
995 = pert->bt[itrack][ixtrack];
996 wave->bg[ixtrack - xtrack0][itrack - track0]
997 = pert->bt[itrack][ixtrack] - pert->pt[itrack][ixtrack];
998 wave->pt[ixtrack - xtrack0][itrack - track0]
999 = pert->pt[itrack][ixtrack];
1000 wave->var[ixtrack - xtrack0][itrack - track0]
1001 = pert->var[itrack][ixtrack];
1002 }
1003}
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::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::var
double var[WX][WY]
Variance [K].
Definition: libiasi.h:356
wave_t::lon
double lon[WX][WY]
Longitude [deg].
Definition: libiasi.h:335
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::time
double time
Time (seconds since 2000-01-01T00:00Z).
Definition: libiasi.h:329

◆ read_pert()

void read_pert ( char *  filename,
char *  pertname,
pert_t *  pert 
)

Read radiance perturbation data.

Definition at line 1007 of file libiasi.c.

1010 {
1011
1012 static char varname[LEN];
1013
1014 static int dimid[2], ncid, varid;
1015
1016 static size_t ntrack, nxtrack, start[2] = { 0, 0 }, count[2] = { 1, 1 };
1017
1018 /* Write info... */
1019 LOG(1, "Read perturbation data: %s", filename);
1020
1021 /* Open netCDF file... */
1022 NC(nc_open(filename, NC_NOWRITE, &ncid));
1023
1024 /* Get dimensions... */
1025 NC(nc_inq_dimid(ncid, "NTRACK", &dimid[0]));
1026 NC(nc_inq_dimid(ncid, "NXTRACK", &dimid[1]));
1027 NC(nc_inq_dimlen(ncid, dimid[0], &ntrack));
1028 NC(nc_inq_dimlen(ncid, dimid[1], &nxtrack));
1029 if (nxtrack > PERT_NXTRACK)
1030 ERRMSG("Too many tracks!");
1031 if (ntrack > PERT_NTRACK)
1032 ERRMSG("Too many scans!");
1033 pert->ntrack = (int) ntrack;
1034 pert->nxtrack = (int) nxtrack;
1035 count[1] = nxtrack;
1036
1037 /* Read data... */
1038 NC(nc_inq_varid(ncid, "time", &varid));
1039 for (size_t itrack = 0; itrack < ntrack; itrack++) {
1040 start[0] = itrack;
1041 NC(nc_get_vara_double(ncid, varid, start, count, pert->time[itrack]));
1042 }
1043
1044 NC(nc_inq_varid(ncid, "lon", &varid));
1045 for (size_t itrack = 0; itrack < ntrack; itrack++) {
1046 start[0] = itrack;
1047 NC(nc_get_vara_double(ncid, varid, start, count, pert->lon[itrack]));
1048 }
1049
1050 NC(nc_inq_varid(ncid, "lat", &varid));
1051 for (size_t itrack = 0; itrack < ntrack; itrack++) {
1052 start[0] = itrack;
1053 NC(nc_get_vara_double(ncid, varid, start, count, pert->lat[itrack]));
1054 }
1055
1056 NC(nc_inq_varid(ncid, "bt_8mu", &varid));
1057 for (size_t itrack = 0; itrack < ntrack; itrack++) {
1058 start[0] = itrack;
1059 NC(nc_get_vara_double(ncid, varid, start, count, pert->dc[itrack]));
1060 }
1061
1062 sprintf(varname, "bt_%s", pertname);
1063 NC(nc_inq_varid(ncid, varname, &varid));
1064 for (size_t itrack = 0; itrack < ntrack; itrack++) {
1065 start[0] = itrack;
1066 NC(nc_get_vara_double(ncid, varid, start, count, pert->bt[itrack]));
1067 }
1068
1069 sprintf(varname, "bt_%s_pt", pertname);
1070 NC(nc_inq_varid(ncid, varname, &varid));
1071 for (size_t itrack = 0; itrack < ntrack; itrack++) {
1072 start[0] = itrack;
1073 NC(nc_get_vara_double(ncid, varid, start, count, pert->pt[itrack]));
1074 }
1075
1076 sprintf(varname, "bt_%s_var", pertname);
1077 NC(nc_inq_varid(ncid, varname, &varid));
1078 for (size_t itrack = 0; itrack < ntrack; itrack++) {
1079 start[0] = itrack;
1080 NC(nc_get_vara_double(ncid, varid, start, count, pert->var[itrack]));
1081 }
1082
1083 /* Close file... */
1084 NC(nc_close(ncid));
1085}
PERT_NXTRACK
#define PERT_NXTRACK
Across-track size of perturbation data.
Definition: libiasi.h:135
PERT_NTRACK
#define PERT_NTRACK
Along-track size of perturbation data.
Definition: libiasi.h:132
pert_t::dc
double dc[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (8 micron) [K].
Definition: libiasi.h:242

◆ variance()

void variance ( wave_t *  wave,
double  dh 
)

Compute local variance.

Definition at line 1089 of file libiasi.c.

1091 {
1092
1093 double dh2, mu, help;
1094
1095 int dx, dy, ix, ix2, iy, iy2, n;
1096
1097 /* Check parameters... */
1098 if (dh <= 0)
1099 return;
1100
1101 /* Compute squared radius... */
1102 dh2 = gsl_pow_2(dh);
1103
1104 /* Get sampling distances... */
1105 dx =
1106 (int) (dh / fabs(wave->x[wave->nx - 1] - wave->x[0]) *
1107 (wave->nx - 1.0) + 1);
1108 dy =
1109 (int) (dh / fabs(wave->y[wave->ny - 1] - wave->y[0]) *
1110 (wave->ny - 1.0) + 1);
1111
1112 /* Loop over data points... */
1113 for (ix = 0; ix < wave->nx; ix++)
1114 for (iy = 0; iy < wave->ny; iy++) {
1115
1116 /* Init... */
1117 mu = help = 0;
1118 n = 0;
1119
1120 /* Get data... */
1121 for (ix2 = GSL_MAX(ix - dx, 0); ix2 <= GSL_MIN(ix + dx, wave->nx - 1);
1122 ix2++)
1123 for (iy2 = GSL_MAX(iy - dy, 0);
1124 iy2 <= GSL_MIN(iy + dy, wave->ny - 1); iy2++)
1125 if ((gsl_pow_2(wave->x[ix] - wave->x[ix2])
1126 + gsl_pow_2(wave->y[iy] - wave->y[iy2])) <= dh2)
1127 if (gsl_finite(wave->pt[ix2][iy2])) {
1128 mu += wave->pt[ix2][iy2];
1129 help += gsl_pow_2(wave->pt[ix2][iy2]);
1130 n++;
1131 }
1132
1133 /* Compute local variance... */
1134 if (n > 1)
1135 wave->var[ix][iy] = help / n - gsl_pow_2(mu / n);
1136 else
1137 wave->var[ix][iy] = GSL_NAN;
1138 }
1139}

◆ wgs84()

double wgs84 ( double  lat)

Calculate Earth radius according to WGS-84 reference ellipsoid.

Definition at line 1143 of file libiasi.c.

1144 {
1145
1146 const double a = 6378.1370, b = 6356.7523;
1147
1148 double cphi, sphi;
1149
1150 cphi = cos(lat * M_PI / 180.);
1151 sphi = sin(lat * M_PI / 180.);
1152
1153 return sqrt((gsl_pow_2(a * a * cphi) + gsl_pow_2(b * b * sphi))
1154 / (gsl_pow_2(a * cphi) + gsl_pow_2(b * sphi)));
1155}

◆ write_l1()

void write_l1 ( char *  filename,
iasi_l1_t *  l1 
)

Write IASI Level-1 data.

Definition at line 1159 of file libiasi.c.

1161 {
1162
1163 int dimid[10], ncid, time_id, lon_id, lat_id,
1164 sat_z_id, sat_lon_id, sat_lat_id, nu_id, rad_id;
1165
1166 /* Open or create netCDF file... */
1167 LOG(1, "Write IASI Level-1 file: %s", filename);
1168 if (nc_open(filename, NC_WRITE, &ncid) != NC_NOERR) {
1169 NC(nc_create(filename, NC_CLOBBER, &ncid));
1170 } else {
1171 NC(nc_redef(ncid));
1172 }
1173
1174 /* Set dimensions... */
1175 if (nc_inq_dimid(ncid, "L1_NTRACK", &dimid[0]) != NC_NOERR)
1176 NC(nc_def_dim(ncid, "L1_NTRACK", l1->ntrack, &dimid[0]));
1177 if (nc_inq_dimid(ncid, "L1_NXTRACK", &dimid[1]) != NC_NOERR)
1178 NC(nc_def_dim(ncid, "L1_NXTRACK", L1_NXTRACK, &dimid[1]));
1179 if (nc_inq_dimid(ncid, "L1_NCHAN", &dimid[2]) != NC_NOERR)
1180 NC(nc_def_dim(ncid, "L1_NCHAN", L1_NCHAN, &dimid[2]));
1181
1182 /* Add variables... */
1183 add_var(ncid, "l1_time", "s", "time (seconds since 2000-01-01T00:00Z)",
1184 NC_DOUBLE, dimid, &time_id, 2);
1185 add_var(ncid, "l1_lon", "deg", "longitude", NC_DOUBLE, dimid, &lon_id, 2);
1186 add_var(ncid, "l1_lat", "deg", "latitude", NC_DOUBLE, dimid, &lat_id, 2);
1187 add_var(ncid, "l1_sat_z", "km", "satellite altitude",
1188 NC_DOUBLE, dimid, &sat_z_id, 1);
1189 add_var(ncid, "l1_sat_lon", "deg", "(estimated) satellite longitude",
1190 NC_DOUBLE, dimid, &sat_lon_id, 1);
1191 add_var(ncid, "l1_sat_lat", "deg", "(estimated) satellite latitude",
1192 NC_DOUBLE, dimid, &sat_lat_id, 1);
1193 add_var(ncid, "l1_nu", "cm^-1", "channel wavenumber",
1194 NC_DOUBLE, &dimid[2], &nu_id, 1);
1195 add_var(ncid, "l1_rad", "W/(m^2 sr cm^-1)", "channel radiance",
1196 NC_FLOAT, dimid, &rad_id, 3);
1197
1198 /* Leave define mode... */
1199 NC(nc_enddef(ncid));
1200
1201 /* Write data... */
1202 NC(nc_put_var_double(ncid, time_id, l1->time[0]));
1203 NC(nc_put_var_double(ncid, lon_id, l1->lon[0]));
1204 NC(nc_put_var_double(ncid, lat_id, l1->lat[0]));
1205 NC(nc_put_var_double(ncid, sat_z_id, l1->sat_z));
1206 NC(nc_put_var_double(ncid, sat_lon_id, l1->sat_lon));
1207 NC(nc_put_var_double(ncid, sat_lat_id, l1->sat_lat));
1208 NC(nc_put_var_double(ncid, nu_id, l1->nu));
1209 NC(nc_put_var_float(ncid, rad_id, &l1->rad[0][0][0]));
1210
1211 /* Close file... */
1212 NC(nc_close(ncid));
1213}
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
L1_NCHAN
#define L1_NCHAN
Number of IASI radiance channels (don't change).
Definition: libiasi.h:96
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
Here is the call graph for this function:

◆ write_l2()

void write_l2 ( char *  filename,
iasi_l2_t *  l2 
)

Write IASI Level-2 data.

Definition at line 1217 of file libiasi.c.

1219 {
1220
1221 int dimid[10], ncid, time_id, z_id, lon_id, lat_id, p_id, t_id;
1222
1223 /* Create netCDF file... */
1224 LOG(1, "Write IASI Level-2 file: %s", filename);
1225 if (nc_open(filename, NC_WRITE, &ncid) != NC_NOERR) {
1226 NC(nc_create(filename, NC_CLOBBER, &ncid));
1227 } else {
1228 NC(nc_redef(ncid));
1229 }
1230
1231 /* Set dimensions... */
1232 if (nc_inq_dimid(ncid, "L2_NTRACK", &dimid[0]) != NC_NOERR)
1233 NC(nc_def_dim(ncid, "L2_NTRACK", l2->ntrack, &dimid[0]));
1234 if (nc_inq_dimid(ncid, "L2_NXTRACK", &dimid[1]) != NC_NOERR)
1235 NC(nc_def_dim(ncid, "L2_NXTRACK", L2_NXTRACK, &dimid[1]));
1236 if (nc_inq_dimid(ncid, "L2_NLAY", &dimid[2]) != NC_NOERR)
1237 NC(nc_def_dim(ncid, "L2_NLAY", L2_NLAY, &dimid[2]));
1238
1239 /* Add variables... */
1240 add_var(ncid, "l2_time", "s", "time (seconds since 2000-01-01T00:00Z)",
1241 NC_DOUBLE, dimid, &time_id, 2);
1242 add_var(ncid, "l2_z", "km", "altitude", NC_DOUBLE, dimid, &z_id, 3);
1243 add_var(ncid, "l2_lon", "deg", "longitude", NC_DOUBLE, dimid, &lon_id, 2);
1244 add_var(ncid, "l2_lat", "deg", "latitude", NC_DOUBLE, dimid, &lat_id, 2);
1245 add_var(ncid, "l2_press", "hPa", "pressure",
1246 NC_DOUBLE, &dimid[2], &p_id, 1);
1247 add_var(ncid, "l2_temp", "K", "temperature", NC_DOUBLE, dimid, &t_id, 3);
1248
1249 /* Leave define mode... */
1250 NC(nc_enddef(ncid));
1251
1252 /* Write data... */
1253 NC(nc_put_var_double(ncid, time_id, l2->time[0]));
1254 NC(nc_put_var_double(ncid, z_id, l2->z[0][0]));
1255 NC(nc_put_var_double(ncid, lon_id, l2->lon[0]));
1256 NC(nc_put_var_double(ncid, lat_id, l2->lat[0]));
1257 NC(nc_put_var_double(ncid, p_id, l2->p));
1258 NC(nc_put_var_double(ncid, t_id, l2->t[0][0]));
1259
1260 /* Close file... */
1261 NC(nc_close(ncid));
1262}
L2_NXTRACK
#define L2_NXTRACK
Across-track size of IASI retrieval granule (don't change).
Definition: libiasi.h:111
L2_NLAY
#define L2_NLAY
Number of IASI pressure layers (don't change).
Definition: libiasi.h:105
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
Here is the call graph for this function:
Generated by doxygen 1.9.4