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

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

◆ noise()

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

Estimate noise.

Definition at line 840 of file libiasi.c.

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

◆ 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 885 of file libiasi.c.

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

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

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

◆ wgs84()

double wgs84 ( double  lat)

Calculate Earth radius according to WGS-84 reference ellipsoid.

Definition at line 1135 of file libiasi.c.

1136 {
1137
1138 const double a = 6378.1370, b = 6356.7523;
1139
1140 double cphi, sphi;
1141
1142 cphi = cos(lat * M_PI / 180.);
1143 sphi = sin(lat * M_PI / 180.);
1144
1145 return sqrt((gsl_pow_2(a * a * cphi) + gsl_pow_2(b * b * sphi))
1146 / (gsl_pow_2(a * cphi) + gsl_pow_2(b * sphi)));
1147}

◆ write_l1()

void write_l1 ( char *  filename,
iasi_l1_t *  l1 
)

Write IASI Level-1 data.

Definition at line 1151 of file libiasi.c.

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

1211 {
1212
1213 int dimid[10], ncid, time_id, z_id, lon_id, lat_id, p_id, t_id;
1214
1215 /* Create netCDF file... */
1216 LOG(1, "Write IASI Level-2 file: %s", filename);
1217 if (nc_open(filename, NC_WRITE, &ncid) != NC_NOERR) {
1218 NC(nc_create(filename, NC_CLOBBER, &ncid));
1219 } else {
1220 NC(nc_redef(ncid));
1221 }
1222
1223 /* Set dimensions... */
1224 if (nc_inq_dimid(ncid, "L2_NTRACK", &dimid[0]) != NC_NOERR)
1225 NC(nc_def_dim(ncid, "L2_NTRACK", l2->ntrack, &dimid[0]));
1226 if (nc_inq_dimid(ncid, "L2_NXTRACK", &dimid[1]) != NC_NOERR)
1227 NC(nc_def_dim(ncid, "L2_NXTRACK", L2_NXTRACK, &dimid[1]));
1228 if (nc_inq_dimid(ncid, "L2_NLAY", &dimid[2]) != NC_NOERR)
1229 NC(nc_def_dim(ncid, "L2_NLAY", L2_NLAY, &dimid[2]));
1230
1231 /* Add variables... */
1232 add_var(ncid, "l2_time", "s", "time (seconds since 2000-01-01T00:00Z)",
1233 NC_DOUBLE, dimid, &time_id, 2);
1234 add_var(ncid, "l2_z", "km", "altitude", NC_DOUBLE, dimid, &z_id, 3);
1235 add_var(ncid, "l2_lon", "deg", "longitude", NC_DOUBLE, dimid, &lon_id, 2);
1236 add_var(ncid, "l2_lat", "deg", "latitude", NC_DOUBLE, dimid, &lat_id, 2);
1237 add_var(ncid, "l2_press", "hPa", "pressure",
1238 NC_DOUBLE, &dimid[2], &p_id, 1);
1239 add_var(ncid, "l2_temp", "K", "temperature", NC_DOUBLE, dimid, &t_id, 3);
1240
1241 /* Leave define mode... */
1242 NC(nc_enddef(ncid));
1243
1244 /* Write data... */
1245 NC(nc_put_var_double(ncid, time_id, l2->time[0]));
1246 NC(nc_put_var_double(ncid, z_id, l2->z[0][0]));
1247 NC(nc_put_var_double(ncid, lon_id, l2->lon[0]));
1248 NC(nc_put_var_double(ncid, lat_id, l2->lat[0]));
1249 NC(nc_put_var_double(ncid, p_id, l2->p));
1250 NC(nc_put_var_double(ncid, t_id, l2->t[0][0]));
1251
1252 /* Close file... */
1253 NC(nc_close(ncid));
1254}
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