Importobs.cc 6.73 KB
Newer Older
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
6
7
8
9
10
11
12
13
14
15
16
17
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

Ralf Mueller's avatar
Ralf Mueller committed
18
#include <cdi.h>
Oliver Heidmann's avatar
Oliver Heidmann committed
19

20
#include "cdo_options.h"
21
#include "process_int.h"
22
#include "array.h"
23
#include "readline.h"
24
#include "cdo_cdi_wrapper.h"
25
#include <mpim_grid.h>
26
#include "griddes.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
28


29
30
static void
init_vars(int vlistID, int gridID, int zaxisID, int nvars)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
{
32
  const int code[] = { 11, 17, 33, 34, 1, 2 /*, 3*/ };
Uwe Schulzweida's avatar
Uwe Schulzweida committed
33
  const char *name[] = { "temp", "depoint", "u", "v", "height", "pressure" /*, "station"*/ };
34
  const char *units[] = { "Celsius", "", "m/s", "m/s", "m", "hPa" /*, ""*/ };
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35

36
  for (int i = 0; i < nvars; ++i)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
37
    {
38
      const auto varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_VARYING);
39
      vlistDefVarParam(vlistID, varID, cdiEncodeParam(code[i], 255, 255));
40
41
      cdiDefKeyString(vlistID, varID, CDI_KEY_NAME, name[i]);
      cdiDefKeyString(vlistID, varID, CDI_KEY_UNITS, units[i]);
42
      vlistDefVarDatatype(vlistID, varID, CDI_DATATYPE_FLT32);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
44
45
    }
}

46
static void
47
init_data(int vlistID, int nvars, Varray2D<double> &data)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
48
{
49
  for (int varID = 0; varID < nvars; ++varID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
    {
51
52
      const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID, varID));
      const auto missval = vlistInqVarMissval(vlistID, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
53
      for (size_t i = 0; i < gridsize; ++i) data[varID][i] = missval;
54
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
55
56
}

57
static void
58
write_data(CdoStreamID streamID, int vlistID, int nvars, Varray2D<double> &data)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59
{
60
  for (int varID = 0; varID < nvars; ++varID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
61
    {
62
63
      const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID, varID));
      const auto missval = vlistInqVarMissval(vlistID, varID);
64
      const auto nmiss = varrayNumMV(gridsize, data[varID], missval);
65
      cdoDefRecord(streamID, varID, 0);
66
      cdoWriteRecord(streamID, data[varID].data(), nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67
68
69
    }
}

70
71
static int
getDate(const char *name)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
72
{
73
  const char *pname = strchr(name, '_');
74
  const int date = pname ? atoi(pname + 1) : 0;
75
  return date;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
76
77
}

78
#define MAX_LINE_LEN 4096
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79

80
81
void *
Importobs(void *process)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
83
{
  char line[MAX_LINE_LEN];
84
85
86
  size_t i, j;
  long index;
  constexpr int nvars = 6;
87
  int vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
88
89
  char dummy[32], station[32], datetime[32];
  float lat, lon, height1, pressure, height2, value;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90
  double latmin = 90, latmax = -90, lonmin = 360, lonmax = -360;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
92
  int code;

93
  cdoInitialize(process);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94

95
  cdoOperatorAdd("import_obs", 0, 0, "grid description file or name");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96

97
  const auto operatorID = cdoOperatorID();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
98

99
  operatorInputArg(cdoOperatorEnter(operatorID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100

101
  const auto gridID = cdoDefineGrid(cdoOperatorArgv(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102

Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
  if (gridInqType(gridID) != GRID_LONLAT) cdoAbort("Unsupported grid type: %s", gridNamePtr(gridInqType(gridID)));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104

105
  const auto gridsize = gridInqSize(gridID);
106
107
  const auto xsize = gridInqXsize(gridID);
  const auto ysize = gridInqYsize(gridID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108

Uwe Schulzweida's avatar
Uwe Schulzweida committed
109
  Varray<double> xvals(gridsize), yvals(gridsize);
110
111
  gridInqXvals(gridID, xvals.data());
  gridInqYvals(gridID, yvals.data());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112

113
  // Convert lat/lon units if required
114
115
  cdo_grid_to_degree(gridID, CDI_XAXIS, gridsize, xvals.data(), "grid center lon");
  cdo_grid_to_degree(gridID, CDI_YAXIS, gridsize, yvals.data(), "grid center lat");
116

Uwe Schulzweida's avatar
Uwe Schulzweida committed
117
  auto fp = fopen(cdoGetStreamName(0), "r");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118
  if (fp == nullptr)
119
    {
120
      perror(cdoGetStreamName(0));
121
122
      exit(EXIT_FAILURE);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
123

124
  auto vdate = getDate(cdoGetStreamName(0));
125
  if (vdate <= 999999) vdate = vdate * 100 + 1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
126

127
  const auto streamID = cdoOpenWrite(1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128

129
  const auto zaxisID = zaxisCreate(ZAXIS_SURFACE, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130

131
  const auto vlistID = vlistCreate();
132

133
  const auto taxisID = cdoTaxisCreate(TAXIS_ABSOLUTE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
134
135
  vlistDefTaxis(vlistID, taxisID);

136
137
  Varray2D<double> data(nvars);
  for (i = 0; i < nvars; ++i) data[i].resize(gridsize);
138
139
140

  init_vars(vlistID, gridID, zaxisID, nvars);

141
  cdoDefVlist(streamID, vlistID);
142

143
  int64_t vdate0 = 0;
144
  int vtime0 = 0;
145
  // ntime = 0;
146
  int tsID = 0;
147
  while (cdo::readline(fp, line, MAX_LINE_LEN))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
148
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
149
150
      sscanf(line, "%s %s %s %g %g %g %d %g %g %g", dummy, station, datetime, &lat, &lon, &height1, &code, &pressure, &height2,
             &value);
151
152
153
      long vdate_l;
      sscanf(datetime, "%ld_%d", &vdate_l, &vtime);
      vdate = vdate_l;
154

155
      if (vdate != vdate0 || vtime != vtime0)
156
        {
157
          if (tsID > 0) write_data(streamID, vlistID, nvars, data);
158
159
160

          vdate0 = vdate;
          vtime0 = vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
          // printf("%s %d %d %g %g %g %d %g %g %g\n", station, vdate, vtime, lat, lon, height1, code, pressure, height2, value);
162
163
          taxisDefVdate(taxisID, vdate);
          taxisDefVtime(taxisID, vtime);
164
          cdoDefTimestep(streamID, tsID);
165

166
          init_data(vlistID, nvars, data);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167

168
169
          tsID++;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
170

171
172
173
174
175
      if (lon < lonmin) lonmin = lon;
      if (lon > lonmax) lonmax = lon;
      if (lat < latmin) latmin = lat;
      if (lat > latmax) latmax = lat;

176
      const auto dy = yvals[1] - yvals[0];
177
178
      for (j = 0; j < ysize; ++j)
        if (lat >= (yvals[j] - dy / 2) && lat < (yvals[j] + dy / 2)) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
179

180
      const auto dx = xvals[1] - xvals[0];
181
182
183
      if (lon < (xvals[0] - dx / 2) && lon < 0) lon += 360;
      for (i = 0; i < xsize; ++i)
        if (lon >= (xvals[i] - dx / 2) && lon < (xvals[i] + dx / 2)) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184

185
      index = -1;
186
187
188
189
      if (code == 11) index = 0;
      if (code == 17) index = 1;
      if (code == 33) index = 2;
      if (code == 34) index = 3;
190

191
192
      // printf("%d %d %d %g %g %g %g\n", i, j, index, dx, dy, lon, lat);
      if (i < xsize && j < ysize && index >= 0)
193
        {
194
          char *pstation = station;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
195
          while (isalpha(*pstation)) pstation++;
196
          // printf("station %s %d\n", pstation, atoi(pstation));
197
198
199
          data[index][j * xsize + i] = value;
          data[4][j * xsize + i] = height1;
          data[5][j * xsize + i] = pressure;
200
201
202
203
          // data[    6][j*xsize+i] = atoi(pstation);
        }

      /*
204
        printf("%s %d %d %g %g %g %d %g %g %g\n", station, vdate, vtime, lat, lon, height1, code, pressure, height2, value);
205
      */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206
    }
207
208
209

  write_data(streamID, vlistID, nvars, data);

210
  if (Options::cdoVerbose) printf("lonmin=%g, lonmax=%g, latmin=%g, latmax=%g\n", lonmin, lonmax, latmin, latmax);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211

212
  processDefVarNum(vlistNvars(vlistID));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
213

214
  cdoStreamClose(streamID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
215
216
217
218
219
220
221
222
223
224

  fclose(fp);

  vlistDestroy(vlistID);
  gridDestroy(gridID);
  zaxisDestroy(zaxisID);
  taxisDestroy(taxisID);

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
225
  return nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
226
}