Skip to content

GitLab

Commit 548d9b4a authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

added missing value support for cosmo GRIB files

parent 97839498
Hide whitespace changes
Inline Side-by-side
src/cdi_util.c
View file @ 548d9b4a
......@@ -18,8 +18,8 @@ int cdiEncodeParam(int pnum, int pcat, int pdis)
int param;
int inum;
if ( pcat < 0 ) pcat = 255;
if ( pdis < 0 ) pdis = 255;
if ( pcat < 0 || pcat > 255 ) pcat = 255;
if ( pdis < 0 || pdis > 255 ) pdis = 255;
inum = pnum;
if ( pnum < 0 ) inum = -pnum;
......
src/stream_cdf.c
View file @ 548d9b4a
......@@ -5858,23 +5858,21 @@ int cdfInqContents(int streamID)
pname = vlistInqVarNamePtr(vlistID, varID);
len = strlen(pname);
if ( len > 3 && isdigit((int) pname[3]) )
if ( memcmp("var", pname, 3) == 0 )
{
vlistDefVarCode(vlistID, varID, atoi(pname+3));
vlistDestroyVarName(vlistID, varID);
}
}
if ( vlistInqVarCode(vlistID, varID) == -varID-1 )
{
char *pname;
pname = vlistInqVarNamePtr(vlistID, varID);
len = strlen(pname);
if ( len > 4 && isdigit((int) pname[4]) )
if ( memcmp("code", pname, 4) == 0 )
{
vlistDefVarCode(vlistID, varID, atoi(pname+4));
vlistDestroyVarName(vlistID, varID);
}
{
if ( memcmp("var", pname, 3) == 0 )
{
vlistDefVarCode(vlistID, varID, atoi(pname+3));
vlistDestroyVarName(vlistID, varID);
}
}
else if ( len > 4 && isdigit((int) pname[4]) )
{
if ( memcmp("code", pname, 4) == 0 )
{
vlistDefVarCode(vlistID, varID, atoi(pname+4));
vlistDestroyVarName(vlistID, varID);
}
}
}
}
......
src/stream_cgribex.c
View file @ 548d9b4a
......@@ -204,8 +204,8 @@ int cgribexGetTimeUnit(int *isec1)
static
void cgribexAddRecord(int streamID, int param, int *isec1, int *isec2, double *fsec2,
int *isec4, long recsize, off_t position, int prec, int ztype)
void cgribexAddRecord(int streamID, int param, int *isec1, int *isec2, double *fsec2, double *fsec3,
int *isec4, long recsize, off_t position, int prec, int ztype, int lmv)
{
static char func[] = "cgribexAddRecord";
int gridtype;
......@@ -436,6 +436,8 @@ void cgribexAddRecord(int streamID, int param, int *isec1, int *isec2, double *f
varDefZtype(varID, ztype);
if ( lmv ) varDefMissval(varID, FSEC3_MissVal);
if ( varInqInst(varID) == CDI_UNDEFID )
{
int center, subcenter, instID;
......@@ -475,16 +477,63 @@ void cgribexAddRecord(int streamID, int param, int *isec1, int *isec2, double *f
varID, param, zaxistype, gridID, levelID);
}
static
void MCH_get_undef(int *isec1, double *undef_pds, double *undef_eps)
{
/* 2010-01-13: Oliver Fuhrer */
if ( ISEC1_CenterID == 215 ) {
if (isec1[34] != 0 && isec1[34] != 255) {
if (isec1[34] & 2) {
if (isec1[34] & 1) {
*undef_pds = -0.99*pow(10.0,-isec1[35]);
} else {
*undef_pds = +0.99*pow(10.0,-isec1[35]);
}
*undef_eps = pow(10.0,-isec1[35]-1);
} else {
if (isec1[34] & 1) {
*undef_pds = -0.99*pow(10.0,+isec1[35]);
} else {
*undef_pds = +0.99*pow(10.0,+isec1[35]);
}
*undef_eps = pow(10.0,isec1[35]-1);
}
}
}
}
static
void cgribexDecodeHeader(int *isec0, int *isec1, int *isec2, double *fsec2,
int *isec3, double *fsec3, int *isec4, double *fsec4,
int *gribbuffer, int recsize, int *lmv)
{
int iret = 0, ipunp = 0, iword = 0;
gribExDP(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
ipunp, (int *) gribbuffer, recsize, &iword, "J", &iret);
*lmv = 0;
if ( ISEC1_CenterID == 215 && (isec1[34] != 0 && isec1[34] != 255) )
{
double undef_pds, undef_eps;
MCH_get_undef(isec1, &undef_pds, &undef_eps);
FSEC3_MissVal = undef_pds;
*lmv = 1;
}
}
void cgribexScanTimestep1(int streamID)
{
static char func[] = "cgribexScanTimestep1";
int *isec0, *isec1, *isec2, *isec3, *isec4;
double fsec2[512], fsec3[2], *fsec4 = NULL;
int lmv = 0;
off_t recpos = 0;
unsigned char *gribbuffer = NULL;
long buffersize = 0;
int iret = 0, ipunp = 0, iword = 0;
int status;
int fileID;
int param = 0;
......@@ -560,8 +609,8 @@ void cgribexScanTimestep1(int streamID)
}
}
gribExDP(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
ipunp, (int *) gribbuffer, recsize, &iword, "J", &iret);
cgribexDecodeHeader(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
(int *) gribbuffer, recsize, &lmv);
param = cdiEncodeParam(ISEC1_Parameter, ISEC1_CodeTable, 255);
if ( ISEC1_LevelType == 100 ) ISEC1_Level1 *= 100;
......@@ -644,7 +693,8 @@ void cgribexScanTimestep1(int streamID)
if ( CDI_Debug )
Message(func, "%4d %8d %4d %8d %8d %6d", nrecs, (int)recpos, param, level1, vdate, vtime);
cgribexAddRecord(streamID, param, isec1, isec2, fsec2, isec4, recsize, recpos, prec, ztype);
cgribexAddRecord(streamID, param, isec1, isec2, fsec2, fsec3,
isec4, recsize, recpos, prec, ztype, lmv);
}
streamptr->rtsteps = 1;
......@@ -717,10 +767,10 @@ int cgribexScanTimestep2(int streamID)
static char func[] = "cgribexScanTimestep2";
int *isec0, *isec1, *isec2, *isec3, *isec4;
double fsec2[512], fsec3[2], *fsec4 = NULL;
int lmv = 0;
off_t recpos = 0;
unsigned char *gribbuffer = NULL;
long buffersize = 0;
int iret = 0, ipunp = 0, iword = 0;
int status;
int fileID;
int param = 0;
......@@ -815,8 +865,8 @@ int cgribexScanTimestep2(int streamID)
}
}
gribExDP(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
ipunp, (int *) gribbuffer, recsize, &iword, "J", &iret);
cgribexDecodeHeader(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
(int *) gribbuffer, recsize, &lmv);
param = cdiEncodeParam(ISEC1_Parameter, ISEC1_CodeTable, 255);
if ( ISEC1_LevelType == 100 ) ISEC1_Level1 *= 100;
......@@ -984,11 +1034,11 @@ int cgribexScanTimestep(int streamID)
static char func[] = "cgribexScanTimestep";
int *isec0, *isec1, *isec2, *isec3, *isec4;
double fsec2[512], fsec3[2], *fsec4 = NULL;
int lmv = 0;
long recsize = 0;
off_t recpos = 0;
unsigned char *gribbuffer;
long buffersize = 0;
int iret = 0, ipunp = 0, iword = 0;
int status;
int fileID;
int param = 0;
......@@ -1084,8 +1134,8 @@ int cgribexScanTimestep(int streamID)
}
}
gribExDP(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
ipunp, (int *) gribbuffer, recsize, &iword, "J", &iret);
cgribexDecodeHeader(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, fsec4,
(int *) gribbuffer, recsize, &lmv);
param = cdiEncodeParam(ISEC1_Parameter, ISEC1_CodeTable, 255);
if ( ISEC1_LevelType == 100 ) ISEC1_Level1 *= 100;
......@@ -1256,7 +1306,7 @@ int cgribexScanTimestep(int streamID)
int cgribexDecode(unsigned char *gribbuffer, int gribsize, double *data, int gridsize,
int unreduced, int *nmiss, int *zip, double *missval)
int unreduced, int *nmiss, int *zip, double missval)
{
static char func[] = "cgribexDecode";
int status = 0;
......@@ -1273,45 +1323,44 @@ int cgribexDecode(unsigned char *gribbuffer, int gribsize, double *data, int gri
if ( unreduced ) strcpy(hoper, "R");
else strcpy(hoper, "D");
FSEC3_MissVal = cdiDefaultMissval;
*missval = cdiDefaultMissval;
FSEC3_MissVal = missval;
if ( (izip = gribGetZip(gribsize, gribbuffer, &unzipsize)) > 0 )
{
*zip = izip;
if ( izip == 128 ) /* szip */
{
unsigned char *itmpbuffer = NULL;
size_t itmpbuffersize = 0;
*zip = izip;
if ( izip == 128 ) /* szip */
{
unsigned char *itmpbuffer = NULL;
size_t itmpbuffersize = 0;
if ( unzipsize < (long) gribsize )
{
fprintf(stderr, "Decompressed size smaller than compressed size (in %d; out %ld)!\n",
gribsize, unzipsize);
return (status);
}
if ( unzipsize < (long) gribsize )
{
fprintf(stderr, "Decompressed size smaller than compressed size (in %d; out %ld)!\n",
gribsize, unzipsize);
return (status);
}
if ( itmpbuffersize < (size_t) gribsize )
{
itmpbuffersize = gribsize;
itmpbuffer = (unsigned char *) realloc(itmpbuffer, itmpbuffersize);
}
if ( itmpbuffersize < (size_t) gribsize )
{
itmpbuffersize = gribsize;
itmpbuffer = (unsigned char *) realloc(itmpbuffer, itmpbuffersize);
}
memcpy(itmpbuffer, gribbuffer, itmpbuffersize);
memcpy(itmpbuffer, gribbuffer, itmpbuffersize);
unzipsize += 100; /* need 0 to 1 bytes for rounding of bds */
unzipsize += 100; /* need 0 to 1 bytes for rounding of bds */
gribsize = gribUnzip(gribbuffer, unzipsize, itmpbuffer, gribsize);
gribsize = gribUnzip(gribbuffer, unzipsize, itmpbuffer, gribsize);
if ( gribsize <= 0 )
Error(func, "Decompression problem!\n");
free(itmpbuffer);
}
else
{
Error(func, "Decompression for %d not implemented!\n", izip);
}
if ( gribsize <= 0 )
Error(func, "Decompression problem!\n");
free(itmpbuffer);
}
else
{
Error(func, "Decompression for %d not implemented!\n", izip);
}
}
gribExDP(isec0, isec1, isec2, fsec2, isec3, fsec3, isec4, data,
......@@ -1322,36 +1371,20 @@ int cgribexDecode(unsigned char *gribbuffer, int gribsize, double *data, int gri
else
*nmiss = 0;
/* 2010-01-13: Oliver Fuhrer */
if (ISEC1_CenterID == 215 ) {
int i;
double undef_pds, undef_eps;
if ( ISEC1_CenterID == 215 && (isec1[34] != 0 && isec1[34] != 255) )
{
double undef_pds, undef_eps;
int i;
if (isec1[34] != 0 && isec1[34] != 255) {
if (isec1[34] & 2) {
if (isec1[34] & 1) {
undef_pds = -0.99*pow(10.0,-isec1[35]);
} else {
undef_pds = +0.99*pow(10.0,-isec1[35]);
}
undef_eps = pow(10.0,-isec1[35]-1);
} else {
if (isec1[34] & 1) {
undef_pds = -0.99*pow(10.0,+isec1[35]);
} else {
undef_pds = +0.99*pow(10.0,+isec1[35]);
}
undef_eps = pow(10.0,isec1[35]-1);
}
MCH_get_undef(isec1, &undef_pds, &undef_eps);
*nmiss = 0;
for ( i = 0; i < gridsize; i++ )
if ( (abs(data[i]-undef_pds) < undef_eps) || IS_EQUAL(data[i],FSEC3_MissVal) ) {
data[i] = undef_pds;
data[i] = missval;
(*nmiss)++;
}
FSEC3_MissVal = undef_pds;
*missval = undef_pds;
}
}
#else
Error(func, "CGRIBEX support not compiled in!");
#endif
......
src/stream_cgribex.h
View file @ 548d9b4a
......@@ -6,7 +6,7 @@ int cgribexScanTimestep2(int streamID);
int cgribexScanTimestep(int streamID);
int cgribexDecode(unsigned char *gribbuffer, int gribsize, double *data, int gridsize,
int unreduced, int *nmiss, int *zip, double *missval);
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 numavg,
long datasize, const double *data, int nmiss, unsigned char *gribbuffer, size_t gribbuffersize);
......
src/stream_grb.c
View file @ 548d9b4a
......@@ -67,7 +67,7 @@ int grbDefRecord(int streamID)
static
int grbDecode(int filetype, unsigned char *gribbuffer, int gribsize, double *data, int gridsize,
int unreduced, int *nmiss, int *zip, double *missval)
int unreduced, int *nmiss, int *zip, double missval)
{
//static char func[] = "grbDecode";
int status = 0;
......@@ -123,9 +123,10 @@ int grbReadRecord(int streamID, double *data, int *nmiss)
fileRead(fileID, gribbuffer, (size_t) recsize);
grbDecode(filetype, gribbuffer, recsize, data, gridsize, streamptr->unreduced, nmiss, &zip, &missval);
printf("missval = %g\n", vlistInqVarMissval(vlistID, varID));
missval = vlistInqVarMissval(vlistID, varID);
if ( *nmiss > 0 ) vlistDefVarMissval(vlistID, varID, missval);
grbDecode(filetype, gribbuffer, recsize, data, gridsize, streamptr->unreduced, nmiss, &zip, missval);
streamptr->tsteps[tsID].records[recID].zip = zip;
......@@ -310,10 +311,10 @@ void grbReadVarDP(int streamID, int varID, double *data, int *nmiss)
fileRead(fileID, gribbuffer, recsize);
grbDecode(filetype, gribbuffer, recsize, &data[levelID*gridsize], gridsize,
streamptr->unreduced, &imiss, &zip, &missval);
missval = vlistInqVarMissval(vlistID, varID);
if ( imiss > 0 ) vlistDefVarMissval(vlistID, varID, missval);
grbDecode(filetype, gribbuffer, recsize, &data[levelID*gridsize], gridsize,
streamptr->unreduced, &imiss, &zip, missval);
*nmiss += imiss;
......@@ -371,9 +372,9 @@ void grbReadVarSliceDP(int streamID, int varID, int levelID, double *data, int *
fileRead(fileID, gribbuffer, recsize);
grbDecode(filetype, gribbuffer, recsize, data, gridsize, streamptr->unreduced, nmiss, &zip, &missval);
missval = vlistInqVarMissval(vlistID, varID);
if ( *nmiss > 0 ) vlistDefVarMissval(vlistID, varID, missval);
grbDecode(filetype, gribbuffer, recsize, data, gridsize, streamptr->unreduced, nmiss, &zip, missval);
fileSetPos(fileID, currentfilepos, SEEK_SET);
......
src/stream_gribapi.c
View file @ 548d9b4a
......@@ -1521,7 +1521,7 @@ int gribapiScanTimestep(int streamID)
int gribapiDecode(unsigned char *gribbuffer, int gribsize, double *data, int gridsize,
int unreduced, int *nmiss, int *zip, double *missval)
int unreduced, int *nmiss, int *zip, double missval)
{
static char func[] = "gribapiDecode";
int status = 0;
......@@ -1544,8 +1544,7 @@ int gribapiDecode(unsigned char *gribbuffer, int gribsize, double *data, int gri
recsize = gribsize;
gh = grib_handle_new_from_message(NULL, (void *) gribbuffer, recsize);
GRIB_CHECK(grib_set_double(gh, "missingValue", GRIB_MISSVAL), 0);
*missval = GRIB_MISSVAL;
GRIB_CHECK(grib_set_double(gh, "missingValue", missval), 0);
GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);
......
src/stream_gribapi.h
View file @ 548d9b4a
......@@ -6,7 +6,7 @@ int gribapiScanTimestep2(int streamID);
int gribapiScanTimestep(int streamID);
int gribapiDecode(unsigned char *gribbuffer, int gribsize, double *data, int gridsize,
int unreduced, int *nmiss, int *zip, double *missval);
int unreduced, int *nmiss, int *zip, double missval);
size_t gribapiEncode(int varID, int levelID, int vlistID, int gridID, int zaxisID, int date, int time, int numavg,
long datasize, const double *data, int nmiss, unsigned char *gribbuffer, size_t gribbuffersize,
......
src/varscan.c
View file @ 548d9b4a
......@@ -45,6 +45,8 @@ typedef struct
int tableID;
int ztype;
int zlevel;
int lmissval;
double missval;
char *name;
char *longname;
char *units;
......@@ -76,6 +78,8 @@ void paramInitEntry(int varID, int param)
vartable[varID].tableID = UNDEFID;
vartable[varID].ztype = COMPRESS_NONE;
vartable[varID].zlevel = 1;
vartable[varID].lmissval = 0;
vartable[varID].missval = 0;
vartable[varID].name = NULL;
vartable[varID].longname = NULL;
vartable[varID].units = NULL;
......@@ -522,6 +526,8 @@ void cdiGenVars(int streamID)
vlistDefVarTimaccu(vlistID, varID, timaccu);
vlistDefVarZtype(vlistID, varID, ztype);
if ( vartable[varid].lmissval ) vlistDefVarMissval(vlistID, varID, vartable[varid].missval);
if ( vartable[varid].name ) vlistDefVarName(vlistID, varID, vartable[varid].name);
if ( vartable[varid].longname ) vlistDefVarLongname(vlistID, varID, vartable[varid].longname);
if ( vartable[varid].units ) vlistDefVarUnits(vlistID, varID, vartable[varid].units);
......@@ -813,6 +819,13 @@ int varDefZaxis(int vlistID, int zaxistype, int nlevels, double *levels, int lbo
}
void varDefMissval(int varID, double missval)
{
vartable[varID].lmissval = 1;
vartable[varID].missval = missval;
}
void varDefZtype(int varID, int ztype)
{
vartable[varID].ztype = ztype;
......
src/varscan.h
View file @ 548d9b4a
......@@ -18,6 +18,7 @@ int varDefZaxis(int vlistID, int zaxistype, int nlevels, double *levels, int lb
double *levels1, double *levels2, int vctsize, double *vct, char *name,
char *longname, char *units, int prec, int mode, int ltype);
void varDefMissval(int varID, double missval);
void varDefZtype(int varID, int ztype);
void varDefZlevel(int varID, int zlevel);
void varDefInst(int varID, int instID);
......
src/vlist_var.c
View file @ 548d9b4a
......@@ -287,12 +287,17 @@ void vlistDefVarCode(int vlistID, int varID, int code)
{
static char func[] = "vlistDefVarCode";
vlist_t *vlistptr;
int param, pnum, pcat, pdis;
vlistptr = vlist_to_pointer(vlistID);
vlistCheckVarID(func, vlistID, varID);
vlistptr->vars[varID].param = cdiEncodeParam(code, 255, 255);
param = vlistptr->vars[varID].param;
cdiDecodeParam(param, &pnum, &pcat, &pdis);
vlistptr->vars[varID].param = cdiEncodeParam(code, pcat, pdis);
}
......@@ -888,6 +893,18 @@ void vlistDefVarTable(int vlistID, int varID, int tableID)
vlistptr = vlist_to_pointer(vlistID);
vlistptr->vars[varID].tableID = tableID;
{
int param, pnum, pcat, pdis;
int tablenum;
tablenum = tableInqNum(tableID);
param = vlistptr->vars[varID].param;
cdiDecodeParam(param, &pnum, &pcat, &pdis);
vlistptr->vars[varID].param = cdiEncodeParam(pnum, tablenum, pdis);
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Imprint | Privacy Policy