IASI Code Collection
Macros | Functions
perturbation.c File Reference

Extract perturbation data. More...

#include "libiasi.h"

Go to the source code of this file.

Macros

#define N4   152
 
#define N15_LOW   21
 
#define N15_HIGH   2
 

Functions

void addatt (int ncid, int varid, const char *unit, const char *long_name)
 
int main (int argc, char *argv[])
 

Detailed Description

Extract perturbation data.

Definition in file perturbation.c.

Macro Definition Documentation

◆ N4

#define N4   152

Definition at line 33 of file perturbation.c.

◆ N15_LOW

#define N15_LOW   21

Definition at line 36 of file perturbation.c.

◆ N15_HIGH

#define N15_HIGH   2

Definition at line 39 of file perturbation.c.

Function Documentation

◆ addatt()

void addatt ( int  ncid,
int  varid,
const char *  unit,
const char *  long_name 
)

Definition at line 411 of file perturbation.c.

415 {
416
417 /* Set long name... */
418 NC(nc_put_att_text(ncid, varid, "long_name", strlen(long_name), long_name));
419
420 /* Set units... */
421 NC(nc_put_att_text(ncid, varid, "units", strlen(unit), unit));
422}
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: libiasi.h:155

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 56 of file perturbation.c.

58 {
59
60 static iasi_rad_t *iasi_rad;
61
62 static pert_t *pert_4mu, *pert_15mu_low, *pert_15mu_high;
63
64 static wave_t wave;
65
66 static double numean, radmean, var_dh = 100.;
67
68 static int list_4mu[N4]
69 = { 6711, 6712, 6713, 6714, 6715, 6716, 6717, 6718, 6719, 6720,
70 6721, 6722, 6723, 6724, 6725, 6726, 6727, 6728, 6729, 6730, 6731,
71 6732, 6733, 6734, 6735, 6736, 6737, 6738, 6739, 6740, 6741, 6742,
72 6743, 6744, 6745, 6746, 6747, 6748, 6749, 6750, 6751, 6752, 6753,
73 6754, 6755, 6756, 6757, 6758, 6759, 6760, 6761, 6762, 6763, 6764,
74 6765, 6766, 6767, 6768, 6769, 6770, 6771, 6772, 6773, 6774, 6775,
75 6776, 6777, 6778, 6779, 6780, 6781, 6782, 6783, 6784, 6785, 6786,
76 6787, 6788, 6789, 6790, 6791, 6792, 6793, 6794, 6795, 6796, 6797,
77 6798, 6799, 6800, 6801, 6802, 6803, 6804, 6830, 6831, 6832, 6833,
78 6834, 6835, 6836, 6837, 6838, 6839, 6840, 6841, 6842, 6843, 6844,
79 6845, 6846, 6847, 6848, 6849, 6850, 6851, 6852, 6853, 6854, 6855,
80 6856, 6857, 6858, 6859, 6860, 6861, 6862, 6863, 6864, 6865, 6866,
81 6867, 6868, 6869, 6870, 6871, 6872, 6873, 6874, 6875, 6876, 6877,
82 6878, 6879, 6880, 6881, 6882, 6883, 6884, 6885, 6886, 6887
83 };
84
85 static int list_15mu_low[N15_LOW]
86 = { 22, 28, 34, 40, 46, 52, 58, 72, 100, 105, 112, 118, 119,
87 124, 125, 130, 131, 136, 137, 143, 144
88 };
89
90 static int list_15mu_high[N15_HIGH]
91 = { 91, 92 };
92
93 static int ix, iy, dimid[2], i, n, ncid, track, track0, xtrack,
94 time_varid, lon_varid, lat_varid, bt_4mu_varid, bt_4mu_pt_varid,
95 bt_4mu_var_varid, bt_8mu_varid, bt_15mu_low_varid, bt_15mu_low_pt_varid,
96 bt_15mu_low_var_varid, bt_15mu_high_varid, bt_15mu_high_pt_varid,
97 bt_15mu_high_var_varid, iarg, format;
98
99 static size_t start[2], count[2];
100
101 /* Check arguments... */
102 if (argc < 4)
103 ERRMSG("Give parameters: <ctl> <out.nc> <l1b_file1> [<l1b_file2> ...]");
104
105 /* Read control parameters... */
106 format = (int) scan_ctl(argc, argv, "FORMAT", -1, "1", NULL);
107
108 /* Allocate... */
109 ALLOC(iasi_rad, iasi_rad_t, 1);
110 ALLOC(pert_4mu, pert_t, 1);
111 ALLOC(pert_15mu_low, pert_t, 1);
112 ALLOC(pert_15mu_high, pert_t, 1);
113
114 /* ------------------------------------------------------------
115 Read HDF files...
116 ------------------------------------------------------------ */
117
118 /* Loop over HDF files... */
119 for (iarg = 3; iarg < argc; iarg++) {
120
121 /* Read IASI data... */
122 printf("Read IASI Level-1C data file: %s\n", argv[iarg]);
123 iasi_read(format, argv[iarg], iasi_rad);
124
125 /* Save geolocation... */
126 pert_4mu->ntrack += iasi_rad->ntrack;
127 if (pert_4mu->ntrack > PERT_NTRACK)
128 ERRMSG("Too many granules!");
129 pert_4mu->nxtrack = L1_NXTRACK;
130 if (pert_4mu->nxtrack > PERT_NXTRACK)
131 ERRMSG("Too many tracks!");
132 for (track = 0; track < iasi_rad->ntrack; track++)
133 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
134 pert_4mu->time[track0 + track][xtrack]
135 = iasi_rad->Time[track][xtrack];
136 pert_4mu->lon[track0 + track][xtrack]
137 = iasi_rad->Longitude[track][xtrack];
138 pert_4mu->lat[track0 + track][xtrack]
139 = iasi_rad->Latitude[track][xtrack];
140 }
141
142 pert_15mu_low->ntrack += iasi_rad->ntrack;
143 if (pert_15mu_low->ntrack > PERT_NTRACK)
144 ERRMSG("Too many granules!");
145 pert_15mu_low->nxtrack = L1_NXTRACK;
146 if (pert_15mu_low->nxtrack > PERT_NXTRACK)
147 ERRMSG("Too many tracks!");
148 for (track = 0; track < iasi_rad->ntrack; track++)
149 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
150 pert_15mu_low->time[track0 + track][xtrack]
151 = iasi_rad->Time[track][xtrack];
152 pert_15mu_low->lon[track0 + track][xtrack]
153 = iasi_rad->Longitude[track][xtrack];
154 pert_15mu_low->lat[track0 + track][xtrack]
155 = iasi_rad->Latitude[track][xtrack];
156 }
157
158 pert_15mu_high->ntrack += iasi_rad->ntrack;
159 if (pert_15mu_high->ntrack > PERT_NTRACK)
160 ERRMSG("Too many granules!");
161 pert_15mu_high->nxtrack = L1_NXTRACK;
162 if (pert_15mu_high->nxtrack > PERT_NXTRACK)
163 ERRMSG("Too many tracks!");
164 for (track = 0; track < iasi_rad->ntrack; track++)
165 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
166 pert_15mu_high->time[track0 + track][xtrack]
167 = iasi_rad->Time[track][xtrack];
168 pert_15mu_high->lon[track0 + track][xtrack]
169 = iasi_rad->Longitude[track][xtrack];
170 pert_15mu_high->lat[track0 + track][xtrack]
171 = iasi_rad->Latitude[track][xtrack];
172 }
173
174 /* Get 8.1 micron brightness temperature... */
175 for (track = 0; track < iasi_rad->ntrack; track++)
176 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++)
177 pert_4mu->dc[track0 + track][xtrack]
178 = BRIGHT(iasi_rad->Rad[track][xtrack][2345], iasi_rad->freq[2345]);
179
180 /* Get 4.3 micron brightness temperature... */
181 for (track = 0; track < iasi_rad->ntrack; track++)
182 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
183 n = 0;
184 numean = radmean = 0;
185 for (i = 0; i < N4; i++)
186 if (gsl_finite(iasi_rad->Rad[track][xtrack][list_4mu[i]])) {
187 radmean += iasi_rad->Rad[track][xtrack][list_4mu[i]];
188 numean += iasi_rad->freq[list_4mu[i]];
189 n++;
190 }
191 if (n > 0.9 * N4)
192 pert_4mu->bt[track0 + track][xtrack]
193 = BRIGHT(radmean / n, numean / n);
194 else
195 pert_4mu->bt[track0 + track][xtrack] = GSL_NAN;
196 }
197
198 /* Get 15 micron brightness temperature (low altitudes)... */
199 for (track = 0; track < iasi_rad->ntrack; track++)
200 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
201 n = 0;
202 numean = radmean = 0;
203 for (i = 0; i < N15_LOW; i++)
204 if (gsl_finite(iasi_rad->Rad[track][xtrack][list_15mu_low[i]])) {
205 radmean += iasi_rad->Rad[track][xtrack][list_15mu_low[i]];
206 numean += iasi_rad->freq[list_15mu_low[i]];
207 n++;
208 }
209 if (n > 0.9 * N15_LOW)
210 pert_15mu_low->bt[track0 + track][xtrack]
211 = BRIGHT(radmean / n, numean / n);
212 else
213 pert_15mu_low->bt[track0 + track][xtrack] = GSL_NAN;
214 }
215
216 /* Get 15 micron brightness temperature (high altitudes)... */
217 for (track = 0; track < iasi_rad->ntrack; track++)
218 for (xtrack = 0; xtrack < L1_NXTRACK; xtrack++) {
219 n = 0;
220 numean = radmean = 0;
221 for (i = 0; i < N15_HIGH; i++)
222 if (gsl_finite(iasi_rad->Rad[track][xtrack][list_15mu_high[i]])) {
223 radmean += iasi_rad->Rad[track][xtrack][list_15mu_high[i]];
224 numean += iasi_rad->freq[list_15mu_high[i]];
225 n++;
226 }
227 if (n > 0.9 * N15_HIGH)
228 pert_15mu_high->bt[track0 + track][xtrack]
229 = BRIGHT(radmean / n, numean / n);
230 else
231 pert_15mu_high->bt[track0 + track][xtrack] = GSL_NAN;
232 }
233
234 /* Increment track counter... */
235 track0 += iasi_rad->ntrack;
236 }
237
238 /* ------------------------------------------------------------
239 Calculate perturbations and variances...
240 ------------------------------------------------------------ */
241
242 /* Convert to wave analysis struct... */
243 pert2wave(pert_4mu, &wave,
244 0, pert_4mu->ntrack - 1, 0, pert_4mu->nxtrack - 1);
245
246 /* Estimate background... */
247 background_poly(&wave, 5, 0);
248
249 /* Compute variance... */
250 variance(&wave, var_dh);
251
252 /* Copy data... */
253 for (ix = 0; ix < wave.nx; ix++)
254 for (iy = 0; iy < wave.ny; iy++) {
255 pert_4mu->pt[iy][ix] = wave.pt[ix][iy];
256 pert_4mu->var[iy][ix] = wave.var[ix][iy];
257 }
258
259 /* Convert to wave analysis struct... */
260 pert2wave(pert_15mu_low, &wave,
261 0, pert_15mu_low->ntrack - 1, 0, pert_15mu_low->nxtrack - 1);
262
263 /* Estimate background... */
264 background_poly(&wave, 5, 0);
265
266 /* Compute variance... */
267 variance(&wave, var_dh);
268
269 /* Copy data... */
270 for (ix = 0; ix < wave.nx; ix++)
271 for (iy = 0; iy < wave.ny; iy++) {
272 pert_15mu_low->pt[iy][ix] = wave.pt[ix][iy];
273 pert_15mu_low->var[iy][ix] = wave.var[ix][iy];
274 }
275
276 /* Convert to wave analysis struct... */
277 pert2wave(pert_15mu_high, &wave,
278 0, pert_15mu_high->ntrack - 1, 0, pert_15mu_high->nxtrack - 1);
279
280 /* Estimate background... */
281 background_poly(&wave, 5, 0);
282
283 /* Compute variance... */
284 variance(&wave, var_dh);
285
286 /* Copy data... */
287 for (ix = 0; ix < wave.nx; ix++)
288 for (iy = 0; iy < wave.ny; iy++) {
289 pert_15mu_high->pt[iy][ix] = wave.pt[ix][iy];
290 pert_15mu_high->var[iy][ix] = wave.var[ix][iy];
291 }
292
293 /* ------------------------------------------------------------
294 Write to netCDF file...
295 ------------------------------------------------------------ */
296
297 /* Create netCDF file... */
298 printf("Write perturbation data file: %s\n", argv[2]);
299 NC(nc_create(argv[2], NC_CLOBBER, &ncid));
300
301 /* Set dimensions... */
302 NC(nc_def_dim(ncid, "NTRACK", NC_UNLIMITED, &dimid[0]));
303 NC(nc_def_dim(ncid, "NXTRACK", L1_NXTRACK, &dimid[1]));
304
305 /* Add variables... */
306 NC(nc_def_var(ncid, "time", NC_DOUBLE, 2, dimid, &time_varid));
307 addatt(ncid, time_varid, "s", "time (seconds since 2000-01-01T00:00Z)");
308 NC(nc_def_var(ncid, "lon", NC_DOUBLE, 2, dimid, &lon_varid));
309 addatt(ncid, lon_varid, "deg", "footprint longitude");
310 NC(nc_def_var(ncid, "lat", NC_DOUBLE, 2, dimid, &lat_varid));
311 addatt(ncid, lat_varid, "deg", "footprint latitude");
312
313 NC(nc_def_var(ncid, "bt_8mu", NC_FLOAT, 2, dimid, &bt_8mu_varid));
314 addatt(ncid, bt_8mu_varid, "K", "brightness temperature at 8.1 micron");
315
316 NC(nc_def_var(ncid, "bt_4mu", NC_FLOAT, 2, dimid, &bt_4mu_varid));
317 addatt(ncid, bt_4mu_varid, "K", "brightness temperature" " at 4.3 micron");
318 NC(nc_def_var(ncid, "bt_4mu_pt", NC_FLOAT, 2, dimid, &bt_4mu_pt_varid));
319 addatt(ncid, bt_4mu_pt_varid, "K", "brightness temperature perturbation"
320 " at 4.3 micron");
321 NC(nc_def_var(ncid, "bt_4mu_var", NC_FLOAT, 2, dimid, &bt_4mu_var_varid));
322 addatt(ncid, bt_4mu_var_varid, "K^2", "brightness temperature variance"
323 " at 4.3 micron");
324
325 NC(nc_def_var(ncid, "bt_15mu_low", NC_FLOAT, 2, dimid, &bt_15mu_low_varid));
326 addatt(ncid, bt_15mu_low_varid, "K", "brightness temperature"
327 " at 15 micron (low altitudes)");
328 NC(nc_def_var(ncid, "bt_15mu_low_pt", NC_FLOAT, 2, dimid,
329 &bt_15mu_low_pt_varid));
330 addatt(ncid, bt_15mu_low_pt_varid, "K",
331 "brightness temperature perturbation"
332 " at 15 micron (low altitudes)");
333 NC(nc_def_var
334 (ncid, "bt_15mu_low_var", NC_FLOAT, 2, dimid, &bt_15mu_low_var_varid));
335 addatt(ncid, bt_15mu_low_var_varid, "K^2",
336 "brightness temperature variance" " at 15 micron (low altitudes)");
337
338 NC(nc_def_var(ncid, "bt_15mu_high", NC_FLOAT, 2, dimid,
339 &bt_15mu_high_varid));
340 addatt(ncid, bt_15mu_high_varid, "K", "brightness temperature"
341 " at 15 micron (high altitudes)");
342 NC(nc_def_var(ncid, "bt_15mu_high_pt", NC_FLOAT, 2, dimid,
343 &bt_15mu_high_pt_varid));
344 addatt(ncid, bt_15mu_high_pt_varid, "K",
345 "brightness temperature perturbation"
346 " at 15 micron (high altitudes)");
347 NC(nc_def_var
348 (ncid, "bt_15mu_high_var", NC_FLOAT, 2, dimid, &bt_15mu_high_var_varid));
349 addatt(ncid, bt_15mu_high_var_varid, "K^2",
350 "brightness temperature variance" " at 15 micron (high altitudes)");
351
352 /* Leave define mode... */
353 NC(nc_enddef(ncid));
354
355 /* Loop over tracks... */
356 for (track = 0; track < pert_4mu->ntrack; track++) {
357
358 /* Set array sizes... */
359 start[0] = (size_t) track;
360 start[1] = 0;
361 count[0] = 1;
362 count[1] = (size_t) pert_4mu->nxtrack;
363
364 /* Write data... */
365 NC(nc_put_vara_double(ncid, time_varid, start, count,
366 pert_4mu->time[track]));
367 NC(nc_put_vara_double(ncid, lon_varid, start, count,
368 pert_4mu->lon[track]));
369 NC(nc_put_vara_double(ncid, lat_varid, start, count,
370 pert_4mu->lat[track]));
371
372 NC(nc_put_vara_double(ncid, bt_8mu_varid, start, count,
373 pert_4mu->dc[track]));
374
375 NC(nc_put_vara_double(ncid, bt_4mu_varid, start, count,
376 pert_4mu->bt[track]));
377 NC(nc_put_vara_double(ncid, bt_4mu_pt_varid, start, count,
378 pert_4mu->pt[track]));
379 NC(nc_put_vara_double(ncid, bt_4mu_var_varid, start, count,
380 pert_4mu->var[track]));
381
382 NC(nc_put_vara_double(ncid, bt_15mu_low_varid, start, count,
383 pert_15mu_low->bt[track]));
384 NC(nc_put_vara_double(ncid, bt_15mu_low_pt_varid, start, count,
385 pert_15mu_low->pt[track]));
386 NC(nc_put_vara_double(ncid, bt_15mu_low_var_varid, start, count,
387 pert_15mu_low->var[track]));
388
389 NC(nc_put_vara_double(ncid, bt_15mu_high_varid, start, count,
390 pert_15mu_high->bt[track]));
391 NC(nc_put_vara_double(ncid, bt_15mu_high_pt_varid, start, count,
392 pert_15mu_high->pt[track]));
393 NC(nc_put_vara_double(ncid, bt_15mu_high_var_varid, start, count,
394 pert_15mu_high->var[track]));
395 }
396
397 /* Close file... */
398 NC(nc_close(ncid));
399
400 /* Free... */
401 free(iasi_rad);
402 free(pert_4mu);
403 free(pert_15mu_low);
404 free(pert_15mu_high);
405
406 return EXIT_SUCCESS;
407}
iasi_read
void iasi_read(int format, char *filename, iasi_rad_t *iasi_rad)
Read IASI Level-1 data.
Definition: libiasi.c:297
background_poly
void background_poly(wave_t *wave, int dim_x, int dim_y)
Get background based on polynomial fits.
Definition: libiasi.c:114
variance
void variance(wave_t *wave, double dh)
Compute local variance.
Definition: libiasi.c:1011
pert2wave
void pert2wave(pert_t *pert, wave_t *wave, int track0, int track1, int xtrack0, int xtrack1)
Convert radiance perturbation data to wave analysis struct.
Definition: libiasi.c:853
PERT_NXTRACK
#define PERT_NXTRACK
Across-track size of perturbation data.
Definition: libiasi.h:135
L1_NXTRACK
#define L1_NXTRACK
Across-track size of IASI radiance granule (don't change).
Definition: libiasi.h:102
PERT_NTRACK
#define PERT_NTRACK
Along-track size of perturbation data.
Definition: libiasi.h:132
N15_LOW
#define N15_LOW
Definition: perturbation.c:36
N4
#define N4
Definition: perturbation.c:33
addatt
void addatt(int ncid, int varid, const char *unit, const char *long_name)
Definition: perturbation.c:411
N15_HIGH
#define N15_HIGH
Definition: perturbation.c:39
iasi_rad_t
IASI converted Level-1 radiation data.
Definition: libiasi.h:288
iasi_rad_t::Longitude
double Longitude[L1_NTRACK][L1_NXTRACK]
Longitude of the sounder pixel.
Definition: libiasi.h:300
iasi_rad_t::Latitude
double Latitude[L1_NTRACK][L1_NXTRACK]
Latitude of the sounder pixel.
Definition: libiasi.h:303
iasi_rad_t::Time
double Time[L1_NTRACK][L1_NXTRACK]
Seconds since 2000-01-01 for each sounder pixel.
Definition: libiasi.h:297
iasi_rad_t::freq
double freq[IASI_L1_NCHAN]
channel wavenumber [cm^-1]
Definition: libiasi.h:294
iasi_rad_t::ntrack
int ntrack
Number of along-track samples.
Definition: libiasi.h:291
iasi_rad_t::Rad
float Rad[L1_NTRACK][L1_NXTRACK][IASI_L1_NCHAN]
Radiance [W/(m^2 sr cm^-1)].
Definition: libiasi.h:306
pert_t
Perturbation data.
Definition: libiasi.h:224
pert_t::var
double var[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature variance (4 or 15 micron) [K].
Definition: libiasi.h:251
pert_t::pt
double pt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature perturbation (4 or 15 micron) [K].
Definition: libiasi.h:248
pert_t::time
double time[PERT_NTRACK][PERT_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libiasi.h:233
pert_t::dc
double dc[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (8 micron) [K].
Definition: libiasi.h:242
pert_t::ntrack
int ntrack
Number of along-track values.
Definition: libiasi.h:227
pert_t::nxtrack
int nxtrack
Number of across-track values.
Definition: libiasi.h:230
pert_t::lon
double lon[PERT_NTRACK][PERT_NXTRACK]
Longitude [deg].
Definition: libiasi.h:236
pert_t::bt
double bt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (4 or 15 micron) [K].
Definition: libiasi.h:245
pert_t::lat
double lat[PERT_NTRACK][PERT_NXTRACK]
Latitude [deg].
Definition: libiasi.h:239
wave_t
Wave analysis data.
Definition: libiasi.h:320
wave_t::nx
int nx
Number of across-track values.
Definition: libiasi.h:323
wave_t::ny
int ny
Number of along-track values.
Definition: libiasi.h:326
wave_t::var
double var[WX][WY]
Variance [K].
Definition: libiasi.h:356
wave_t::pt
double pt[WX][WY]
Perturbation [K].
Definition: libiasi.h:353
Here is the call graph for this function:
Generated by doxygen 1.9.4