Commit bcba91e4 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

ieg write: added support for projection CDI_PROJ_RLL.

parent 954d3800
2016-07-30 Uwe Schulzweida
* ieg write: added support for projection CDI_PROJ_RLL
2016-07-28 Uwe Schulzweida 2016-07-28 Uwe Schulzweida
* iegWriteVarSliceDP: does not work (bug fix) * iegWriteVarSliceDP: does not work (bug fix)
......
...@@ -142,6 +142,10 @@ extern "C" { ...@@ -142,6 +142,10 @@ extern "C" {
#define GRID_SINUSOIDAL 14 /* Sinusoidal */ #define GRID_SINUSOIDAL 14 /* Sinusoidal */
#define GRID_PROJECTION 15 /* Projected coordinates */ #define GRID_PROJECTION 15 /* Projected coordinates */
#define CDI_PROJ_RLL 11 /* Rotated Latitude Longitude */
#define CDI_PROJ_LCC 12 /* Lambert Conformal Conic */
#define CDI_PROJ_LAEA 13 /* Lambert Azimuthal Equal Area */
/* ZAXIS types */ /* ZAXIS types */
#define ZAXIS_SURFACE 0 /* Surface level */ #define ZAXIS_SURFACE 0 /* Surface level */
...@@ -723,6 +727,9 @@ int gridDuplicate(int gridID); ...@@ -723,6 +727,9 @@ int gridDuplicate(int gridID);
/* gridInqProj: Get the projection ID of a Grid */ /* gridInqProj: Get the projection ID of a Grid */
int gridInqProj(int gridID); int gridInqProj(int gridID);
/* gridInqProjType: Get the projection type */
int gridInqProjType(int gridID);
/* gridInqType: Get the type of a Grid */ /* gridInqType: Get the type of a Grid */
int gridInqType(int gridID); int gridInqType(int gridID);
......
...@@ -94,15 +94,17 @@ void grid_init(grid_t *gridptr) ...@@ -94,15 +94,17 @@ void grid_init(grid_t *gridptr)
gridptr->self = CDI_UNDEFID; gridptr->self = CDI_UNDEFID;
gridptr->type = CDI_UNDEFID; gridptr->type = CDI_UNDEFID;
gridptr->proj = CDI_UNDEFID; gridptr->proj = CDI_UNDEFID;
gridptr->projtype = CDI_UNDEFID;
gridptr->mask = NULL; gridptr->mask = NULL;
gridptr->mask_gme = NULL; gridptr->mask_gme = NULL;
gridptr->x.vals = NULL; gridptr->x.vals = NULL;
gridptr->y.vals = NULL; gridptr->y.vals = NULL;
gridptr->area = NULL;
gridptr->x.bounds = NULL; gridptr->x.bounds = NULL;
gridptr->y.bounds = NULL; gridptr->y.bounds = NULL;
gridptr->area = NULL;
gridptr->rowlon = NULL; gridptr->rowlon = NULL;
gridptr->nrowlon = 0; gridptr->nrowlon = 0;
gridptr->x.first = 0.0; gridptr->x.first = 0.0;
gridptr->x.last = 0.0; gridptr->x.last = 0.0;
gridptr->x.inc = 0.0; gridptr->x.inc = 0.0;
...@@ -979,6 +981,7 @@ void gridInqYunits(int gridID, char *yunits) ...@@ -979,6 +981,7 @@ void gridInqYunits(int gridID, char *yunits)
(void)cdiGridInqKeyStr(gridID, CDI_KEY_YUNITS, CDI_MAX_NAME, yunits); (void)cdiGridInqKeyStr(gridID, CDI_KEY_YUNITS, CDI_MAX_NAME, yunits);
} }
void gridInqYstdname(int gridID, char *ystdname) void gridInqYstdname(int gridID, char *ystdname)
{ {
grid_t *gridptr = gridID2Ptr(gridID); grid_t *gridptr = gridID2Ptr(gridID);
...@@ -988,6 +991,7 @@ void gridInqYstdname(int gridID, char *ystdname) ...@@ -988,6 +991,7 @@ void gridInqYstdname(int gridID, char *ystdname)
ystdname[0] = 0; ystdname[0] = 0;
} }
int gridInqProj(int gridID) int gridInqProj(int gridID)
{ {
grid_t *gridptr = gridID2Ptr(gridID); grid_t *gridptr = gridID2Ptr(gridID);
...@@ -995,6 +999,27 @@ int gridInqProj(int gridID) ...@@ -995,6 +999,27 @@ int gridInqProj(int gridID)
return gridptr->proj; return gridptr->proj;
} }
int gridInqProjType(int gridID)
{
grid_t *gridptr = gridID2Ptr(gridID);
int projtype = gridptr->projtype;
if ( projtype == -1 )
{
char mapping[CDI_MAX_NAME]; mapping[0] = 0;
cdiGridInqKeyStr(gridID, CDI_KEY_MAPPING, CDI_MAX_NAME, mapping);
if ( mapping[0] )
{
if ( strcmp(mapping, "rotated_latitude_longitude") == 0 ) projtype = CDI_PROJ_RLL;
gridptr->projtype = projtype;
}
}
return projtype;
}
/* /*
@Function gridInqType @Function gridInqType
@Title Get the type of a Grid @Title Get the type of a Grid
......
...@@ -105,6 +105,7 @@ struct grid_t { ...@@ -105,6 +105,7 @@ struct grid_t {
int type; /* grid type */ int type; /* grid type */
int prec; /* grid precision */ int prec; /* grid precision */
int proj; /* grid projection */ int proj; /* grid projection */
int projtype; /* grid projection type */
mask_t *mask; mask_t *mask;
mask_t *mask_gme; mask_t *mask_gme;
double *area; double *area;
......
...@@ -1613,7 +1613,7 @@ static ...@@ -1613,7 +1613,7 @@ static
void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridID) void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridID)
{ {
static bool lwarning = true; static bool lwarning = true;
bool lrotatedll = false; bool lrotated = false;
bool lcurvi = false; bool lcurvi = false;
memset(isec2, 0, 16*sizeof(int)); memset(isec2, 0, 16*sizeof(int));
...@@ -1661,15 +1661,10 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI ...@@ -1661,15 +1661,10 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI
lcurvi = true; lcurvi = true;
} }
if ( gridtype == GRID_PROJECTION ) if ( gridtype == GRID_PROJECTION && gridInqProjType(gridID) == CDI_PROJ_RLL )
{ {
char mapping[CDI_MAX_NAME]; mapping[0] = 0; gridtype = GRID_LONLAT;
cdiGridInqKeyStr(gridID, CDI_KEY_MAPPING, CDI_MAX_NAME, mapping); lrotated = true;
if ( mapping[0] && strcmp(mapping, "rotated_latitude_longitude") == 0 )
{
gridtype = GRID_LONLAT;
lrotatedll = true;
}
} }
ISEC2_Reduced = FALSE; ISEC2_Reduced = FALSE;
...@@ -1688,7 +1683,7 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI ...@@ -1688,7 +1683,7 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI
if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_GAUSSIAN_REDUCED ) if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_GAUSSIAN_REDUCED )
ISEC2_GridType = GRIB1_GTYPE_GAUSSIAN; ISEC2_GridType = GRIB1_GTYPE_GAUSSIAN;
else if ( gridtype == GRID_LONLAT && (isRotated || lrotatedll) ) else if ( gridtype == GRID_LONLAT && (isRotated || lrotated) )
ISEC2_GridType = GRIB1_GTYPE_LATLON_ROT; ISEC2_GridType = GRIB1_GTYPE_LATLON_ROT;
else else
ISEC2_GridType = GRIB1_GTYPE_LATLON; ISEC2_GridType = GRIB1_GTYPE_LATLON;
...@@ -1764,7 +1759,7 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI ...@@ -1764,7 +1759,7 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI
ISEC2_ResFlag = ( ISEC2_LatIncr == 0 || ISEC2_LonIncr == 0 ) ? 0 : 128; ISEC2_ResFlag = ( ISEC2_LatIncr == 0 || ISEC2_LonIncr == 0 ) ? 0 : 128;
if ( isRotated || lrotatedll ) if ( isRotated || lrotated )
{ {
double xpole = 0, ypole = 0, angle = 0; double xpole = 0, ypole = 0, angle = 0;
if ( isRotated ) if ( isRotated )
......
...@@ -246,13 +246,19 @@ calc_resfac(double xfirst, double xlast, double xinc, double yfirst, double ylas ...@@ -246,13 +246,19 @@ calc_resfac(double xfirst, double xlast, double xinc, double yfirst, double ylas
static static
void iegDefGrid(int *gdb, int gridID) void iegDefGrid(int *gdb, int gridID)
{ {
int projID = gridInqProj(gridID);
if ( projID != CDI_UNDEFID && gridInqProjType(projID) == CDI_PROJ_RLL ) gridID = projID;
int gridtype = gridInqType(gridID); int gridtype = gridInqType(gridID);
int projtype = CDI_UNDEFID;
if ( gridtype == GRID_PROJECTION && gridInqProjType(gridID) == CDI_PROJ_RLL ) projtype = CDI_PROJ_RLL;
int xsize = gridInqXsize(gridID);
int ysize = gridInqYsize(gridID);
if ( gridtype == GRID_GENERIC ) if ( gridtype == GRID_GENERIC )
{ {
int xsize = gridInqXsize(gridID);
int ysize = gridInqYsize(gridID);
if ( (ysize == 32 || ysize == 48 || ysize == 64 || if ( (ysize == 32 || ysize == 48 || ysize == 64 ||
ysize == 96 || ysize == 160) && ysize == 96 || ysize == 160) &&
(xsize == 2*ysize || xsize == 1) ) (xsize == 2*ysize || xsize == 1) )
...@@ -276,39 +282,32 @@ void iegDefGrid(int *gdb, int gridID) ...@@ -276,39 +282,32 @@ void iegDefGrid(int *gdb, int gridID)
gridtype = GRID_LONLAT; gridtype = GRID_LONLAT;
} }
if ( gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN ) bool lrotated = ((gridtype == GRID_LONLAT && gridIsRotated(gridID)) || projtype == CDI_PROJ_RLL);
if ( gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN || projtype == CDI_PROJ_RLL )
{ {
double xfirst = 0, xlast = 0, xinc = 0; double xfirst = 0, xlast = 0, xinc = 0;
double yfirst = 0, ylast = 0, yinc = 0; double yfirst = 0, ylast = 0, yinc = 0;
int nlon = gridInqXsize(gridID), if ( xsize == 0 ) xsize = 1;
nlat = gridInqYsize(gridID);
if ( nlon == 0 )
{
nlon = 1;
}
else else
{ {
xfirst = gridInqXval(gridID, 0); xfirst = gridInqXval(gridID, 0);
xlast = gridInqXval(gridID, nlon-1); xlast = gridInqXval(gridID, xsize-1);
xinc = gridInqXinc(gridID); xinc = gridInqXinc(gridID);
} }
if ( nlat == 0 ) if ( ysize == 0 ) ysize = 1;
{
nlat = 1;
}
else else
{ {
yfirst = gridInqYval(gridID, 0); yfirst = gridInqYval(gridID, 0);
ylast = gridInqYval(gridID, nlat-1); ylast = gridInqYval(gridID, ysize-1);
yinc = gridInqYinc(gridID); yinc = gridInqYinc(gridID);
} }
if ( gridtype == GRID_GAUSSIAN ) if ( gridtype == GRID_GAUSSIAN )
IEG_G_GridType(gdb) = 4; IEG_G_GridType(gdb) = 4;
else if ( gridtype == GRID_LONLAT && gridIsRotated(gridID) ) else if ( lrotated )
IEG_G_GridType(gdb) = 10; IEG_G_GridType(gdb) = 10;
else else
IEG_G_GridType(gdb) = 0; IEG_G_GridType(gdb) = 0;
...@@ -319,8 +318,8 @@ void iegDefGrid(int *gdb, int gridID) ...@@ -319,8 +318,8 @@ void iegDefGrid(int *gdb, int gridID)
IEG_G_ResFac(gdb) = iresfac; IEG_G_ResFac(gdb) = iresfac;
IEG_G_NumLon(gdb) = nlon; IEG_G_NumLon(gdb) = xsize;
IEG_G_NumLat(gdb) = nlat; IEG_G_NumLat(gdb) = ysize;
IEG_G_FirstLat(gdb) = (int)lround(yfirst*resfac); IEG_G_FirstLat(gdb) = (int)lround(yfirst*resfac);
IEG_G_LastLat(gdb) = (int)lround(ylast*resfac); IEG_G_LastLat(gdb) = (int)lround(ylast*resfac);
IEG_G_FirstLon(gdb) = (int)lround(xfirst*resfac); IEG_G_FirstLon(gdb) = (int)lround(xfirst*resfac);
...@@ -330,7 +329,7 @@ void iegDefGrid(int *gdb, int gridID) ...@@ -330,7 +329,7 @@ void iegDefGrid(int *gdb, int gridID)
IEG_G_LonIncr(gdb) = 0; IEG_G_LonIncr(gdb) = 0;
if ( gridtype == GRID_GAUSSIAN ) if ( gridtype == GRID_GAUSSIAN )
IEG_G_LatIncr(gdb) = nlat/2; IEG_G_LatIncr(gdb) = ysize/2;
else else
{ {
IEG_G_LatIncr(gdb) = (int)lround(yinc*resfac); IEG_G_LatIncr(gdb) = (int)lround(yinc*resfac);
...@@ -346,15 +345,22 @@ void iegDefGrid(int *gdb, int gridID) ...@@ -346,15 +345,22 @@ void iegDefGrid(int *gdb, int gridID)
if ( IEG_G_NumLon(gdb) == 1 && IEG_G_NumLat(gdb) > 1 ) if ( IEG_G_NumLon(gdb) == 1 && IEG_G_NumLat(gdb) > 1 )
if ( IEG_G_LonIncr(gdb) == 0 && IEG_G_LatIncr(gdb) != 0 ) IEG_G_LonIncr(gdb) = IEG_G_LatIncr(gdb); if ( IEG_G_LonIncr(gdb) == 0 && IEG_G_LatIncr(gdb) != 0 ) IEG_G_LonIncr(gdb) = IEG_G_LatIncr(gdb);
if ( IEG_G_LatIncr(gdb) == 0 || IEG_G_LonIncr(gdb) == 0 ) IEG_G_ResFlag(gdb) = (IEG_G_LatIncr(gdb) == 0 || IEG_G_LonIncr(gdb) == 0) ? 0 : 128;
IEG_G_ResFlag(gdb) = 0;
else
IEG_G_ResFlag(gdb) = 128;
if ( gridIsRotated(gridID) ) if ( lrotated )
{ {
IEG_G_LatSP(gdb) = - (int)lround(gridInqYpole(gridID) * resfac); double xpole = 0, ypole = 0, angle = 0;
IEG_G_LonSP(gdb) = (int)lround((gridInqXpole(gridID) + 180) * resfac); if ( projtype == CDI_PROJ_RLL )
gridInqProjParamRLL(gridID, &xpole, &ypole, &angle);
else
{
xpole = gridInqXpole(gridID);
ypole = gridInqYpole(gridID);
angle = gridInqAngle(gridID);
}
IEG_G_LatSP(gdb) = - (int)lround(ypole * resfac);
IEG_G_LonSP(gdb) = (int)lround((xpole + 180) * resfac);
IEG_G_Size(gdb) = 42; IEG_G_Size(gdb) = 42;
} }
else else
......
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