Skip to content
GitLab
Commit 6ade6d78 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

added support of timerange for GRIB2

parent 23d7687b
Hide whitespace changes
Inline Side-by-side
src/stream_cgribex.c
View file @ 6ade6d78
......@@ -926,15 +926,14 @@ int cgribexScanTimestep2(int streamID)
if ( taxisInqType(taxisID) == TAXIS_RELATIVE )
{
taxis->type = TAXIS_RELATIVE;
taxis->unit = cgribexGetTimeUnit(isec1);
taxis->rdate = gribRefDate(isec1);
taxis->rtime = gribRefTime(isec1);
}
else
{
taxis->type = TAXIS_ABSOLUTE;
taxis->unit = cgribexGetTimeUnit(isec1);
}
taxis->unit = cgribexGetTimeUnit(isec1);
taxis->vdate = vdate;
taxis->vtime = vtime;
......@@ -1209,15 +1208,14 @@ int cgribexScanTimestep(int streamID)
if ( taxisInqType(taxisID) == TAXIS_RELATIVE )
{
taxis->type = TAXIS_RELATIVE;
taxis->unit = cgribexGetTimeUnit(isec1);
taxis->rdate = gribRefDate(isec1);
taxis->rtime = gribRefTime(isec1);
}
else
{
taxis->type = TAXIS_ABSOLUTE;
taxis->unit = cgribexGetTimeUnit(isec1);
}
taxis->unit = cgribexGetTimeUnit(isec1);
taxis->vdate = vdate;
taxis->vtime = vtime;
......@@ -1514,101 +1512,131 @@ void cgribexDefParam(int *isec1, int param)
}
static
void cgribexDefTime(int *isec1, int date, int time, int tsteptype, int numavg, int taxisID)
int cgribexDefTimerange(int tsteptype, int factor, int calendar,
int rdate, int rtime, int vdate, int vtime, int *pip1, int *pip2)
{
int timerange = -1;
int year, month, day, hour, minute, second;
int julday1, secofday1, julday2, secofday2, days, secs;
int ip, ip1 = 0, ip2 = 0;
cdiDecodeDate(rdate, &year, &month, &day);
cdiDecodeTime(rtime, &hour, &minute, &second);
encode_juldaysec(calendar, year, month, day, hour, minute, &julday1, &secofday1);
cdiDecodeDate(vdate, &year, &month, &day);
cdiDecodeTime(vtime, &hour, &minute, &second);
encode_juldaysec(calendar, year, month, day, hour, minute, &julday2, &secofday2);
(void) julday_sub(julday1, secofday1, julday2, secofday2, &days, &secs);
if ( !(int) fmod(days*86400.0 + secs, factor) )
{
ip = (int) ((days*86400.0 + secs)/factor);
switch ( tsteptype )
{
case TSTEP_INSTANT: timerange = 0; ip1 = ip; ip2 = 0; break;
case TSTEP_INSTANT2: timerange = 1; ip1 = 0; ip2 = 0; break;
case TSTEP_RANGE: timerange = 2; ip1 = 0; ip2 = ip; break;
case TSTEP_AVG: timerange = 3; ip1 = 0; ip2 = ip; break;
case TSTEP_ACCUM: timerange = 4; ip1 = 0; ip2 = ip; break;
case TSTEP_DIFF: timerange = 5; ip1 = 0; ip2 = ip; break;
case TSTEP_INSTANT3:
default: timerange = 10; ip1 = ip; ip2 = 0; break;
}
}
*pip1 = ip1;
*pip2 = ip2;
return (timerange);
}
static
int cgribexDefDateTime(int *isec1, int timeunit, int date, int time)
{
int year, month, day, hour, minute, second;
int century = 0;
int factor = 1;
cdiDecodeDate(date, &year, &month, &day);
cdiDecodeTime(time, &hour, &minute, &second);
century = year / 100;
ISEC1_Year = year - century*100;
if ( year < 0 )
{
century = -century;
ISEC1_Year = -ISEC1_Year;
}
if ( ISEC1_Year == 0 )
{
ISEC1_Year = 100;
century -= 1;
}
century += 1;
if ( year < 0 ) century = -century;
ISEC1_Month = month;
ISEC1_Day = day;
ISEC1_Hour = hour;
ISEC1_Minute = minute;
ISEC1_Century = century;
switch (timeunit)
{
case TUNIT_MINUTE: factor = 60; ISEC1_TimeUnit = ISEC1_TABLE4_MINUTE; break;
case TUNIT_QUARTER: factor = 900; ISEC1_TimeUnit = ISEC1_TABLE4_QUARTER; break;
case TUNIT_HOUR: factor = 3600; ISEC1_TimeUnit = ISEC1_TABLE4_HOUR; break;
case TUNIT_DAY: factor = 86400; ISEC1_TimeUnit = ISEC1_TABLE4_DAY; break;
default: factor = 3600; ISEC1_TimeUnit = ISEC1_TABLE4_HOUR; break;
}
return (factor);
}
static
void cgribexDefTime(int *isec1, int vdate, int vtime, int tsteptype, int numavg, int taxisID)
{
int timetype = -1;
int timerange = 0;
int timeunit;
if ( taxisID != -1 ) timetype = taxisInqType(taxisID);
timeunit = taxisInqTunit(taxisID);
if ( timetype == TAXIS_RELATIVE )
{
int factor = 1;
int rdate, rtime;
int julday1, secofday1, julday2, secofday2, days, secs;
int ip, ip1 = 0, ip2 = 0;
int ip1 = 0, ip2 = 0;
int calendar;
calendar = taxisInqCalendar(taxisID);
rdate = taxisInqRdate(taxisID);
rtime = taxisInqRtime(taxisID);
cdiDecodeDate(rdate, &year, &month, &day);
cdiDecodeTime(rtime, &hour, &minute, &second);
century = year / 100;
ISEC1_Year = year - century*100;
if ( year < 0 )
{
century = -century;
ISEC1_Year = -ISEC1_Year;
}
if ( ISEC1_Year == 0 )
{
ISEC1_Year = 100;
century -= 1;
}
century += 1;
if ( year < 0 ) century = -century;
ISEC1_Month = month;
ISEC1_Day = day;
ISEC1_Hour = hour;
ISEC1_Minute = minute;
/*
printf("year %d, month %d, day %d, hour %d, minute %d\n", year, month, day, hour, minute);
*/
encode_juldaysec(calendar, year, month, day, hour, minute, &julday1, &secofday1);
switch (taxisInqTunit(taxisID))
{
case TUNIT_MINUTE: factor = 60; ISEC1_TimeUnit = ISEC1_TABLE4_MINUTE; break;
case TUNIT_QUARTER: factor = 900; ISEC1_TimeUnit = ISEC1_TABLE4_QUARTER; break;
case TUNIT_HOUR: factor = 3600; ISEC1_TimeUnit = ISEC1_TABLE4_HOUR; break;
case TUNIT_DAY: factor = 86400; ISEC1_TimeUnit = ISEC1_TABLE4_DAY; break;
default: factor = 3600; ISEC1_TimeUnit = ISEC1_TABLE4_HOUR; break;
}
cdiDecodeDate(date, &year, &month, &day);
cdiDecodeTime(time, &hour, &minute, &second);
/*
printf("date time %d %d\n", date, time);
printf("year %d, month %d, day %d, hour %d, minute %d\n", year, month, day, hour, minute);
*/
encode_juldaysec(calendar, year, month, day, hour, minute, &julday2, &secofday2);
(void) julday_sub(julday1, secofday1, julday2, secofday2, &days, &secs);
/* ip = (int) ((days*86400.0 + secs)/factor); */
if ( (int) fmod(days*86400.0 + secs, factor) )
ip = -1;
else
ip = (int) ((days*86400.0 + secs)/factor);
factor = cgribexDefDateTime(isec1, timeunit, rdate, rtime);
switch ( tsteptype )
{
case TSTEP_INSTANT: timerange = 0; ip1 = ip; ip2 = 0; break;
case TSTEP_INSTANT2: timerange = 1; ip1 = 0; ip2 = 0; break;
case TSTEP_RANGE: timerange = 2; ip1 = 0; ip2 = ip; break;
case TSTEP_AVG: timerange = 3; ip1 = 0; ip2 = ip; break;
case TSTEP_ACCUM: timerange = 4; ip1 = 0; ip2 = ip; break;
case TSTEP_DIFF: timerange = 5; ip1 = 0; ip2 = ip; break;
case TSTEP_INSTANT3:
default: timerange = 10; ip1 = ip; ip2 = 0; break;
}
// printf("timerange: %d %d %d\n", timerange, ip1, ip2);
timerange = cgribexDefTimerange(tsteptype, factor, calendar,
rdate, rtime, vdate, vtime, &ip1, &ip2);
if ( timerange == 10 )
{
if ( ip < 0 || ip > 0xffff ) timetype = TAXIS_ABSOLUTE;
if ( ip1 < 0 || ip1 > 0xffff ) timetype = TAXIS_ABSOLUTE;
if ( ip2 < 0 || ip2 > 0xffff ) timetype = TAXIS_ABSOLUTE;
}
else
{
if ( ip < 0 || ip > 0xff ) timetype = TAXIS_ABSOLUTE;
if ( ip1 < 0 || ip1 > 0xff ) timetype = TAXIS_ABSOLUTE;
if ( ip2 < 0 || ip2 > 0xff ) timetype = TAXIS_ABSOLUTE;
}
ISEC1_TimeRange = timerange;
......@@ -1618,42 +1646,7 @@ void cgribexDefTime(int *isec1, int date, int time, int tsteptype, int numavg, i
if ( timetype == TAXIS_ABSOLUTE )
{
cdiDecodeDate(date, &year, &month, &day);
cdiDecodeTime(time, &hour, &minute, &second);
century = year / 100;
ISEC1_Year = year - century*100;
if ( year < 0 )
{
century = -century;
ISEC1_Year = -ISEC1_Year;
}
if ( ISEC1_Year == 0 )
{
ISEC1_Year = 100;
century -= 1;
}
century += 1;
if ( year < 0 ) century = -century;
ISEC1_Month = month;
ISEC1_Day = day;
ISEC1_Hour = hour;
ISEC1_Minute = minute;
/* ISEC1_TimeUnit = 0; */
switch (taxisInqTunit(taxisID))
{
case TUNIT_MINUTE: ISEC1_TimeUnit = ISEC1_TABLE4_MINUTE; break;
case TUNIT_QUARTER: ISEC1_TimeUnit = ISEC1_TABLE4_QUARTER; break;
case TUNIT_HOUR: ISEC1_TimeUnit = ISEC1_TABLE4_HOUR; break;
case TUNIT_DAY: ISEC1_TimeUnit = ISEC1_TABLE4_DAY; break;
default: ISEC1_TimeUnit = ISEC1_TABLE4_MINUTE; break;
}
(void) cgribexDefDateTime(isec1, timeunit, vdate, vtime);
/*
if ( numavg > 0 )
......@@ -1661,14 +1654,12 @@ void cgribexDefTime(int *isec1, int date, int time, int tsteptype, int numavg, i
else
*/
ISEC1_TimeRange = 10;
ISEC1_TimePeriod1 = 0;
ISEC1_TimePeriod2 = 0;
}
ISEC1_AvgNum = numavg;
ISEC1_AvgMiss = 0;
ISEC1_Century = century;
ISEC1_DecScaleFactor = 0;
}
......@@ -2094,7 +2085,7 @@ void cgribexDefaultSec4(int *isec4)
size_t cgribexEncode(int varID, int levelID, int vlistID, int gridID, int zaxisID,
int date, int time, int tsteptype, int numavg,
int vdate, int vtime, int tsteptype, int numavg,
long datasize, const double *data, int nmiss, unsigned char *gribbuffer, size_t gribbuffersize)
{
static char func[] = "cgribexEncode";
......@@ -2123,7 +2114,7 @@ size_t cgribexEncode(int varID, int levelID, int vlistID, int gridID, int zaxisI
datatype = vlistInqVarDatatype(vlistID, varID);
cgribexDefParam(isec1, param);
cgribexDefTime(isec1, date, time, tsteptype, numavg, vlistInqTaxis(vlistID));
cgribexDefTime(isec1, vdate, vtime, tsteptype, numavg, vlistInqTaxis(vlistID));
cgribexDefGrid(isec1, isec2, gridID);
cgribexDefLevel(isec1, isec2, fsec2, zaxisID, levelID);
cgribexDefMask(isec3);
......
src/stream_cgribex.h
View file @ 6ade6d78
......@@ -9,7 +9,7 @@ int cgribexDecode(unsigned char *gribbuffer, int gribsize, double *data, int gri
int unreduced, int *nmiss, int *zip, double missval);
size_t cgribexEncode(int varID, int levelID, int vlistID, int gridID, int zaxisID,
int date, int time, int tsteptype, int numavg,
int vdate, int vtime, int tsteptype, int numavg,
long datasize, const double *data, int nmiss, unsigned char *gribbuffer, size_t gribbuffersize);
#endif /* _STREAM_CGRIBEX_H */
src/stream_grb.c
View file @ 6ade6d78
......@@ -390,12 +390,14 @@ size_t grbEncode(int filetype, int varID, int levelID, int vlistID, int gridID,
if ( filetype == FILETYPE_GRB )
{
nbytes = cgribexEncode(varID, levelID, vlistID, gridID, zaxisID, date, time, tsteptype, numavg,
nbytes = cgribexEncode(varID, levelID, vlistID, gridID, zaxisID,
date, time, tsteptype, numavg,
datasize, data, nmiss, gribbuffer, gribbuffersize);
}
else
{
nbytes = gribapiEncode(varID, levelID, vlistID, gridID, zaxisID, date, time, numavg,
nbytes = gribapiEncode(varID, levelID, vlistID, gridID, zaxisID,
date, time, tsteptype, numavg,
datasize, data, nmiss, gribbuffer, gribbuffersize, ljpeg, gribHandle);
}
......
src/stream_gribapi.c
View file @ 6ade6d78
......@@ -174,6 +174,114 @@ int gribapiGetZaxisType(int grb_ltype)
}
#endif
#if defined (HAVE_LIBGRIB_API)
static
int gribapiGetTimeUnits(grib_handle *gh)
{
static char func[] = "gribapiGetTimeUnits";
int timeunits = -1;
int stepunits;
long lpar;
static int lprint = TRUE;
GRIB_CHECK(grib_get_long(gh, "stepUnits", &lpar), 0);
stepunits = (int) lpar;
switch ( stepunits )
{
case 0: timeunits = TUNIT_MINUTE; break;
case 1: timeunits = TUNIT_HOUR; break;
case 2: timeunits = TUNIT_DAY; break;
case 3: timeunits = TUNIT_MONTH; break;
case 4: timeunits = TUNIT_YEAR; break;
default:
if ( lprint )
{
Message(func, "Step units %d unsupported", stepunits);
lprint = FALSE;
}
}
return (timeunits);
}
#endif
#if defined (HAVE_LIBGRIB_API)
static
int gribapiTimeIsFC(grib_handle *gh)
{
int isFC = TRUE;
long sigofrtime;
GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
if ( sigofrtime == 3 ) isFC = FALSE;
return (isFC);
}
#endif
#if defined (HAVE_LIBGRIB_API)
static
int gribapiGetTsteptype(grib_handle *gh)
{
static char func[] = "gribapiGetTsteptype";
int tsteptype = 0;
int timerange;
long lpar;
static int lprint = TRUE;
if ( gribapiTimeIsFC(gh) )
{
GRIB_CHECK(grib_get_long(gh, "stepType", &lpar), 0);
timerange = (int) lpar;
// printf("timerange %d\n", timerange);
switch ( timerange )
{
case 0: tsteptype = TSTEP_AVG; break;
case 1: tsteptype = TSTEP_ACCUM; break;
case 2: tsteptype = TSTEP_MIN; break;
case 3: tsteptype = TSTEP_MAX; break;
case 4: tsteptype = TSTEP_DIFF; break;
default:
if ( lprint )
{
Message(func, "Time range %d unsupported", timerange);
lprint = FALSE;
}
}
}
return (tsteptype);
}
#endif
#if defined (HAVE_LIBGRIB_API)
void gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
{
long lpar;
long sigofrtime;
GRIB_CHECK(grib_get_long(gh, "significanceOfReferenceTime", &sigofrtime), 0);
if ( sigofrtime == 3 )
{
GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
*vdate = (int) lpar;
GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);
*vtime = (int) lpar*100;
}
else
{
GRIB_CHECK(grib_get_long(gh, "validityDate", &lpar), 0);
*vdate = (int) lpar;
GRIB_CHECK(grib_get_long(gh, "validityTime", &lpar), 0);
*vtime = (int) lpar*100;
}
}
#endif
#if defined (HAVE_LIBGRIB_API)
static
void gribapiAddRecord(int streamID, int param, grib_handle *gh,
......@@ -188,6 +296,7 @@ void gribapiAddRecord(int streamID, int param, grib_handle *gh,
int tsID, recID;
int level1, level2;
int numavg;
int tsteptype;
int lbounds = 0;
record_t *record;
grid_t grid;
......@@ -209,6 +318,7 @@ void gribapiAddRecord(int streamID, int param, grib_handle *gh,
recID = recordNewEntry(streamID, tsID);
record = &streamptr->tsteps[tsID].records[recID];
tsteptype = gribapiGetTsteptype(gh);
// numavg = ISEC1_AvgNum;
numavg = 0;
/*
......@@ -516,11 +626,14 @@ void gribapiAddRecord(int streamID, int param, grib_handle *gh,
{
int modelID;
long processID;
GRIB_CHECK(grib_get_long(gh, "generatingProcessIdentifier", &processID), 0);
modelID = modelInq(varInqInst(varID), processID, NULL);
if ( modelID == CDI_UNDEFID )
modelID = modelDef(varInqInst(varID), processID, NULL);
varDefModel(varID, modelID);
status = grib_get_long(gh, "generatingProcessIdentifier", &processID);
if ( status == 0 )
{
modelID = modelInq(varInqInst(varID), processID, NULL);
if ( modelID == CDI_UNDEFID )
modelID = modelDef(varInqInst(varID), processID, NULL);
varDefModel(varID, modelID);
}
}
if ( varInqTable(varID) == CDI_UNDEFID )
......@@ -562,7 +675,8 @@ void gribapiScanTimestep1(int streamID)
int status;
int fileID;
int rtabnum = 0;
int rcode = 0, level1 = 0, level2 = 0, vdate = 0, vtime = 0;
int rcode = 0, level1 = 0, level2 = 0;
int vdate = 0, vtime = 0;
int param = 0;
DateTime datetime, datetime0;
int tsID;
......@@ -701,10 +815,10 @@ void gribapiScanTimestep1(int streamID)
level1 = (int) dlevel;
level2 = 0;
GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
vdate = (int) lpar;
GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);
vtime = (int) lpar*100;
gribapiGetValidityDateTime(gh, &vdate, &vtime);
/*
printf("%d %d %d.%d.%d\n", vdate, vtime, pnum, pcat, pdis);
*/
GRIB_CHECK(grib_get_long(gh,"bitsPerValue", &lpar),0);
bitsPerValue = (int) lpar;
......@@ -717,16 +831,12 @@ void gribapiScanTimestep1(int streamID)
{
datetime0.date = vdate;
datetime0.time = vtime;
rdate = 0;
rtime = 0;
tunit = 0;
fcast = 0;
/*
rdate = gribRefDate(isec1);
rtime = gribRefTime(isec1);
tunit = cgribexGetTimeUnit(isec1);
fcast = gribTimeIsFC(isec1);
*/
GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
rdate = (int) lpar;
GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);
rtime = (int) lpar*100;
fcast = gribapiTimeIsFC(gh);
if ( fcast ) tunit = gribapiGetTimeUnits(gh);
}
else
{
......@@ -812,7 +922,6 @@ void gribapiScanTimestep1(int streamID)
{
taxisID = taxisCreate(TAXIS_ABSOLUTE);
taxis->type = TAXIS_ABSOLUTE;
taxis->unit = tunit;
}
taxis->vdate = datetime0.date;
......@@ -883,6 +992,7 @@ int gribapiScanTimestep2(int streamID)
int nrecords, nrecs, recID, rindex;
long recsize = 0;
int warn_numavg = TRUE;
int tsteptype;
int taxisID = -1;
TAXIS *taxis;
int vlistID;
......@@ -1026,24 +1136,18 @@ int gribapiScanTimestep2(int streamID)
level1 = (int) dlevel;
level2 = 0;
GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
vdate = (int) lpar;
GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);
vtime = (int) lpar*100;
gribapiGetValidityDateTime(gh, &vdate, &vtime);
if ( rindex == 0 )
{
if ( taxisInqType(taxisID) == TAXIS_RELATIVE )
{
taxis->type = TAXIS_RELATIVE;
taxis->rdate = 0;
taxis->rtime = 0;
taxis->unit = 0;
/*
taxis->rdate = gribRefDate(isec1);
taxis->rtime = gribRefTime(isec1);
taxis->unit = cgribexGetTimeUnit(isec1);
*/
GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
taxis->rdate = (int) lpar;
GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar), 0);
taxis->rtime = (int) lpar*100;
taxis->unit = gribapiGetTimeUnits(gh);
}
else
{
......@@ -1055,6 +1159,8 @@ int gribapiScanTimestep2(int streamID)
datetime0.date = vdate;
datetime0.time = vtime;
}
tsteptype = gribapiGetTsteptype(gh);
/*
if ( ISEC1_AvgNum )
{
......@@ -1089,7 +1195,7 @@ int gribapiScanTimestep2(int streamID)
{
char paramstr[32];
cdiParamToString(param, paramstr, sizeof(paramstr));
Warning(func, "Param=%s level=%d not defined at timestep 1!", param, level1);
Warning(func, "Param=%s level=%d not defined at timestep 1!", paramstr, level1);
return (CDI_EUFSTRUCT);
}
......@@ -1154,6 +1260,9 @@ int gribapiScanTimestep2(int streamID)
gridChangeType(gridID, GRID_TRAJECTORY);
}
*/
if ( tsteptype != vlistInqVarTsteptype(vlistID, varID) )
vlistDefVarTsteptype(vlistID, varID, tsteptype);
rindex++;
}
......@@ -1358,10 +1467,7 @@ int gribapiScanTimestep(int streamID)
level1 = (int) dlevel;
level2 = 0;
GRIB_CHECK(grib_get_long(gh, "dataDate", &lpar), 0);
vdate = (int) lpar;
GRIB_CHECK(grib_get_long(gh, "dataTime", &lpar),