MPTRAC
atm_dist.c
Go to the documentation of this file.
1/*
2 This file is part of MPTRAC.
3
4 MPTRAC is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
8
9 MPTRAC is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with MPTRAC. If not, see <http://www.gnu.org/licenses/>.
16
17 Copyright (C) 2013-2026 Forschungszentrum Juelich GmbH
18*/
19
25#include "mptrac.h"
26
27/* ------------------------------------------------------------
28 Functions...
29 ------------------------------------------------------------ */
30
32void usage(
33 void);
34
35/* ------------------------------------------------------------
36 Main...
37 ------------------------------------------------------------ */
38
39int main(
40 int argc,
41 char *argv[]) {
42
43 ctl_t ctl;
44
45 atm_t *atm1, *atm2;
46
47 FILE *out;
48
49 double *ahtd, *aqtd, *avtd, ahtdm, aqtdm[NQ], avtdm, *lat1_old, *lat2_old,
50 *lh1, *lh2, *lon1_old, *lon2_old, *lv1, *lv2, *rhtd, *rqtd, *rvtd, rhtdm,
51 rqtdm[NQ], rvtdm, t0 =
52 0, x0[3], x1[3], x2[3], z1, *z1_old, z2, *z2_old, *work;
53
54 int init = 0, np;
55
56 /* Allocate... */
57 ALLOC(atm1, atm_t, 1);
58 ALLOC(atm2, atm_t, 1);
59 ALLOC(lon1_old, double,
60 NP);
61 ALLOC(lat1_old, double,
62 NP);
63 ALLOC(z1_old, double,
64 NP);
65 ALLOC(lh1, double,
66 NP);
67 ALLOC(lv1, double,
68 NP);
69 ALLOC(lon2_old, double,
70 NP);
71 ALLOC(lat2_old, double,
72 NP);
73 ALLOC(z2_old, double,
74 NP);
75 ALLOC(lh2, double,
76 NP);
77 ALLOC(lv2, double,
78 NP);
79 ALLOC(ahtd, double,
80 NP);
81 ALLOC(avtd, double,
82 NP);
83 ALLOC(aqtd, double,
84 NP * NQ);
85 ALLOC(rhtd, double,
86 NP);
87 ALLOC(rvtd, double,
88 NP);
89 ALLOC(rqtd, double,
90 NP * NQ);
91 ALLOC(work, double,
92 NP);
93
94 /* Print usage information... */
95 USAGE;
96
97 /* Check arguments... */
98 if (argc < 6)
99 ERRMSG("Missing or invalid command-line arguments.\n\n"
100 "Usage: atm_dist <ctl> <dist.tab> <param> <atm1a> <atm1b> [<atm2a> <atm2b> ...]\n\n"
101 "Use -h for full help.");
102
103 /* Read control parameters... */
104 mptrac_read_ctl(argv[1], argc, argv, &ctl);
105 const int ens =
106 (int) scan_ctl(argv[1], argc, argv, "DIST_ENS", -1, "-999", NULL);
107 const double p0 =
108 P(scan_ctl(argv[1], argc, argv, "DIST_Z0", -1, "-1000", NULL));
109 const double p1 =
110 P(scan_ctl(argv[1], argc, argv, "DIST_Z1", -1, "1000", NULL));
111 const double lat0 =
112 scan_ctl(argv[1], argc, argv, "DIST_LAT0", -1, "-1000", NULL);
113 const double lat1 =
114 scan_ctl(argv[1], argc, argv, "DIST_LAT1", -1, "1000", NULL);
115 const double lon0 =
116 scan_ctl(argv[1], argc, argv, "DIST_LON0", -1, "-1000", NULL);
117 const double lon1 =
118 scan_ctl(argv[1], argc, argv, "DIST_LON1", -1, "1000", NULL);
119 const double zscore =
120 scan_ctl(argv[1], argc, argv, "DIST_ZSCORE", -1, "-999", NULL);
121
122 /* Write info... */
123 LOG(1, "Write transport deviations: %s", argv[2]);
124
125 /* Create output file... */
126 if (!(out = fopen(argv[2], "w")))
127 ERRMSG("Cannot create file!");
128
129 /* Write header... */
130 fprintf(out,
131 "# $1 = time [s]\n"
132 "# $2 = time difference [s]\n"
133 "# $3 = absolute horizontal distance (%s) [km]\n"
134 "# $4 = relative horizontal distance (%s) [%%]\n"
135 "# $5 = absolute vertical distance (%s) [km]\n"
136 "# $6 = relative vertical distance (%s) [%%]\n",
137 argv[3], argv[3], argv[3], argv[3]);
138 for (int iq = 0; iq < ctl.nq; iq++)
139 fprintf(out,
140 "# $%d = %s absolute difference (%s) [%s]\n"
141 "# $%d = %s relative difference (%s) [%%]\n",
142 7 + 2 * iq, ctl.qnt_name[iq], argv[3], ctl.qnt_unit[iq],
143 8 + 2 * iq, ctl.qnt_name[iq], argv[3]);
144 fprintf(out, "# $%d = number of particles\n\n", 7 + 2 * ctl.nq);
145
146 /* Loop over file pairs... */
147 for (int f = 4; f < argc; f += 2) {
148
149 /* Read atmopheric data... */
150 if (!mptrac_read_atm(argv[f], &ctl, atm1)
151 || !mptrac_read_atm(argv[f + 1], &ctl, atm2))
152 continue;
153
154 /* Check if structs match... */
155 if (atm1->np != atm2->np)
156 ERRMSG("Different numbers of particles!");
157
158 /* Get time from filename... */
159 const double t = time_from_filename(argv[f], ctl.atm_type < 2 ? 20 : 19);
160
161 /* Save initial time... */
162 if (!init) {
163 init = 1;
164 t0 = t;
165 }
166
167 /* Init... */
168 np = 0;
169 for (int ip = 0; ip < atm1->np; ip++) {
170 ahtd[ip] = avtd[ip] = rhtd[ip] = rvtd[ip] = 0;
171 for (int iq = 0; iq < ctl.nq; iq++)
172 aqtd[iq * NP + ip] = rqtd[iq * NP + ip] = 0;
173 }
174
175 /* Loop over air parcels... */
176 for (int ip = 0; ip < atm1->np; ip++) {
177
178 /* Check air parcel index... */
179 if (ctl.qnt_idx > 0
180 && (atm1->q[ctl.qnt_idx][ip] != atm2->q[ctl.qnt_idx][ip]))
181 ERRMSG("Air parcel index does not match!");
182
183 /* Check ensemble index... */
184 if (ctl.qnt_ens > 0
185 && (atm1->q[ctl.qnt_ens][ip] != ens
186 || atm2->q[ctl.qnt_ens][ip] != ens))
187 continue;
188
189 /* Check time... */
190 if (!isfinite(atm1->time[ip]) || !isfinite(atm2->time[ip]))
191 continue;
192
193 /* Check spatial range... */
194 if (atm1->p[ip] > p0 || atm1->p[ip] < p1
195 || atm1->lon[ip] < lon0 || atm1->lon[ip] > lon1
196 || atm1->lat[ip] < lat0 || atm1->lat[ip] > lat1)
197 continue;
198 if (atm2->p[ip] > p0 || atm2->p[ip] < p1
199 || atm2->lon[ip] < lon0 || atm2->lon[ip] > lon1
200 || atm2->lat[ip] < lat0 || atm2->lat[ip] > lat1)
201 continue;
202
203 /* Convert coordinates... */
204 geo2cart(0, atm1->lon[ip], atm1->lat[ip], x1);
205 geo2cart(0, atm2->lon[ip], atm2->lat[ip], x2);
206 z1 = Z(atm1->p[ip]);
207 z2 = Z(atm2->p[ip]);
208
209 /* Calculate absolute transport deviations... */
210 ahtd[np] = DIST(x1, x2);
211 avtd[np] = z1 - z2;
212 for (int iq = 0; iq < ctl.nq; iq++)
213 aqtd[iq * NP + np] = atm1->q[iq][ip] - atm2->q[iq][ip];
214
215 /* Calculate relative transport deviations... */
216 if (f > 4) {
217
218 /* Get trajectory lengths... */
219 geo2cart(0, lon1_old[ip], lat1_old[ip], x0);
220 lh1[ip] += DIST(x0, x1);
221 lv1[ip] += fabs(z1_old[ip] - z1);
222
223 geo2cart(0, lon2_old[ip], lat2_old[ip], x0);
224 lh2[ip] += DIST(x0, x2);
225 lv2[ip] += fabs(z2_old[ip] - z2);
226
227 /* Get relative transport deviations... */
228 if (lh1[ip] + lh2[ip] > 0)
229 rhtd[np] = 200. * DIST(x1, x2) / (lh1[ip] + lh2[ip]);
230 if (lv1[ip] + lv2[ip] > 0)
231 rvtd[np] = 200. * (z1 - z2) / (lv1[ip] + lv2[ip]);
232 }
233
234 /* Get relative transport deviations... */
235 for (int iq = 0; iq < ctl.nq; iq++)
236 rqtd[iq * NP + np] = 200. * (atm1->q[iq][ip] - atm2->q[iq][ip])
237 / (fabs(atm1->q[iq][ip]) + fabs(atm2->q[iq][ip]));
238
239 /* Save positions of air parcels... */
240 lon1_old[ip] = atm1->lon[ip];
241 lat1_old[ip] = atm1->lat[ip];
242 z1_old[ip] = z1;
243
244 lon2_old[ip] = atm2->lon[ip];
245 lat2_old[ip] = atm2->lat[ip];
246 z2_old[ip] = z2;
247
248 /* Increment air parcel counter... */
249 np++;
250 }
251
252 /* Filter data... */
253 if (zscore > 0 && np > 1) {
254
255 /* Get means and standard deviations of transport deviations... */
256 const size_t n = (size_t) np;
257 const double muh = gsl_stats_mean(ahtd, 1, n);
258 const double muv = gsl_stats_mean(avtd, 1, n);
259 const double sigh = gsl_stats_sd(ahtd, 1, n);
260 const double sigv = gsl_stats_sd(avtd, 1, n);
261
262 /* Filter data... */
263 np = 0;
264 for (size_t i = 0; i < n; i++)
265 if (fabs((ahtd[i] - muh) / sigh) < zscore
266 && fabs((avtd[i] - muv) / sigv) < zscore) {
267 ahtd[np] = ahtd[i];
268 rhtd[np] = rhtd[i];
269 avtd[np] = avtd[i];
270 rvtd[np] = rvtd[i];
271 for (int iq = 0; iq < ctl.nq; iq++) {
272 aqtd[iq * NP + np] = aqtd[iq * NP + (int) i];
273 rqtd[iq * NP + np] = rqtd[iq * NP + (int) i];
274 }
275 np++;
276 }
277 }
278
279 /* Get statistics... */
280 if (strcasecmp(argv[3], "mean") == 0) {
281 ahtdm = gsl_stats_mean(ahtd, 1, (size_t) np);
282 rhtdm = gsl_stats_mean(rhtd, 1, (size_t) np);
283 avtdm = gsl_stats_mean(avtd, 1, (size_t) np);
284 rvtdm = gsl_stats_mean(rvtd, 1, (size_t) np);
285 for (int iq = 0; iq < ctl.nq; iq++) {
286 aqtdm[iq] = gsl_stats_mean(&aqtd[iq * NP], 1, (size_t) np);
287 rqtdm[iq] = gsl_stats_mean(&rqtd[iq * NP], 1, (size_t) np);
288 }
289 } else if (strcasecmp(argv[3], "stddev") == 0) {
290 ahtdm = gsl_stats_sd(ahtd, 1, (size_t) np);
291 rhtdm = gsl_stats_sd(rhtd, 1, (size_t) np);
292 avtdm = gsl_stats_sd(avtd, 1, (size_t) np);
293 rvtdm = gsl_stats_sd(rvtd, 1, (size_t) np);
294 for (int iq = 0; iq < ctl.nq; iq++) {
295 aqtdm[iq] = gsl_stats_sd(&aqtd[iq * NP], 1, (size_t) np);
296 rqtdm[iq] = gsl_stats_sd(&rqtd[iq * NP], 1, (size_t) np);
297 }
298 } else if (strcasecmp(argv[3], "min") == 0) {
299 ahtdm = gsl_stats_min(ahtd, 1, (size_t) np);
300 rhtdm = gsl_stats_min(rhtd, 1, (size_t) np);
301 avtdm = gsl_stats_min(avtd, 1, (size_t) np);
302 rvtdm = gsl_stats_min(rvtd, 1, (size_t) np);
303 for (int iq = 0; iq < ctl.nq; iq++) {
304 aqtdm[iq] = gsl_stats_min(&aqtd[iq * NP], 1, (size_t) np);
305 rqtdm[iq] = gsl_stats_min(&rqtd[iq * NP], 1, (size_t) np);
306 }
307 } else if (strcasecmp(argv[3], "max") == 0) {
308 ahtdm = gsl_stats_max(ahtd, 1, (size_t) np);
309 rhtdm = gsl_stats_max(rhtd, 1, (size_t) np);
310 avtdm = gsl_stats_max(avtd, 1, (size_t) np);
311 rvtdm = gsl_stats_max(rvtd, 1, (size_t) np);
312 for (int iq = 0; iq < ctl.nq; iq++) {
313 aqtdm[iq] = gsl_stats_max(&aqtd[iq * NP], 1, (size_t) np);
314 rqtdm[iq] = gsl_stats_max(&rqtd[iq * NP], 1, (size_t) np);
315 }
316 } else if (strcasecmp(argv[3], "skew") == 0) {
317 ahtdm = gsl_stats_skew(ahtd, 1, (size_t) np);
318 rhtdm = gsl_stats_skew(rhtd, 1, (size_t) np);
319 avtdm = gsl_stats_skew(avtd, 1, (size_t) np);
320 rvtdm = gsl_stats_skew(rvtd, 1, (size_t) np);
321 for (int iq = 0; iq < ctl.nq; iq++) {
322 aqtdm[iq] = gsl_stats_skew(&aqtd[iq * NP], 1, (size_t) np);
323 rqtdm[iq] = gsl_stats_skew(&rqtd[iq * NP], 1, (size_t) np);
324 }
325 } else if (strcasecmp(argv[3], "kurt") == 0) {
326 ahtdm = gsl_stats_kurtosis(ahtd, 1, (size_t) np);
327 rhtdm = gsl_stats_kurtosis(rhtd, 1, (size_t) np);
328 avtdm = gsl_stats_kurtosis(avtd, 1, (size_t) np);
329 rvtdm = gsl_stats_kurtosis(rvtd, 1, (size_t) np);
330 for (int iq = 0; iq < ctl.nq; iq++) {
331 aqtdm[iq] = gsl_stats_kurtosis(&aqtd[iq * NP], 1, (size_t) np);
332 rqtdm[iq] = gsl_stats_kurtosis(&rqtd[iq * NP], 1, (size_t) np);
333 }
334 } else if (strcasecmp(argv[3], "absdev") == 0) {
335 ahtdm = gsl_stats_absdev_m(ahtd, 1, (size_t) np, 0.0);
336 rhtdm = gsl_stats_absdev_m(rhtd, 1, (size_t) np, 0.0);
337 avtdm = gsl_stats_absdev_m(avtd, 1, (size_t) np, 0.0);
338 rvtdm = gsl_stats_absdev_m(rvtd, 1, (size_t) np, 0.0);
339 for (int iq = 0; iq < ctl.nq; iq++) {
340 aqtdm[iq] = gsl_stats_absdev_m(&aqtd[iq * NP], 1, (size_t) np, 0.0);
341 rqtdm[iq] = gsl_stats_absdev_m(&rqtd[iq * NP], 1, (size_t) np, 0.0);
342 }
343 } else if (strcasecmp(argv[3], "median") == 0) {
344 ahtdm = gsl_stats_median(ahtd, 1, (size_t) np);
345 rhtdm = gsl_stats_median(rhtd, 1, (size_t) np);
346 avtdm = gsl_stats_median(avtd, 1, (size_t) np);
347 rvtdm = gsl_stats_median(rvtd, 1, (size_t) np);
348 for (int iq = 0; iq < ctl.nq; iq++) {
349 aqtdm[iq] = gsl_stats_median(&aqtd[iq * NP], 1, (size_t) np);
350 rqtdm[iq] = gsl_stats_median(&rqtd[iq * NP], 1, (size_t) np);
351 }
352 } else if (strcasecmp(argv[3], "mad") == 0) {
353 ahtdm = gsl_stats_mad0(ahtd, 1, (size_t) np, work);
354 rhtdm = gsl_stats_mad0(rhtd, 1, (size_t) np, work);
355 avtdm = gsl_stats_mad0(avtd, 1, (size_t) np, work);
356 rvtdm = gsl_stats_mad0(rvtd, 1, (size_t) np, work);
357 for (int iq = 0; iq < ctl.nq; iq++) {
358 aqtdm[iq] = gsl_stats_mad0(&aqtd[iq * NP], 1, (size_t) np, work);
359 rqtdm[iq] = gsl_stats_mad0(&rqtd[iq * NP], 1, (size_t) np, work);
360 }
361 } else
362 ERRMSG("Unknown parameter!");
363
364 /* Write output... */
365 fprintf(out, "%.2f %.2f %g %g %g %g", t, t - t0,
366 ahtdm, rhtdm, avtdm, rvtdm);
367 for (int iq = 0; iq < ctl.nq; iq++) {
368 fprintf(out, " ");
369 fprintf(out, ctl.qnt_format[iq], aqtdm[iq]);
370 fprintf(out, " ");
371 fprintf(out, ctl.qnt_format[iq], rqtdm[iq]);
372 }
373 fprintf(out, " %d\n", np);
374 }
375
376 /* Close file... */
377 fclose(out);
378
379 /* Free... */
380 free(atm1);
381 free(atm2);
382 free(lon1_old);
383 free(lat1_old);
384 free(z1_old);
385 free(lh1);
386 free(lv1);
387 free(lon2_old);
388 free(lat2_old);
389 free(z2_old);
390 free(lh2);
391 free(lv2);
392 free(ahtd);
393 free(avtd);
394 free(aqtd);
395 free(rhtd);
396 free(rvtd);
397 free(rqtd);
398 free(work);
399
400 return EXIT_SUCCESS;
401}
402
403/*****************************************************************************/
404
406void usage(
407 void) {
408
409 printf("\nMPTRAC atm_dist tool.\n\n");
410 printf
411 ("Calculate transport deviations between pairs of trajectory data sets.\n");
412 printf("\n");
413 printf("Usage:\n");
414 printf
415 (" atm_dist <ctl> <dist.tab> <param> <atm1a> <atm1b> [<atm2a> <atm2b> ...]\n");
416 printf("\n");
417 printf("Arguments:\n");
418 printf(" <ctl> Control file.\n");
419 printf(" <dist.tab> Output table.\n");
420 printf
421 (" <param> Statistic: mean, stddev, min, max, skew, kurt, absdev,\n");
422 printf(" median, or mad.\n");
423 printf(" <atm*a/b> Atmospheric input files to compare pairwise.\n");
424 printf("\nFurther information:\n");
425 printf(" Manual: https://slcs-jsc.github.io/mptrac/\n");
426}
main
int main(int argc, char *argv[])
Definition: atm_dist.c:39
usage
void usage(void)
Print command-line help.
Definition: atm_dist.c:406
scan_ctl
double scan_ctl(const char *filename, int argc, char *argv[], const char *varname, const int arridx, const char *defvalue, char *value)
Scans a control file or command-line arguments for a specified variable.
Definition: mptrac.c:10745
time_from_filename
double time_from_filename(const char *filename, const int offset)
Extracts and converts a timestamp from a filename to Julian seconds.
Definition: mptrac.c:11017
mptrac_read_atm
int mptrac_read_atm(const char *filename, const ctl_t *ctl, atm_t *atm)
Reads air parcel data from a specified file into the given atmospheric structure.
Definition: mptrac.c:5117
mptrac_read_ctl
void mptrac_read_ctl(const char *filename, int argc, char *argv[], ctl_t *ctl)
Reads control parameters from a configuration file and populates the given structure.
Definition: mptrac.c:5248
geo2cart
void geo2cart(const double z, const double lon, const double lat, double *x)
Converts geographic coordinates (longitude, latitude, altitude) to Cartesian coordinates.
Definition: mptrac.c:1674
mptrac.h
MPTRAC library declarations.
ERRMSG
#define ERRMSG(...)
Print an error message with contextual information and terminate the program.
Definition: mptrac.h:2102
USAGE
#define USAGE
Print usage information on -h or --help.
Definition: mptrac.h:1909
Z
#define Z(p)
Convert pressure to altitude.
Definition: mptrac.h:1939
P
#define P(z)
Compute pressure at given altitude.
Definition: mptrac.h:1480
NQ
#define NQ
Maximum number of quantities per data point.
Definition: mptrac.h:360
ALLOC
#define ALLOC(ptr, type, n)
Allocate memory for a pointer with error handling.
Definition: mptrac.h:453
NP
#define NP
Maximum number of atmospheric data points.
Definition: mptrac.h:355
LOG
#define LOG(level,...)
Print a log message with a specified logging level.
Definition: mptrac.h:2032
DIST
#define DIST(a, b)
Calculate the distance between two points in Cartesian coordinates.
Definition: mptrac.h:704
atm_t
Air parcel data.
Definition: mptrac.h:3241
atm_t::time
double time[NP]
Time [s].
Definition: mptrac.h:3247
atm_t::lat
double lat[NP]
Latitude [deg].
Definition: mptrac.h:3256
atm_t::lon
double lon[NP]
Longitude [deg].
Definition: mptrac.h:3253
atm_t::np
int np
Number of air parcels.
Definition: mptrac.h:3244
atm_t::q
double q[NQ][NP]
Quantity data (for various, user-defined attributes).
Definition: mptrac.h:3259
atm_t::p
double p[NP]
Pressure [hPa].
Definition: mptrac.h:3250
ctl_t
Control parameters.
Definition: mptrac.h:2190
ctl_t::qnt_format
char qnt_format[NQ][LEN]
Quantity output format.
Definition: mptrac.h:2209
ctl_t::atm_type
int atm_type
Type of atmospheric data files (0=ASCII, 1=binary, 2=netCDF, 3=CLaMS_traj, 4=CLaMS_pos).
Definition: mptrac.h:3007
ctl_t::qnt_unit
char qnt_unit[NQ][LEN]
Quantity units.
Definition: mptrac.h:2206
ctl_t::qnt_name
char qnt_name[NQ][LEN]
Quantity names.
Definition: mptrac.h:2200
ctl_t::qnt_ens
int qnt_ens
Quantity array index for ensemble IDs.
Definition: mptrac.h:2215
ctl_t::qnt_idx
int qnt_idx
Quantity array index for air parcel IDs.
Definition: mptrac.h:2212
ctl_t::nq
int nq
Number of quantities.
Definition: mptrac.h:2197
Generated by doxygen 1.9.4