1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Uwe Schulzweida
5
6 */
7
8 #include <cdi.h>
9
10 #include "cdo_options.h"
11 #include "process_int.h"
12 #include "varray.h"
13 #include "readline.h"
14 #include "cdo_cdi_wrapper.h"
15 #include <mpim_grid.h>
16 #include "griddes.h"
17
18 static void
init_vars(int vlistID,int gridID,int zaxisID,int nvars)19 init_vars(int vlistID, int gridID, int zaxisID, int nvars)
20 {
21 const int code[] = { 11, 17, 33, 34, 1, 2 /*, 3*/ };
22 const char *name[] = { "temp", "depoint", "u", "v", "height", "pressure" /*, "station"*/ };
23 const char *units[] = { "Celsius", "", "m/s", "m/s", "m", "hPa" /*, ""*/ };
24
25 for (int i = 0; i < nvars; ++i)
26 {
27 const auto varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_VARYING);
28 vlistDefVarParam(vlistID, varID, cdiEncodeParam(code[i], 255, 255));
29 cdiDefKeyString(vlistID, varID, CDI_KEY_NAME, name[i]);
30 cdiDefKeyString(vlistID, varID, CDI_KEY_UNITS, units[i]);
31 vlistDefVarDatatype(vlistID, varID, CDI_DATATYPE_FLT32);
32 }
33 }
34
35 static void
init_data(int vlistID,int nvars,Varray2D<double> & data)36 init_data(int vlistID, int nvars, Varray2D<double> &data)
37 {
38 for (int varID = 0; varID < nvars; ++varID)
39 {
40 const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID, varID));
41 const auto missval = vlistInqVarMissval(vlistID, varID);
42 for (size_t i = 0; i < gridsize; ++i) data[varID][i] = missval;
43 }
44 }
45
46 static void
write_data(CdoStreamID streamID,int vlistID,int nvars,Varray2D<double> & data)47 write_data(CdoStreamID streamID, int vlistID, int nvars, Varray2D<double> &data)
48 {
49 for (int varID = 0; varID < nvars; ++varID)
50 {
51 const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID, varID));
52 const auto missval = vlistInqVarMissval(vlistID, varID);
53 const auto nmiss = varray_num_mv(gridsize, data[varID], missval);
54 cdo_def_record(streamID, varID, 0);
55 cdo_write_record(streamID, data[varID].data(), nmiss);
56 }
57 }
58
59 static int
getDate(const char * name)60 getDate(const char *name)
61 {
62 const char *pname = strchr(name, '_');
63 const int date = pname ? atoi(pname + 1) : 0;
64 return date;
65 }
66
67 #define MAX_LINE_LEN 4096
68
69 void *
Importobs(void * process)70 Importobs(void *process)
71 {
72 char line[MAX_LINE_LEN];
73 size_t i, j;
74 long index;
75 constexpr int nvars = 6;
76 int vtime = 0;
77 char dummy[32], station[32], datetime[32];
78 float lat, lon, height1, pressure, height2, value;
79 double latmin = 90, latmax = -90, lonmin = 360, lonmax = -360;
80 int code;
81
82 cdo_initialize(process);
83
84 cdo_operator_add("import_obs", 0, 0, "grid description file or name");
85
86 const auto operatorID = cdo_operator_id();
87
88 operator_input_arg(cdo_operator_enter(operatorID));
89
90 const auto gridID = cdo_define_grid(cdo_operator_argv(0));
91
92 if (gridInqType(gridID) != GRID_LONLAT) cdo_abort("Unsupported grid type: %s", gridNamePtr(gridInqType(gridID)));
93
94 const auto gridsize = gridInqSize(gridID);
95 const auto xsize = gridInqXsize(gridID);
96 const auto ysize = gridInqYsize(gridID);
97
98 Varray<double> xvals(gridsize), yvals(gridsize);
99 gridInqXvals(gridID, xvals.data());
100 gridInqYvals(gridID, yvals.data());
101
102 // Convert lat/lon units if required
103 cdo_grid_to_degree(gridID, CDI_XAXIS, gridsize, xvals.data(), "grid center lon");
104 cdo_grid_to_degree(gridID, CDI_YAXIS, gridsize, yvals.data(), "grid center lat");
105
106 auto fp = fopen(cdo_get_stream_name(0), "r");
107 if (fp == nullptr)
108 {
109 perror(cdo_get_stream_name(0));
110 exit(EXIT_FAILURE);
111 }
112
113 auto vdate = getDate(cdo_get_stream_name(0));
114 if (vdate <= 999999) vdate = vdate * 100 + 1;
115
116 const auto streamID = cdo_open_write(1);
117
118 const auto zaxisID = zaxisCreate(ZAXIS_SURFACE, 1);
119
120 const auto vlistID = vlistCreate();
121
122 const auto taxisID = cdo_taxis_create(TAXIS_ABSOLUTE);
123 vlistDefTaxis(vlistID, taxisID);
124
125 Varray2D<double> data(nvars);
126 for (i = 0; i < nvars; ++i) data[i].resize(gridsize);
127
128 init_vars(vlistID, gridID, zaxisID, nvars);
129
130 cdo_def_vlist(streamID, vlistID);
131
132 int64_t vdate0 = 0;
133 int vtime0 = 0;
134 // ntime = 0;
135 int tsID = 0;
136 while (cdo::readline(fp, line, MAX_LINE_LEN))
137 {
138 sscanf(line, "%s %s %s %g %g %g %d %g %g %g", dummy, station, datetime, &lat, &lon, &height1, &code, &pressure, &height2,
139 &value);
140 long vdate_l;
141 sscanf(datetime, "%ld_%d", &vdate_l, &vtime);
142 vdate = vdate_l;
143
144 if (vdate != vdate0 || vtime != vtime0)
145 {
146 if (tsID > 0) write_data(streamID, vlistID, nvars, data);
147
148 vdate0 = vdate;
149 vtime0 = vtime;
150 // printf("%s %d %d %g %g %g %d %g %g %g\n", station, vdate, vtime, lat, lon, height1, code, pressure, height2, value);
151 taxisDefVdate(taxisID, vdate);
152 taxisDefVtime(taxisID, vtime);
153 cdo_def_timestep(streamID, tsID);
154
155 init_data(vlistID, nvars, data);
156
157 tsID++;
158 }
159
160 if (lon < lonmin) lonmin = lon;
161 if (lon > lonmax) lonmax = lon;
162 if (lat < latmin) latmin = lat;
163 if (lat > latmax) latmax = lat;
164
165 const auto dy = yvals[1] - yvals[0];
166 for (j = 0; j < ysize; ++j)
167 if (lat >= (yvals[j] - dy / 2) && lat < (yvals[j] + dy / 2)) break;
168
169 const auto dx = xvals[1] - xvals[0];
170 if (lon < (xvals[0] - dx / 2) && lon < 0) lon += 360;
171 for (i = 0; i < xsize; ++i)
172 if (lon >= (xvals[i] - dx / 2) && lon < (xvals[i] + dx / 2)) break;
173
174 index = -1;
175 if (code == 11) index = 0;
176 if (code == 17) index = 1;
177 if (code == 33) index = 2;
178 if (code == 34) index = 3;
179
180 // printf("%d %d %d %g %g %g %g\n", i, j, index, dx, dy, lon, lat);
181 if (i < xsize && j < ysize && index >= 0)
182 {
183 char *pstation = station;
184 while (isalpha(*pstation)) pstation++;
185 // printf("station %s %d\n", pstation, atoi(pstation));
186 data[index][j * xsize + i] = value;
187 data[4][j * xsize + i] = height1;
188 data[5][j * xsize + i] = pressure;
189 // data[ 6][j*xsize+i] = atoi(pstation);
190 }
191
192 /*
193 printf("%s %d %d %g %g %g %d %g %g %g\n", station, vdate, vtime, lat, lon, height1, code, pressure, height2, value);
194 */
195 }
196
197 write_data(streamID, vlistID, nvars, data);
198
199 if (Options::cdoVerbose) printf("lonmin=%g, lonmax=%g, latmin=%g, latmax=%g\n", lonmin, lonmax, latmin, latmax);
200
201 process_def_var_num(vlistNvars(vlistID));
202
203 cdo_stream_close(streamID);
204
205 fclose(fp);
206
207 vlistDestroy(vlistID);
208 gridDestroy(gridID);
209 zaxisDestroy(zaxisID);
210 taxisDestroy(taxisID);
211
212 cdo_finish();
213
214 return nullptr;
215 }
216