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

added function gribapiGetGrid

parent 6f773b52
Hide whitespace changes
Inline Side-by-side
src/stream_gribapi.c
View file @ bc258201
......@@ -237,6 +237,7 @@ int gribapiGetTsteptype(grib_handle *gh)
#endif
#if defined (HAVE_LIBGRIB_API)
static
void gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
{
long lpar;
......@@ -269,112 +270,16 @@ void gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
#if defined (HAVE_LIBGRIB_API)
static
double grib1GetLevel(grib_handle *gh, int leveltype)
{
double dlevel;
GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0);
if ( leveltype == 100 ) dlevel *= 100;
if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
return (dlevel);
}
static
double grib2GetLevel(grib_handle *gh, int leveltype)
{
double dlevel;
GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0);
if ( leveltype == 100 ) dlevel *= 100;
if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
return (dlevel);
}
static
void gribapiAddRecord(int streamID, int param, grib_handle *gh,
long recsize, off_t position, int prec, int ztype)
void gribapiGetGrid(grib_handle *gh, grid_t *grid)
{
long editionNumber;
int gridtype;
int zaxistype;
int gridID = CDI_UNDEFID, varID;
int levelID = 0;
int tsID, recID;
int level1, level2;
int numavg;
int tsteptype;
int lbounds = 0;
record_t *record;
grid_t grid;
int vlistID;
stream_t *streamptr;
int leveltype;
double dlevel;
long lpar;
int status;
long numberOfPoints;
size_t datasize;
char name[256], longname[256], units[256];
size_t vlen;
streamptr = stream_to_pointer(streamID);
vlistID = streamInqVlist(streamID);
tsID = streamptr->curTsID;
recID = recordNewEntry(streamID, tsID);
record = &streamptr->tsteps[tsID].records[recID];
tsteptype = gribapiGetTsteptype(gh);
// numavg = ISEC1_AvgNum;
numavg = 0;
long numberOfPoints;
long lpar;
GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);
if ( editionNumber <= 1 )
{
status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
if ( status == 0 )
{
leveltype = (int) lpar;
dlevel = grib1GetLevel(gh, leveltype);
if ( leveltype == 99 ) leveltype++;
}
else
{
leveltype = 0;
dlevel = 0;
}
}
else
{
status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
if ( status == 0 )
{
leveltype = (int) lpar;
dlevel = grib2GetLevel(gh, leveltype);
if ( leveltype == 99 ) leveltype++;
}
else
{
leveltype = 0;
dlevel = 0;
}
}
level1 = (int) dlevel;
level2 = 0;
// fprintf(stderr, "param %d %d %d %d\n", param, level1, level2, leveltype);
(*record).size = recsize;
(*record).position = position;
(*record).param = param;
(*record).ilevel = level1;
(*record).ilevel2 = level2;
(*record).ltype = leveltype;
gridtype = gribapiGetGridType(gh);
/*
if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED )
......@@ -384,7 +289,7 @@ void gribapiAddRecord(int streamID, int param, grib_handle *gh,
ISEC4_NumValues = ISEC2_NumLon*ISEC2_NumLat;
}
*/
memset(&grid, 0, sizeof(grid_t));
memset(grid, 0, sizeof(grid_t));
GRIB_CHECK(grib_get_size(gh, "values", &datasize), 0);
GRIB_CHECK(grib_get_long(gh, "numberOfPoints", &numberOfPoints), 0);
......@@ -404,55 +309,55 @@ void gribapiAddRecord(int streamID, int param, grib_handle *gh,
if ( numberOfPoints != nlon*nlat )
Error("numberOfPoints (%d) and gridSize (%d) differ!",
(int)numberOfPoints, nlon*nlat);
grid.size = numberOfPoints;
grid.xsize = nlon;
grid.ysize = nlat;
grid.xinc = 0;
grid.yinc = 0;
grid.xdef = 0;
GRIB_CHECK(grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &grid.xfirst), 0);
GRIB_CHECK(grib_get_double(gh, "longitudeOfLastGridPointInDegrees", &grid.xlast), 0);
GRIB_CHECK(grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &grid.yfirst), 0);
GRIB_CHECK(grib_get_double(gh, "latitudeOfLastGridPointInDegrees", &grid.ylast), 0);
GRIB_CHECK(grib_get_double(gh, "iDirectionIncrementInDegrees", &grid.xinc), 0);
grid->size = numberOfPoints;
grid->xsize = nlon;
grid->ysize = nlat;
grid->xinc = 0;
grid->yinc = 0;
grid->xdef = 0;
GRIB_CHECK(grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &grid->xfirst), 0);
GRIB_CHECK(grib_get_double(gh, "longitudeOfLastGridPointInDegrees", &grid->xlast), 0);
GRIB_CHECK(grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &grid->yfirst), 0);
GRIB_CHECK(grib_get_double(gh, "latitudeOfLastGridPointInDegrees", &grid->ylast), 0);
GRIB_CHECK(grib_get_double(gh, "iDirectionIncrementInDegrees", &grid->xinc), 0);
if ( gridtype == GRID_LONLAT )
GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &grid.yinc), 0);
GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &grid->yinc), 0);
if ( IS_EQUAL(grid.xinc, GRIB_MISSING_DOUBLE) ) grid.xinc = 0;
if ( IS_EQUAL(grid->xinc, GRIB_MISSING_DOUBLE) ) grid->xinc = 0;
/* if ( IS_NOT_EQUAL(grid.xfirst, 0) || IS_NOT_EQUAL(grid.xlast, 0) ) */
/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
{
if ( grid.xsize > 1 )
if ( grid->xsize > 1 )
{
if ( (grid.xfirst > grid.xlast) && (grid.xfirst >= 180) ) grid.xfirst -= 360;
if ( (grid->xfirst > grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
if ( editionNumber <= 1 )
{
/* correct xinc if necessary */
if ( IS_EQUAL(grid.xfirst, 0) && grid.xlast > 354 )
if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
{
double xinc = 360. / grid.xsize;
double xinc = 360. / grid->xsize;
if ( fabs(grid.xinc-xinc) > 0.0 )
if ( fabs(grid->xinc-xinc) > 0.0 )
{
grid.xinc = xinc;
if ( CDI_Debug ) Message("set xinc to %g", grid.xinc);
grid->xinc = xinc;
if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
}
}
}
}
grid.xdef = 2;
grid->xdef = 2;
}
grid.ydef = 0;
/* if ( IS_NOT_EQUAL(grid.yfirst, 0) || IS_NOT_EQUAL(grid.ylast, 0) ) */
grid->ydef = 0;
/* if ( IS_NOT_EQUAL(grid->yfirst, 0) || IS_NOT_EQUAL(grid->ylast, 0) ) */
{
if ( grid.ysize > 1 )
if ( grid->ysize > 1 )
{
if ( editionNumber <= 1 )
{
}
}
grid.ydef = 2;
grid->ydef = 2;
}
break;
}
......@@ -465,60 +370,60 @@ void gribapiAddRecord(int streamID, int param, grib_handle *gh,
GRIB_CHECK(grib_get_long(gh, "Nj", &lpar), 0);
nlat = lpar;
grid.size = numberOfPoints;
grid->size = numberOfPoints;
grid.rowlon = (int *) malloc(nlat*sizeof(int));
grid->rowlon = (int *) malloc(nlat*sizeof(int));
pl = (long *) malloc(nlat*sizeof(long));
dummy = nlat;
GRIB_CHECK(grib_get_long_array(gh, "pl", pl, &dummy), 0);
for ( i = 0; i < nlat; ++i ) grid.rowlon[i] = pl[i];
for ( i = 0; i < nlat; ++i ) grid->rowlon[i] = pl[i];
free(pl);
grid.ysize = nlat;
grid.xinc = 0;
grid.yinc = 0;
grid.xdef = 0;
GRIB_CHECK(grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &grid.xfirst), 0);
GRIB_CHECK(grib_get_double(gh, "longitudeOfLastGridPointInDegrees", &grid.xlast), 0);
GRIB_CHECK(grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &grid.yfirst), 0);
GRIB_CHECK(grib_get_double(gh, "latitudeOfLastGridPointInDegrees", &grid.ylast), 0);
GRIB_CHECK(grib_get_double(gh, "iDirectionIncrementInDegrees", &grid.xinc), 0);
grid->ysize = nlat;
grid->xinc = 0;
grid->yinc = 0;
grid->xdef = 0;
GRIB_CHECK(grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &grid->xfirst), 0);
GRIB_CHECK(grib_get_double(gh, "longitudeOfLastGridPointInDegrees", &grid->xlast), 0);
GRIB_CHECK(grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &grid->yfirst), 0);
GRIB_CHECK(grib_get_double(gh, "latitudeOfLastGridPointInDegrees", &grid->ylast), 0);
GRIB_CHECK(grib_get_double(gh, "iDirectionIncrementInDegrees", &grid->xinc), 0);
if ( IS_EQUAL(grid.xinc, GRIB_MISSING_DOUBLE) ) grid.xinc = 0;
if ( IS_EQUAL(grid->xinc, GRIB_MISSING_DOUBLE) ) grid->xinc = 0;
/* if ( IS_NOT_EQUAL(grid.xfirst, 0) || IS_NOT_EQUAL(grid.xlast, 0) ) */
/* if ( IS_NOT_EQUAL(grid->xfirst, 0) || IS_NOT_EQUAL(grid->xlast, 0) ) */
{
if ( grid.xsize > 1 )
if ( grid->xsize > 1 )
{
if ( (grid.xfirst > grid.xlast) && (grid.xfirst >= 180) ) grid.xfirst -= 360;
if ( (grid->xfirst > grid->xlast) && (grid->xfirst >= 180) ) grid->xfirst -= 360;
if ( editionNumber <= 1 )
{
/* correct xinc if necessary */
if ( IS_EQUAL(grid.xfirst, 0) && grid.xlast > 354 )
if ( IS_EQUAL(grid->xfirst, 0) && grid->xlast > 354 )
{
double xinc = 360. / grid.xsize;
double xinc = 360. / grid->xsize;
if ( fabs(grid.xinc-xinc) > 0.0 )
if ( fabs(grid->xinc-xinc) > 0.0 )
{
grid.xinc = xinc;
if ( CDI_Debug ) Message("set xinc to %g", grid.xinc);
grid->xinc = xinc;
if ( CDI_Debug ) Message("set xinc to %g", grid->xinc);
}
}
}
}
grid.xdef = 2;
grid->xdef = 2;
}
grid.ydef = 0;
/* if ( IS_NOT_EQUAL(grid.yfirst, 0) || IS_NOT_EQUAL(grid.ylast, 0) ) */
grid->ydef = 0;
/* if ( IS_NOT_EQUAL(grid->yfirst, 0) || IS_NOT_EQUAL(grid->ylast, 0) ) */
{
if ( grid.ysize > 1 )
if ( grid->ysize > 1 )
{
if ( editionNumber <= 1 )
{
}
}
grid.ydef = 2;
grid->ydef = 2;
}
break;
}
......@@ -529,22 +434,22 @@ void gribapiAddRecord(int streamID, int param, grib_handle *gh,
Error("numberOfPoints (%d) and gridSize (%d) differ!",
ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);
grid.size = ISEC4_NumValues;
grid.xsize = ISEC2_NumLon;
grid.ysize = ISEC2_NumLat;
grid->size = ISEC4_NumValues;
grid->xsize = ISEC2_NumLon;
grid->ysize = ISEC2_NumLat;
grid.lcc_xinc = ISEC2_Lambert_dx;
grid.lcc_yinc = ISEC2_Lambert_dy;
grid.lcc_originLon = ISEC2_FirstLon * 0.001;
grid.lcc_originLat = ISEC2_FirstLat * 0.001;
grid.lcc_lonParY = ISEC2_Lambert_Lov * 0.001;
grid.lcc_lat1 = ISEC2_Lambert_LatS1 * 0.001;
grid.lcc_lat2 = ISEC2_Lambert_LatS2 * 0.001;
grid.lcc_projflag = ISEC2_Lambert_ProjFlag;
grid.lcc_scanflag = ISEC2_ScanFlag;
grid->lcc_xinc = ISEC2_Lambert_dx;
grid->lcc_yinc = ISEC2_Lambert_dy;
grid->lcc_originLon = ISEC2_FirstLon * 0.001;
grid->lcc_originLat = ISEC2_FirstLat * 0.001;
grid->lcc_lonParY = ISEC2_Lambert_Lov * 0.001;
grid->lcc_lat1 = ISEC2_Lambert_LatS1 * 0.001;
grid->lcc_lat2 = ISEC2_Lambert_LatS2 * 0.001;
grid->lcc_projflag = ISEC2_Lambert_ProjFlag;
grid->lcc_scanflag = ISEC2_ScanFlag;
grid.xdef = 0;
grid.ydef = 0;
grid->xdef = 0;
grid->ydef = 0;
break;
}
......@@ -554,22 +459,22 @@ void gribapiAddRecord(int streamID, int param, grib_handle *gh,
size_t len = 256;
char typeOfPacking[256];
GRIB_CHECK(grib_get_string(gh, "packingType", typeOfPacking, &len), 0);
grid.lcomplex = 0;
if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid.lcomplex = 1;
grid->lcomplex = 0;
if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid->lcomplex = 1;
grid.size = datasize;
grid->size = datasize;
GRIB_CHECK(grib_get_long(gh, "J", &lpar), 0);
grid.trunc = lpar;
grid->trunc = lpar;
break;
}
case GRID_GME:
{
grid.size = numberOfPoints;
if ( grib_get_long(gh, "nd", &lpar) == 0 ) grid.nd = lpar;
if ( grib_get_long(gh, "Ni", &lpar) == 0 ) grid.ni = lpar;
if ( grib_get_long(gh, "n2", &lpar) == 0 ) grid.ni2 = lpar;
if ( grib_get_long(gh, "n3", &lpar) == 0 ) grid.ni3 = lpar;
grid->size = numberOfPoints;
if ( grib_get_long(gh, "nd", &lpar) == 0 ) grid->nd = lpar;
if ( grib_get_long(gh, "Ni", &lpar) == 0 ) grid->ni = lpar;
if ( grib_get_long(gh, "n2", &lpar) == 0 ) grid->ni2 = lpar;
if ( grib_get_long(gh, "n3", &lpar) == 0 ) grid->ni3 = lpar;
break;
}
......@@ -579,15 +484,15 @@ void gribapiAddRecord(int streamID, int param, grib_handle *gh,
size_t len = sizeof(reference_link);
reference_link[0] = 0;
grid.size = numberOfPoints;
grid->size = numberOfPoints;
if ( grib_get_long(gh, "numberOfGridUsed", &lpar) == 0 )
{
grid.number = lpar;
if ( grib_get_long(gh, "numberOfGridInReference", &lpar) == 0 ) grid.position = lpar;
grid->number = lpar;
if ( grib_get_long(gh, "numberOfGridInReference", &lpar) == 0 ) grid->position = lpar;
if ( grib_get_string(gh, "gridDescriptionFile", reference_link, &len) == 0 )
{
if ( strncmp(reference_link, "file://", 7) == 0 )
grid.reference = strdupx(reference_link);
grid->reference = strdupx(reference_link);
}
}
......@@ -600,17 +505,17 @@ void gribapiAddRecord(int streamID, int param, grib_handle *gh,
if ( grib_get_long(gh, "Ni", &lpar) == 0 ) nlon = lpar;
if ( grib_get_long(gh, "Nj", &lpar) == 0 ) nlat = lpar;
grid.size = numberOfPoints;
grid->size = numberOfPoints;
if ( nlon > 0 && nlat > 0 && nlon*nlat == grid.size )
if ( nlon > 0 && nlat > 0 && nlon*nlat == grid->size )
{
grid.xsize = nlon;
grid.ysize = nlat;
grid->xsize = nlon;
grid->ysize = nlat;
}
else
{
grid.xsize = 0;
grid.ysize = 0;
grid->xsize = 0;
grid->ysize = 0;
}
break;
......@@ -622,21 +527,130 @@ void gribapiAddRecord(int streamID, int param, grib_handle *gh,
}
}
grid.isRotated = FALSE;
grid->isRotated = FALSE;
if ( gribapiGetIsRotated(gh) )
{
grid.isRotated = TRUE;
GRIB_CHECK(grib_get_double(gh, "latitudeOfSouthernPoleInDegrees", &grid.ypole), 0);
GRIB_CHECK(grib_get_double(gh, "longitudeOfSouthernPoleInDegrees", &grid.xpole), 0);
GRIB_CHECK(grib_get_double(gh, "angleOfRotation", &grid.angle), 0);
grid->isRotated = TRUE;
GRIB_CHECK(grib_get_double(gh, "latitudeOfSouthernPoleInDegrees", &grid->ypole), 0);
GRIB_CHECK(grib_get_double(gh, "longitudeOfSouthernPoleInDegrees", &grid->xpole), 0);
GRIB_CHECK(grib_get_double(gh, "angleOfRotation", &grid->angle), 0);
/* change from south to north pole */
grid.ypole = -grid.ypole;
grid.xpole = grid.xpole - 180;
grid->ypole = -grid->ypole;
grid->xpole = grid->xpole - 180;
}
grid->xvals = NULL;
grid->yvals = NULL;
grid->type = gridtype;
}
#endif
#if defined (HAVE_LIBGRIB_API)
static
double grib1GetLevel(grib_handle *gh, int leveltype)
{
double dlevel;
GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0);
if ( leveltype == 100 ) dlevel *= 100;
if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
return (dlevel);
}
static
double grib2GetLevel(grib_handle *gh, int leveltype)
{
double dlevel;
GRIB_CHECK(grib_get_double(gh, "level", &dlevel), 0);
if ( leveltype == 100 ) dlevel *= 100;
if ( dlevel < -2.e9 || dlevel > 2.e9 ) dlevel = 0;
return (dlevel);
}
static
void gribapiAddRecord(int streamID, int param, grib_handle *gh,
long recsize, off_t position, int prec, int ztype)
{
long editionNumber;
int zaxistype;
int gridID = CDI_UNDEFID, varID;
int levelID = 0;
int tsID, recID;
int level1, level2;
int numavg;
int tsteptype;
int lbounds = 0;
record_t *record;
grid_t grid;
int vlistID;
stream_t *streamptr;
int leveltype;
double dlevel;
long lpar;
int status;
char name[256], longname[256], units[256];
size_t vlen;
streamptr = stream_to_pointer(streamID);
vlistID = streamInqVlist(streamID);
tsID = streamptr->curTsID;
recID = recordNewEntry(streamID, tsID);
record = &streamptr->tsteps[tsID].records[recID];
tsteptype = gribapiGetTsteptype(gh);
// numavg = ISEC1_AvgNum;
numavg = 0;
GRIB_CHECK(grib_get_long(gh, "editionNumber", &editionNumber), 0);
if ( editionNumber <= 1 )
{
status = grib_get_long(gh, "indicatorOfTypeOfLevel", &lpar);
if ( status == 0 )
{
leveltype = (int) lpar;
dlevel = grib1GetLevel(gh, leveltype);
if ( leveltype == 99 ) leveltype++;
}
else
{
leveltype = 0;
dlevel = 0;
}
}
else
{
status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar);
if ( status == 0 )
{
leveltype = (int) lpar;
dlevel = grib2GetLevel(gh, leveltype);
if ( leveltype == 99 ) leveltype++;
}
else
{
leveltype = 0;
dlevel = 0;
}
}
grid.xvals = NULL;
grid.yvals = NULL;
grid.type = gridtype;
level1 = (int) dlevel;
level2 = 0;
// fprintf(stderr, "param %d %d %d %d\n", param, level1, level2, leveltype);
(*record).size = recsize;
(*record).position = position;
(*record).param = param;
(*record).ilevel = level1;
(*record).ilevel2 = level2;
(*record).ltype = leveltype;
gribapiGetGrid(gh, &grid);
gridID = varDefGrid(vlistID, grid, 0);
......
Supports Markdown
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