Commit 22f46e00 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

ieg read: added support for projection CDI_PROJ_RLL.

parent f903c72d
2016-08-01 Uwe Schulzweida
* ieg read: added support for projection CDI_PROJ_RLL
2016-07-30 Uwe Schulzweida
* ieg write: added support for projection CDI_PROJ_RLL
......
......@@ -846,6 +846,7 @@ int gridIsCircular(int gridID);
int gridIsRotated(int gridID);
void gridInqProjParamRLL(int gridID, double *xpole, double *ypole, double *angle);
void gridDefProjParamRLL(int gridID, double xpole, double ypole, double angle);
void gridDefXpole(int gridID, double xpole);
double gridInqXpole(int gridID);
......
......@@ -1826,6 +1826,15 @@ void gridInqProjParamRLL(int gridID, double *xpole, double *ypole, double *angle
Warning("rotated_latitude_longitude mapping parameter missing!");
}
void gridDefProjParamRLL(int gridID, double xpole, double ypole, double angle)
{
cdiGridDefKeyStr(gridID, CDI_KEY_MAPPING, CDI_MAX_NAME, "rotated_latitude_longitude");
vlistDefAttFlt(gridID, CDI_GLOBAL, "grid_north_pole_longitude", DATATYPE_FLT64, 1, &xpole);
vlistDefAttFlt(gridID, CDI_GLOBAL, "grid_north_pole_latitude", DATATYPE_FLT64, 1, &ypole);
if ( IS_NOT_EQUAL(angle, 0) ) vlistDefAttFlt(gridID, CDI_GLOBAL, "north_pole_grid_longitude", DATATYPE_FLT64, 1, &angle);
}
/*
@Function
@Title
......
......@@ -588,9 +588,14 @@ void iegAddRecord(stream_t *streamptr, int param, int *pdb, int *gdb, double *vc
record->ilevel2 = level2;
record->ltype = IEG_P_LevelType(pdb);
int gridtype =
( IEG_G_GridType(gdb) == 0 || IEG_G_GridType(gdb) == 10 ) ? GRID_LONLAT :
( IEG_G_GridType(gdb) == 4 ) ? GRID_GAUSSIAN : GRID_GENERIC;
#ifdef TEST_PROJECTION
int gridtype = (IEG_G_GridType(gdb) == 0) ? GRID_LONLAT :
(IEG_G_GridType(gdb) == 10) ? GRID_PROJECTION :
(IEG_G_GridType(gdb) == 4) ? GRID_GAUSSIAN : GRID_GENERIC;
#else
int gridtype = (IEG_G_GridType(gdb) == 0 || IEG_G_GridType(gdb) == 10) ? GRID_LONLAT :
(IEG_G_GridType(gdb) == 4) ? GRID_GAUSSIAN : GRID_GENERIC;
#endif
grid_t *grid = (grid_t *)Malloc(sizeof (*grid));
grid_init(grid);
......@@ -648,27 +653,33 @@ void iegAddRecord(stream_t *streamptr, int param, int *pdb, int *gdb, double *vc
}
grid->isRotated = FALSE;
double xpole = 0, ypole = 0;
if ( IEG_G_GridType(gdb) == 10 )
{
grid->isRotated = TRUE;
grid->rll.ypole = - IEG_G_LatSP(gdb) * resfac;
grid->rll.xpole = IEG_G_LonSP(gdb) * resfac - 180;
grid->rll.angle = 0;
xpole = IEG_G_LonSP(gdb) * resfac - 180;
ypole = - IEG_G_LatSP(gdb) * resfac;
if ( gridtype == GRID_LONLAT )
{
grid->isRotated = TRUE;
grid->rll.xpole = xpole;
grid->rll.ypole = ypole;
grid->rll.angle = 0;
}
}
struct addIfNewRes gridAdded = cdiVlistAddGridIfNew(vlistID, grid, 0);
int gridID = gridAdded.Id;
if (!gridAdded.isNew) Free(grid);
if ( !gridAdded.isNew ) Free(grid);
else if ( gridtype == GRID_PROJECTION ) gridDefProjParamRLL(gridID, xpole, ypole, 0);
int leveltype = iegGetZaxisType(IEG_P_LevelType(pdb));
if ( leveltype == ZAXIS_HYBRID )
{
double tmpvct[100];
size_t vctsize = (size_t)IEG_G_NumVCP(gdb);
for (size_t i = 0; i < vctsize/2; i++ ) tmpvct[i] = vct[i];
for (size_t i = 0; i < vctsize/2; i++ ) tmpvct[i+vctsize/2] = vct[i+50];
for ( size_t i = 0; i < vctsize/2; i++ ) tmpvct[i] = vct[i];
for (size_t i = 0; i < vctsize/2; i++ ) tmpvct[i+vctsize/2] = vct[i+50];
varDefVCT(vctsize, tmpvct);
}
......
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