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 "process_int.h"
11 #include "cdo_cdi_wrapper.h"
12
13 #define NLON 1440
14 #define NLAT 720
15 #define MAX_VARS 6
16
17 static void
init_amsr_day(int vlistID,int gridID,int zaxisID,int nvars)18 init_amsr_day(int vlistID, int gridID, int zaxisID, int nvars)
19 {
20 /*
21 Version-5 RSS AMSR-E or AMSR-J daily files
22
23 filename with path in form satname_yyyymmdd_v5.gz
24 where satname = name of satellite (amsre or amsr)
25 yyyy = year
26 mm = month
27 dd = day of month
28
29 1:time time of measurement in fractional hours GMT
30 2:sst sea surface temperature in deg Celcius
31 3:wind 10m surface wind in meters/second
32 4:vapor columnar water vapor in millimeters
33 5:cloud cloud liquid water in millimeters
34 6:rain rain rate in millimeters/hour
35 */
36 const char *name[] = { "hours", "sst", "wind", "vapor", "cloud", "rain" };
37 const char *units[] = { "h", "deg Celcius", "m/s", "mm", "mm", "mm/h" };
38 const double xscale[] = { 0.1, 0.15, 0.2, 0.3, 0.01, 0.1 };
39 const double xminval[] = { 0., -3., 0., 0., 0., 0. };
40
41 for (int i = 0; i < nvars; ++i)
42 {
43 const int varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_VARYING);
44 cdiDefKeyString(vlistID, varID, CDI_KEY_NAME, name[i]);
45 cdiDefKeyString(vlistID, varID, CDI_KEY_UNITS, units[i]);
46 vlistDefVarDatatype(vlistID, varID, CDI_DATATYPE_INT16);
47 vlistDefVarMissval(vlistID, varID, 254);
48 vlistDefVarScalefactor(vlistID, varID, xscale[i]);
49 vlistDefVarAddoffset(vlistID, varID, xminval[i]);
50 }
51 }
52
53 static void
init_amsr_averaged(int vlistID,int gridID,int zaxisID,int nvars)54 init_amsr_averaged(int vlistID, int gridID, int zaxisID, int nvars)
55 {
56 /*
57 Version-5 AMSR-E or AMSR-J time-averaged files including:
58 3-day (average of 3 days ending on file date)
59 weekly (average of 7 days ending on Saturday of file date)
60 monthly (average of all days in month)
61
62
63 filename
64 format of file names are:
65 3-day satname_yyyymmddv5_d3d.gz
66 weekly satname_yyyymmddv5.gz
67 monthly satname_yyyymmv5.gz
68
69 where satname =name of satellite (amsre or amsr)
70 yyyy =year
71 mm =month
72 dd =day of month
73
74 1:sst sea surface temperature in deg Celcius
75 2:wind 10m surface wind in meters/second
76 3:vapor columnar water vapor in millimeters
77 4:cloud cloud liquid water in millimeters
78 5:rain rain rate in millimeters/hour
79 */
80 const char *name[] = { "sst", "wind", "vapor", "cloud", "rain" };
81 const char *units[] = { "deg Celcius", "m/s", "mm", "mm", "mm/h" };
82 const double xscale[] = { 0.15, 0.2, 0.3, 0.01, 0.1 };
83 const double xminval[] = { -3., 0., 0., 0., 0. };
84
85 for (int i = 0; i < nvars; ++i)
86 {
87 /* varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_CONSTANT); */
88 const int varID = vlistDefVar(vlistID, gridID, zaxisID, TIME_VARYING);
89 cdiDefKeyString(vlistID, varID, CDI_KEY_NAME, name[i]);
90 cdiDefKeyString(vlistID, varID, CDI_KEY_UNITS, units[i]);
91 vlistDefVarDatatype(vlistID, varID, CDI_DATATYPE_INT16);
92 vlistDefVarMissval(vlistID, varID, 254);
93 vlistDefVarScalefactor(vlistID, varID, xscale[i]);
94 vlistDefVarAddoffset(vlistID, varID, xminval[i]);
95 }
96 }
97
98 static void
read_amsr(FILE * fp,int vlistID,int nvars,Varray2D<double> & data,size_t * nmiss)99 read_amsr(FILE *fp, int vlistID, int nvars, Varray2D<double> &data, size_t *nmiss)
100 {
101 int varID, i, gridsize;
102 std::vector<unsigned char> amsr_data;
103 double xminval, xscale, missval;
104 size_t size;
105
106 for (varID = 0; varID < nvars; ++varID)
107 {
108 gridsize = gridInqSize(vlistInqVarGrid(vlistID, varID));
109 amsr_data.resize(gridsize);
110 size = fread(amsr_data.data(), 1, gridsize, fp);
111 if ((int) size != gridsize) cdo_abort("Read error!");
112
113 missval = vlistInqVarMissval(vlistID, varID);
114 xminval = vlistInqVarAddoffset(vlistID, varID);
115 xscale = vlistInqVarScalefactor(vlistID, varID);
116
117 nmiss[varID] = 0;
118 for (i = 0; i < gridsize; ++i)
119 {
120 if (amsr_data[i] <= 250)
121 {
122 data[varID][i] = amsr_data[i] * xscale + xminval;
123 }
124 else
125 {
126 data[varID][i] = missval;
127 nmiss[varID]++;
128 }
129 }
130 }
131 }
132
133 static void
write_data(CdoStreamID streamID,int nvars,Varray2D<double> & data,size_t * nmiss)134 write_data(CdoStreamID streamID, int nvars, Varray2D<double> &data, size_t *nmiss)
135 {
136 for (int varID = 0; varID < nvars; ++varID)
137 {
138 cdo_def_record(streamID, varID, 0);
139 cdo_write_record(streamID, data[varID].data(), nmiss[varID]);
140 }
141 }
142
143 static int
getDate(const char * name)144 getDate(const char *name)
145 {
146 const char *pname = strchr(name, '_');
147 const int date = pname ? atoi(pname + 1) : 0;
148 return date;
149 }
150
151 void *
Importamsr(void * process)152 Importamsr(void *process)
153 {
154 int tsID;
155 int nvars;
156 int vtime = 0;
157 double xvals[NLON], yvals[NLAT];
158 Varray2D<double> data(MAX_VARS);
159 size_t nmiss[MAX_VARS];
160
161 cdo_initialize(process);
162
163 operator_check_argc(0);
164
165 auto fp = fopen(cdo_get_stream_name(0), "r");
166 if (fp == nullptr)
167 {
168 perror(cdo_get_stream_name(0));
169 exit(EXIT_FAILURE);
170 }
171
172 fseek(fp, 0L, SEEK_END);
173 const size_t fsize = (size_t) ftell(fp);
174 fseek(fp, 0L, SEEK_SET);
175
176 auto vdate = getDate(cdo_get_stream_name(0));
177 if (vdate <= 999999) vdate = vdate * 100 + 1;
178
179 const auto streamID = cdo_open_write(1);
180
181 /*
182 Longitude is 0.25*xdim-0.125 degrees east
183 Latitude is 0.25*ydim-90.125
184 */
185 const size_t gridsize = NLON * NLAT;
186 const auto gridID = gridCreate(GRID_LONLAT, gridsize);
187 gridDefXsize(gridID, NLON);
188 gridDefYsize(gridID, NLAT);
189 for (int i = 0; i < NLON; ++i) xvals[i] = 0.25 * (i + 1) - 0.125;
190 for (int i = 0; i < NLAT; ++i) yvals[i] = 0.25 * (i + 1) - 90.125;
191 gridDefXvals(gridID, xvals);
192 gridDefYvals(gridID, yvals);
193
194 const auto zaxisID = zaxisCreate(ZAXIS_SURFACE, 1);
195
196 const auto vlistID = vlistCreate();
197 const auto taxisID = cdo_taxis_create(TAXIS_ABSOLUTE);
198 vlistDefTaxis(vlistID, taxisID);
199
200 if (fsize == 12441600)
201 {
202 nvars = 6;
203 for (int i = 0; i < nvars; ++i) data[i].resize(gridsize);
204
205 init_amsr_day(vlistID, gridID, zaxisID, nvars);
206
207 cdo_def_vlist(streamID, vlistID);
208
209 vtime = 13000; /* 1:30:00 */
210 for (tsID = 0; tsID < 2; ++tsID)
211 {
212 taxisDefVdate(taxisID, vdate);
213 taxisDefVtime(taxisID, vtime);
214 vtime += 120000; /* 13:30:00 */
215 cdo_def_timestep(streamID, tsID);
216
217 read_amsr(fp, vlistID, nvars, data, nmiss);
218
219 write_data(streamID, nvars, data, nmiss);
220 }
221 }
222 else if (fsize == 5184000)
223 {
224 nvars = 5;
225 for (int i = 0; i < nvars; ++i) data[i].resize(gridsize);
226
227 init_amsr_averaged(vlistID, gridID, zaxisID, nvars);
228
229 /* vlistDefNtsteps(vlistID, 0);*/
230 cdo_def_vlist(streamID, vlistID);
231
232 taxisDefVdate(taxisID, vdate);
233 taxisDefVtime(taxisID, vtime);
234 tsID = 0;
235 cdo_def_timestep(streamID, tsID);
236
237 read_amsr(fp, vlistID, nvars, data, nmiss);
238
239 write_data(streamID, nvars, data, nmiss);
240 }
241 else
242 cdo_abort("Unexpected file size for AMSR data!");
243
244 process_def_var_num(vlistNvars(vlistID));
245
246 cdo_stream_close(streamID);
247
248 fclose(fp);
249
250 vlistDestroy(vlistID);
251 gridDestroy(gridID);
252 zaxisDestroy(zaxisID);
253 taxisDestroy(taxisID);
254
255 cdo_finish();
256
257 return nullptr;
258 }
259