Commit 1c3e51bc authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

grb2 read: added support for projection CDI_PROJ_RLL.

parent 55a27d6e
2016-08-01 Uwe Schulzweida
* grb write: added support for projection CDI_PROJ_RLL
* ieg read: added support for projection CDI_PROJ_RLL
* grb, grb2, ieg read: added support for projection CDI_PROJ_RLL
2016-07-30 Uwe Schulzweida
......
......@@ -429,7 +429,11 @@ int gribapiGetGridType(grib_handle *gh)
gridtype = ( gribGetLong(gh, "Ni") == (long) GRIB_MISSING_LONG ) ? GRID_GAUSSIAN_REDUCED : GRID_GAUSSIAN;
break;
#ifdef TEST_PROJECTION
case GRIB2_GTYPE_LATLON_ROT: gridtype = GRID_PROJECTION; break;
#else
case GRIB2_GTYPE_LATLON_ROT: gridtype = GRID_LONLAT; break;
#endif
case GRIB2_GTYPE_LCC: gridtype = GRID_LCC; break;
case GRIB2_GTYPE_SPECTRAL: gridtype = GRID_SPECTRAL; break;
case GRIB2_GTYPE_GME: gridtype = GRID_GME; break;
......@@ -450,6 +454,7 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
{
long editionNumber = gribEditionNumber(gh);
int gridtype = gribapiGetGridType(gh);
int projtype = (gridtype == GRID_PROJECTION && gribapiGetIsRotated(gh)) ? CDI_PROJ_RLL : CDI_UNDEFID;
/*
if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED )
{
......@@ -466,10 +471,9 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
long numberOfPoints;
FAIL_ON_GRIB_ERROR(grib_get_long, gh, "numberOfPoints", &numberOfPoints);
switch (gridtype)
{
case GRID_LONLAT:
case GRID_GAUSSIAN:
long lpar;
if ( gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN || projtype == CDI_PROJ_RLL )
{
long nlon, nlat;
FAIL_ON_GRIB_ERROR(grib_get_long, gh, "Ni", &nlon);
......@@ -477,7 +481,6 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
if ( gridtype == GRID_GAUSSIAN )
{
long lpar;
FAIL_ON_GRIB_ERROR(grib_get_long, gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar);
grid->np = (int)lpar;
}
......@@ -485,7 +488,6 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
if ( numberOfPoints != nlon*nlat )
Error("numberOfPoints (%ld) and gridSize (%ld) differ!", numberOfPoints, nlon*nlat);
/* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
grid->size = (int)numberOfPoints;
grid->x.size = (int)nlon;
grid->y.size = (int)nlat;
......@@ -518,7 +520,6 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
if ( IS_EQUAL(grid->x.first, 0) && grid->x.last > 354 )
{
double xinc = 360. / grid->x.size;
if ( fabs(grid->x.inc-xinc) > 0.0 )
{
grid->x.inc = xinc;
......@@ -540,20 +541,15 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
}
grid->y.flag = 2;
}
break;
}
case GRID_GAUSSIAN_REDUCED:
else if ( gridtype == GRID_GAUSSIAN_REDUCED )
{
long lpar;
FAIL_ON_GRIB_ERROR(grib_get_long, gh, "numberOfParallelsBetweenAPoleAndTheEquator", &lpar);
/* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
grid->np = (int)lpar;
FAIL_ON_GRIB_ERROR(grib_get_long, gh, "Nj", &lpar);
/* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
int nlat = (int)lpar;
/* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
grid->size = (int)numberOfPoints;
grid->nrowlon = nlat;
......@@ -561,7 +557,6 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
long *pl = (long *) Malloc((size_t)nlat * sizeof (long));
size_t dummy = (size_t)nlat;
FAIL_ON_GRIB_ERROR(grib_get_long_array, gh, "pl", pl, &dummy);
/* FIXME: assert(pl[i] >= INT_MIN && pl[i] <= INT_MIN) */
for ( int i = 0; i < nlat; ++i ) grid->rowlon[i] = (int)pl[i];
Free(pl);
......@@ -589,7 +584,6 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
if ( IS_EQUAL(grid->x.first, 0) && grid->x.last > 354 )
{
double xinc = 360. / grid->x.size;
if ( fabs(grid->x.inc-xinc) > 0.0 )
{
grid->x.inc = xinc;
......@@ -611,13 +605,10 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
}
grid->y.flag = 2;
}
break;
}
case GRID_LCC:
else if ( gridtype == GRID_LCC )
{
int nlon, nlat;
long lpar;
FAIL_ON_GRIB_ERROR(grib_get_long, gh, "Nx", &lpar);
nlon = (int)lpar;
FAIL_ON_GRIB_ERROR(grib_get_long, gh, "Ny", &lpar);
......@@ -648,10 +639,8 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
grid->x.flag = 0;
grid->y.flag = 0;
break;
}
case GRID_SPECTRAL:
else if ( gridtype == GRID_SPECTRAL )
{
size_t len = 256;
char typeOfPacking[256];
......@@ -659,32 +648,20 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
grid->lcomplex = 0;
if ( strncmp(typeOfPacking, "spectral_complex", len) == 0 ) grid->lcomplex = 1;
/* FIXME: assert(datasize >= INT_MIN && datasize <= INT_MAX) */
grid->size = (int)datasize;
long lpar;
FAIL_ON_GRIB_ERROR(grib_get_long, gh, "J", &lpar);
/* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
grid->trunc = (int)lpar;
break;
}
case GRID_GME:
else if ( gridtype == GRID_GME )
{
/* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
grid->size = (int)numberOfPoints;
long lpar;
/* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
if ( grib_get_long(gh, "nd", &lpar) == 0 ) grid->gme.nd = (int)lpar;
/* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
if ( grib_get_long(gh, "Ni", &lpar) == 0 ) grid->gme.ni = (int)lpar;
/* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
if ( grib_get_long(gh, "n2", &lpar) == 0 ) grid->gme.ni2 = (int)lpar;
/* FIXME: assert(lpar >= INT_MIN && lpar <= INT_MAX) */
if ( grib_get_long(gh, "n3", &lpar) == 0 ) grid->gme.ni3 = (int)lpar;
break;
}
case GRID_UNSTRUCTURED:
else if ( gridtype == GRID_UNSTRUCTURED )
{
unsigned char uuid[CDI_UUID_SIZE];
/*
......@@ -692,15 +669,11 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
size_t len = sizeof(reference_link);
reference_link[0] = 0;
*/
/* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
grid->size = (int)numberOfPoints;
long lpar;
if ( grib_get_long(gh, "numberOfGridUsed", &lpar) == 0 )
{
/* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
grid->number = (int)lpar;
/* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
if ( grib_get_long(gh, "numberOfGridInReference", &lpar) == 0 )
grid->position = (int)lpar;
/*
......@@ -716,18 +689,13 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
memcpy(grid->uuid, uuid, CDI_UUID_SIZE);
}
}
break;
}
case GRID_GENERIC:
else if ( gridtype == GRID_GENERIC )
{
int nlon = 0, nlat = 0;
long lpar;
/* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
if ( grib_get_long(gh, "Ni", &lpar) == 0 ) nlon = (int)lpar;
/* FIXME: assert(lpar <= INT_MAX && lpar >= INT_MIN) */
if ( grib_get_long(gh, "Nj", &lpar) == 0 ) nlat = (int)lpar;
/* FIXME: assert(numberOfPoints <= INT_MAX && numberOfPoints >= INT_MIN) */
grid->size = (int)numberOfPoints;
if ( nlon > 0 && nlat > 0 && nlon*nlat == grid->size )
......@@ -740,27 +708,22 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
grid->x.size = 0;
grid->y.size = 0;
}
break;
}
default:
else
{
Error("Unsupported grid type: %s", gridNamePtr(gridtype));
break;
}
}
grid->isRotated = FALSE;
if ( gribapiGetIsRotated(gh) )
{
double xpole, ypole, angle;
double xpole = 0, ypole = 0, angle = 0;
grid->isRotated = TRUE;
FAIL_ON_GRIB_ERROR(grib_get_double, gh, "latitudeOfSouthernPoleInDegrees", &ypole);
FAIL_ON_GRIB_ERROR(grib_get_double, gh, "longitudeOfSouthernPoleInDegrees", &xpole);
FAIL_ON_GRIB_ERROR(grib_get_double, gh, "angleOfRotation", &angle);
/* change from south to north pole */
if ( fabs(ypole) > 0 ) ypole = -ypole;
xpole = xpole - 180;
if ( fabs(ypole) > 0 ) ypole = -ypole; // change from south to north pole
if ( fabs(angle) > 0 ) angle = -angle;
grid->rll.xpole = xpole;
grid->rll.ypole = ypole;
......@@ -768,6 +731,7 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
}
grid->type = gridtype;
grid->projtype = projtype;
}
#endif
/*
......
......@@ -164,6 +164,7 @@ void cgribexGetGrid(stream_t *streamptr, int *isec2, double *fsec2, int *isec4,
grid_init(grid);
cdiGridTypeInit(grid, gridtype, 0);
if ( gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN || projtype == CDI_PROJ_RLL )
{
if ( ISEC4_NumValues != ISEC2_NumLon*ISEC2_NumLat )
......
......@@ -459,12 +459,24 @@ void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
// I. e. kick the fixed size array and allocate enough space, whatever that may be.
strncpy(record->varname, varname, sizeof(record->varname));
grid_t *grid = (grid_t *)Malloc(sizeof (*grid));
grid_t *grid = (grid_t *)Malloc(sizeof(*grid));
gribapiGetGrid(gh, grid);
struct addIfNewRes gridAdded = cdiVlistAddGridIfNew(vlistID, grid, 0);
int gridID = gridAdded.Id;
if (!gridAdded.isNew) Free(grid);
if ( !gridAdded.isNew ) Free(grid);
else if ( grid->projtype == CDI_PROJ_RLL )
{
double xpole = 0, ypole = 0, angle = 0;
grib_get_double(gh, "latitudeOfSouthernPoleInDegrees", &ypole);
grib_get_double(gh, "longitudeOfSouthernPoleInDegrees", &xpole);
grib_get_double(gh, "angleOfRotation", &angle);
xpole = xpole - 180;
if ( fabs(ypole) > 0 ) ypole = -ypole; // change from south to north pole
if ( fabs(angle) > 0 ) angle = -angle;
gridDefProjParamRLL(gridID, xpole, ypole, angle);
}
int zaxistype = gribapiGetZaxisType(gribEditionNumber(gh), leveltype1);
......
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