AIRS Code Collection
island.c
Go to the documentation of this file.
1/*
2 This file is part of the AIRS Code Collection.
3
4 the AIRS Code Collections is free software: you can redistribute it
5 and/or modify it under the terms of the GNU General Public License
6 as published by the Free Software Foundation, either version 3 of
7 the License, or (at your option) any later version.
8
9 The AIRS Code Collection is distributed in the hope that it will be
10 useful, but WITHOUT ANY WARRANTY; without even the implied warranty
11 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with the AIRS Code Collection. If not, see
16 <http://www.gnu.org/licenses/>.
17
18 Copyright (C) 2019-2025 Forschungszentrum Juelich GmbH
19*/
20
26#include "libairs.h"
27
28int main(
29 int argc,
30 char *argv[]) {
31
32 static pert_t *pert;
33
34 static wave_t *wave;
35
36 static FILE *in, *out;
37
38 static char pertname[LEN], ncfile[LEN];
39
40 static double ebt, emu, enoise, evar, wbt, wmu, wnoise, wvar, etime, wtime,
41 aux;
42
43 static int orb_old =
44 -1, en, wn, ncid, time_varid, track_varid, np_east_varid,
45 var_east_varid, np_west_varid, var_west_varid, year_varid, doy_varid,
46 track, year, mon, day, doy, iaux;
47
48 static size_t count[2] = { 1, 1 }, start[2];
49
50 /* Check arguments... */
51 if (argc < 4)
52 ERRMSG("Give parameters: <ctl> <var.tab> <pert1.nc> [<pert2.nc> ...]");
53
54 /* Get control parameters... */
55 scan_ctl(argc, argv, "PERTNAME", -1, "4mu", pertname);
56 const double lon0 = scan_ctl(argc, argv, "LON0", -1, "", NULL);
57 const double lat0 = scan_ctl(argc, argv, "LAT0", -1, "", NULL);
58 const double dlon = scan_ctl(argc, argv, "DLON", -1, "", NULL);
59 const double dlat = scan_ctl(argc, argv, "DLAT", -1, "", NULL);
60 const double offset = scan_ctl(argc, argv, "OFFSET", -1, "1", NULL);
61 const int bg_poly_x =
62 (int) scan_ctl(argc, argv, "BG_POLY_X", -1, "0", NULL);
63 const int bg_poly_y =
64 (int) scan_ctl(argc, argv, "BG_POLY_Y", -1, "0", NULL);
65 const int bg_smooth_x =
66 (int) scan_ctl(argc, argv, "BG_SMOOTH_X", -1, "0", NULL);
67 const int bg_smooth_y =
68 (int) scan_ctl(argc, argv, "BG_SMOOTH_Y", -1, "0", NULL);
69 const double gauss_fwhm = scan_ctl(argc, argv, "GAUSS_FWHM", -1, "0", NULL);
70 const double var_dh = scan_ctl(argc, argv, "VAR_DH", -1, "0", NULL);
71 const double orblat = scan_ctl(argc, argv, "ORBLAT", -1, "0", NULL);
72 const double dt230 = scan_ctl(argc, argv, "DT230", -1, "0.16", NULL);
73 const double nu = scan_ctl(argc, argv, "NU", -1, "2345.0", NULL);
74 scan_ctl(argc, argv, "NCFILE", -1, "-", ncfile);
75
76 /* Allocate... */
77 ALLOC(pert, pert_t, 1);
78
79 /* Create file... */
80 printf("Write variance statistics: %s\n", argv[2]);
81 if (!(out = fopen(argv[2], "w")))
82 ERRMSG("Cannot create file!");
83
84 /* Write header... */
85 fprintf(out,
86 "# $1 = time [s]\n"
87 "# $2 = orbit number\n"
88 "# $3 = eastern box: number of footprints\n"
89 "# $4 = eastern box: variance [K^2]\n"
90 "# $5 = eastern box: mean background temperature [K]\n"
91 "# $6 = eastern box: noise estimate [K]\n"
92 "# $7 = western box: number of footprints\n"
93 "# $8 = western box: variance [K^2]\n"
94 "# $9 = western box: mean background temperature [K]\n"
95 "# $10 = western box: noise estimate [K]\n\n");
96
97 /* Create netCDF file... */
98 if (ncfile[0] != '-') {
99
100 /* Create file... */
101 printf("Write variance statistics: %s\n", ncfile);
102 NC(nc_create(ncfile, NC_CLOBBER, &ncid));
103
104 /* Set dimensions... */
105 static int dimid[2];
106 NC(nc_def_dim(ncid, "NP", NC_UNLIMITED, &dimid[0]));
107
108 /* Add attributes... */
109 aux = lon0;
110 nc_put_att_double(ncid, NC_GLOBAL, "box_east_lon0", NC_DOUBLE, 1, &aux);
111 aux = lon0 + dlon;
112 nc_put_att_double(ncid, NC_GLOBAL, "box_east_lon1", NC_DOUBLE, 1, &aux);
113 aux = lat0 - 0.5 * dlat;
114 nc_put_att_double(ncid, NC_GLOBAL, "box_east_lat0", NC_DOUBLE, 1, &aux);
115 aux = lat0 + 0.5 * dlat;
116 nc_put_att_double(ncid, NC_GLOBAL, "box_east_lat1", NC_DOUBLE, 1, &aux);
117 aux = lon0 - dlon - offset;
118 nc_put_att_double(ncid, NC_GLOBAL, "box_west_lon0", NC_DOUBLE, 1, &aux);
119 aux = lon0 - offset;
120 nc_put_att_double(ncid, NC_GLOBAL, "box_west_lon1", NC_DOUBLE, 1, &aux);
121 aux = lat0 - 0.5 * dlat;
122 nc_put_att_double(ncid, NC_GLOBAL, "box_west_lat0", NC_DOUBLE, 1, &aux);
123 aux = lat0 + 0.5 * dlat;
124 nc_put_att_double(ncid, NC_GLOBAL, "box_west_lat1", NC_DOUBLE, 1, &aux);
125
126 /* Add variables... */
127 NC(nc_def_var(ncid, "time", NC_DOUBLE, 1, dimid, &time_varid));
128 add_att(ncid, time_varid, "s", "time (seconds since 2000-01-01T00:00Z)");
129 NC(nc_def_var(ncid, "year", NC_INT, 1, dimid, &year_varid));
130 add_att(ncid, year_varid, "1", "year");
131 NC(nc_def_var(ncid, "doy", NC_INT, 1, dimid, &doy_varid));
132 add_att(ncid, doy_varid, "1", "day of year");
133 NC(nc_def_var(ncid, "track", NC_INT, 1, dimid, &track_varid));
134 add_att(ncid, track_varid, "1", "along-track index");
135 NC(nc_def_var(ncid, "var_east", NC_DOUBLE, 1, dimid, &var_east_varid));
136 add_att(ncid, var_east_varid, "K^2", "BT variance (east)");
137 NC(nc_def_var(ncid, "var_west", NC_DOUBLE, 1, dimid, &var_west_varid));
138 add_att(ncid, var_west_varid, "K^2", "BT variance (west)");
139 NC(nc_def_var(ncid, "np_east", NC_INT, 1, dimid, &np_east_varid));
140 add_att(ncid, np_east_varid, "1", "number of footprints (east)");
141 NC(nc_def_var(ncid, "np_west", NC_INT, 1, dimid, &np_west_varid));
142 add_att(ncid, np_west_varid, "1", "number of footprints (west)");
143
144 /* Leave define mode... */
145 NC(nc_enddef(ncid));
146 }
147
148 /* Loop over perturbation files... */
149 for (int iarg = 3; iarg < argc; iarg++) {
150
151 /* Check filename... */
152 if (!strcmp(argv[iarg], ncfile))
153 continue;
154
155 /* Initialize... */
156 int orb = 0;
157
158 /* Read perturbation data... */
159 if (!(in = fopen(argv[iarg], "r")))
160 continue;
161 else {
162 fclose(in);
163 read_pert(argv[iarg], pertname, pert);
164 }
165
166 /* Recalculate background and perturbations... */
167 if (bg_poly_x > 0 || bg_poly_y > 0 ||
168 bg_smooth_x > 0 || bg_smooth_y > 0 || gauss_fwhm > 0 || var_dh > 0) {
169
170 /* Allocate... */
171 ALLOC(wave, wave_t, 1);
172
173 /* Convert to wave analysis struct... */
174 pert2wave(pert, wave, 0, pert->ntrack - 1, 0, pert->nxtrack - 1);
175
176 /* Estimate background... */
177 background_poly(wave, bg_poly_x, bg_poly_y);
178 background_smooth(wave, bg_smooth_x, bg_smooth_y);
179
180 /* Gaussian filter... */
181 gauss(wave, gauss_fwhm);
182
183 /* Compute variance... */
184 variance(wave, var_dh);
185
186 /* Copy data... */
187 for (int ix = 0; ix < wave->nx; ix++)
188 for (int iy = 0; iy < wave->ny; iy++) {
189 pert->pt[iy][ix] = wave->pt[ix][iy];
190 pert->var[iy][ix] = wave->var[ix][iy];
191 }
192
193 /* Free... */
194 free(wave);
195 }
196
197 /* Detection... */
198 for (int itrack = 0; itrack < pert->ntrack; itrack++)
199 for (int ixtrack = 0; ixtrack < pert->nxtrack; ixtrack++) {
200
201 /* Check data... */
202 if (pert->time[itrack][ixtrack] < 0
203 || pert->lon[itrack][ixtrack] < -180
204 || pert->lon[itrack][ixtrack] > 180
205 || pert->lat[itrack][ixtrack] < -90
206 || pert->lat[itrack][ixtrack] > 90
207 || pert->pt[itrack][ixtrack] < -100
208 || pert->pt[itrack][ixtrack] > 100
209 || !gsl_finite(pert->bt[itrack][ixtrack])
210 || !gsl_finite(pert->pt[itrack][ixtrack])
211 || !gsl_finite(pert->var[itrack][ixtrack])
212 || !gsl_finite(pert->dc[itrack][ixtrack]))
213 continue;
214
215 /* Count orbits... */
216 if (itrack > 0 && ixtrack == pert->nxtrack / 2)
217 if (pert->lat[itrack - 1][ixtrack] <= orblat
218 && pert->lat[itrack][ixtrack] >= orblat)
219 orb++;
220 if (orb != orb_old) {
221
222 /* Set orbit index... */
223 orb_old = orb;
224
225 /* Write output... */
226 if (en > 0 && wn > 0) {
227
228 /* Estimate noise... */
229 if (dt230 > 0) {
230 const double nesr =
231 PLANCK(230.0 + dt230, nu) - PLANCK(230.0, nu);
232 enoise = BRIGHT(PLANCK(ebt / en, nu) + nesr, nu) - ebt / en;
233 wnoise = BRIGHT(PLANCK(wbt / wn, nu) + nesr, nu) - wbt / wn;
234 }
235
236 /* Write output... */
237 fprintf(out, "%.2f %d %d %g %g %g %d %g %g %g\n", etime / en, orb,
238 en, evar / en - gsl_pow_2(emu / en), ebt / en, enoise,
239 wn, wvar / wn - gsl_pow_2(wmu / wn), wbt / wn, wnoise);
240
241 /* Write to netCDF file... */
242 if (ncfile[0] != '-') {
243
244 /* Get year and doy... */
245 jsec2time(etime / en, &year, &mon, &day, &iaux, &iaux, &iaux,
246 &aux);
247 day2doy(year, mon, day, &doy);
248
249 /* Find along-track index... */
250 track = 0;
251 for (int itrack2 = 0; itrack2 < pert->ntrack; itrack2++)
252 if (fabs(pert->time[itrack2][0] - etime / en)
253 < fabs(pert->time[track][0] - etime / en))
254 track = itrack2;
255
256 /* Write data... */
257 aux = etime / en;
258 NC(nc_put_vara_double(ncid, time_varid, start, count, &aux));
259 NC(nc_put_vara_int(ncid, year_varid, start, count, &year));
260 NC(nc_put_vara_int(ncid, doy_varid, start, count, &doy));
261 NC(nc_put_vara_int(ncid, track_varid, start, count, &track));
262 NC(nc_put_vara_int(ncid, np_east_varid, start, count, &en));
263 aux = evar / en - gsl_pow_2(emu / en) - gsl_pow_2(enoise);
264 NC(nc_put_vara_double
265 (ncid, var_east_varid, start, count, &aux));
266 NC(nc_put_vara_int(ncid, np_west_varid, start, count, &wn));
267 aux = wvar / wn - gsl_pow_2(wmu / wn) - gsl_pow_2(wnoise);
268 NC(nc_put_vara_double
269 (ncid, var_west_varid, start, count, &aux));
270
271 /* Increment data point counter... */
272 start[0]++;
273 }
274 }
275
276 /* Initialize... */
277 etime = wtime = 0;
278 evar = wvar = 0;
279 emu = wmu = 0;
280 ebt = wbt = 0;
281 en = wn = 0;
282 }
283
284 /* Check if footprint is in eastern box... */
285 if (pert->lon[itrack][ixtrack] >= lon0
286 && pert->lon[itrack][ixtrack] <= lon0 + dlon
287 && pert->lat[itrack][ixtrack] >= lat0 - dlat / 2.
288 && pert->lat[itrack][ixtrack] <= lat0 + dlat / 2.) {
289
290 etime += pert->time[itrack][ixtrack];
291 emu += pert->pt[itrack][ixtrack];
292 evar += gsl_pow_2(pert->pt[itrack][ixtrack]);
293 ebt += pert->bt[itrack][ixtrack];
294 en++;
295 }
296
297 /* Check if footprint is in western box... */
298 if (pert->lon[itrack][ixtrack] >= lon0 - offset - dlon
299 && pert->lon[itrack][ixtrack] <= lon0 - offset
300 && pert->lat[itrack][ixtrack] >= lat0 - dlat / 2.
301 && pert->lat[itrack][ixtrack] <= lat0 + dlat / 2.) {
302
303 wtime += pert->time[itrack][ixtrack];
304 wmu += pert->pt[itrack][ixtrack];
305 wvar += gsl_pow_2(pert->pt[itrack][ixtrack]);
306 wbt += pert->bt[itrack][ixtrack];
307 wn++;
308 }
309 }
310
311 /* Write output for last orbit... */
312 if (en > 0 && wn > 0) {
313
314 /* Estimate noise... */
315 if (dt230 > 0) {
316 const double nesr = PLANCK(230.0 + dt230, nu) - PLANCK(230.0, nu);
317 enoise = BRIGHT(PLANCK(ebt / en, nu) + nesr, nu) - ebt / en;
318 wnoise = BRIGHT(PLANCK(wbt / wn, nu) + nesr, nu) - wbt / wn;
319 }
320
321 /* Write output... */
322 fprintf(out, "%.2f %d %d %g %g %g %d %g %g %g\n", etime / en, orb,
323 en, evar / en - gsl_pow_2(emu / en), ebt / en, enoise,
324 wn, wvar / wn - gsl_pow_2(wmu / wn), wbt / wn, wnoise);
325
326 /* Write to netCDF file... */
327 if (ncfile[0] != '-') {
328
329 /* Get year and doy... */
330 jsec2time(etime / en, &year, &mon, &day, &iaux, &iaux, &iaux, &aux);
331 day2doy(year, mon, day, &doy);
332
333 /* Find along-track index... */
334 track = 0;
335 for (int itrack2 = 0; itrack2 < pert->ntrack; itrack2++)
336 if (fabs(pert->time[itrack2][0] - etime / en)
337 < fabs(pert->time[track][0] - etime / en))
338 track = itrack2;
339
340 /* Write data... */
341 aux = etime / en;
342 NC(nc_put_vara_double(ncid, time_varid, start, count, &aux));
343 NC(nc_put_vara_int(ncid, year_varid, start, count, &year));
344 NC(nc_put_vara_int(ncid, doy_varid, start, count, &doy));
345 NC(nc_put_vara_int(ncid, track_varid, start, count, &track));
346 NC(nc_put_vara_int(ncid, np_east_varid, start, count, &en));
347 aux = evar / en - gsl_pow_2(emu / en) - gsl_pow_2(enoise);
348 NC(nc_put_vara_double(ncid, var_east_varid, start, count, &aux));
349 NC(nc_put_vara_int(ncid, np_west_varid, start, count, &wn));
350 aux = wvar / wn - gsl_pow_2(wmu / wn) - gsl_pow_2(wnoise);
351 NC(nc_put_vara_double(ncid, var_west_varid, start, count, &aux));
352
353 /* Increment data point counter... */
354 start[0]++;
355 }
356 }
357 }
358
359 /* Close file... */
360 fclose(out);
361
362 /* Close file... */
363 if (ncfile[0] != '-')
364 NC(nc_close(ncid));
365
366 /* Free... */
367 free(pert);
368
369 return EXIT_SUCCESS;
370}
NC
#define NC(cmd)
Execute netCDF library command and check result.
Definition: diff_apr.c:35
main
int main(int argc, char *argv[])
Definition: island.c:28
jsec2time
void jsec2time(const double jsec, int *year, int *mon, int *day, int *hour, int *min, int *sec, double *remain)
Convert seconds to date.
Definition: jurassic.c:3971
scan_ctl
double scan_ctl(int argc, char *argv[], const char *varname, int arridx, const char *defvalue, char *value)
Search control parameter file for variable entry.
Definition: jurassic.c:5114
LEN
#define LEN
Maximum length of ASCII data lines.
Definition: jurassic.h:383
BRIGHT
#define BRIGHT(rad, nu)
Compute brightness temperature.
Definition: jurassic.h:126
ERRMSG
#define ERRMSG(...)
Print error message and quit program.
Definition: jurassic.h:237
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory.
Definition: jurassic.h:121
PLANCK
#define PLANCK(T, nu)
Compute Planck function.
Definition: jurassic.h:183
background_smooth
void background_smooth(wave_t *wave, int npts_x, int npts_y)
Smooth background.
Definition: libairs.c:176
add_att
void add_att(int ncid, int varid, const char *unit, const char *long_name)
Add variable attributes to netCDF file.
Definition: libairs.c:30
day2doy
void day2doy(int year, int mon, int day, int *doy)
Get day of year from date.
Definition: libairs.c:304
background_poly
void background_poly(wave_t *wave, int dim_x, int dim_y)
Get background based on polynomial fits.
Definition: libairs.c:126
variance
void variance(wave_t *wave, double dh)
Compute local variance.
Definition: libairs.c:1584
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: libairs.c:999
read_pert
void read_pert(char *filename, char *pertname, pert_t *pert)
Read radiance perturbation data.
Definition: libairs.c:1137
gauss
void gauss(wave_t *wave, double fwhm)
Apply Gaussian filter to perturbations...
Definition: libairs.c:579
libairs.h
AIRS Code Collection library declarations.
pert_t
Perturbation data.
Definition: libairs.h:205
pert_t::var
double var[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature variance (4 or 15 micron) [K].
Definition: libairs.h:232
pert_t::pt
double pt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature perturbation (4 or 15 micron) [K].
Definition: libairs.h:229
pert_t::time
double time[PERT_NTRACK][PERT_NXTRACK]
Time (seconds since 2000-01-01T00:00Z).
Definition: libairs.h:214
pert_t::dc
double dc[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (8 micron) [K].
Definition: libairs.h:223
pert_t::ntrack
int ntrack
Number of along-track values.
Definition: libairs.h:208
pert_t::nxtrack
int nxtrack
Number of across-track values.
Definition: libairs.h:211
pert_t::lon
double lon[PERT_NTRACK][PERT_NXTRACK]
Longitude [deg].
Definition: libairs.h:217
pert_t::bt
double bt[PERT_NTRACK][PERT_NXTRACK]
Brightness temperature (4 or 15 micron) [K].
Definition: libairs.h:226
pert_t::lat
double lat[PERT_NTRACK][PERT_NXTRACK]
Latitude [deg].
Definition: libairs.h:220
wave_t
Wave analysis data.
Definition: libairs.h:287
wave_t::nx
int nx
Number of across-track values.
Definition: libairs.h:290
wave_t::ny
int ny
Number of along-track values.
Definition: libairs.h:293
wave_t::var
double var[WX][WY]
Variance [K].
Definition: libairs.h:323
wave_t::pt
double pt[WX][WY]
Perturbation [K].
Definition: libairs.h:320
Generated on Wed Feb 5 2025 15:19:54 for AIRS Code Collection by doxygen 1.9.4