Commit 9c57b031 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Added function gridInqProjParamRLL().

parent b1aeefc2
......@@ -175,7 +175,7 @@ void cdfDefineAttributes(int cdiID, int varID, int fileID, int ncvarID)
int natts;
vlistInqNatts(cdiID, varID, &natts);
for ( int iatt = 0; iatt < natts; iatt++ )
for ( int iatt = 0; iatt < natts; ++iatt )
{
vlistInqAtt(cdiID, varID, iatt, attname, &atttype, &attlen);
......
......@@ -837,6 +837,9 @@ double gridInqYinc(int gridID);
int gridIsCircular(int gridID);
int gridIsRotated(int gridID);
void gridInqProjParamRLL(int gridID, double *xpole, double *ypole, double *angle);
void gridDefXpole(int gridID, double xpole);
double gridInqXpole(int gridID);
void gridDefYpole(int gridID, double ypole);
......
......@@ -1776,6 +1776,41 @@ double gridInqYinc(int gridID)
return yinc;
}
void gridInqProjParamRLL(int gridID, double *xpole, double *ypole, double *angle)
{
*xpole = 0; *ypole = 0; *angle = 0;
char mapping[CDI_MAX_NAME]; mapping[0] = 0;
cdiGridInqKeyStr(gridID, CDI_KEY_MAPPING, CDI_MAX_NAME, mapping);
if ( mapping[0] && strcmp(mapping, "rotated_latitude_longitude") == 0 )
{
int atttype, attlen;
char attname[CDI_MAX_NAME+1];
int natts;
vlistInqNatts(gridID, CDI_GLOBAL, &natts);
for ( int iatt = 0; iatt < natts; ++iatt )
{
vlistInqAtt(gridID, CDI_GLOBAL, iatt, attname, &atttype, &attlen);
if ( attlen != 1 ) continue;
if ( atttype == DATATYPE_FLT32 || atttype == DATATYPE_FLT64 )
{
double attflt;
vlistInqAttFlt(gridID, CDI_GLOBAL, attname, attlen, &attflt);
if ( strcmp(attname, "grid_north_pole_longitude") == 0 ) *xpole = attflt;
else if ( strcmp(attname, "grid_north_pole_latitude") == 0 ) *ypole = attflt;
else if ( strcmp(attname, "north_pole_grid_longitude") == 0 ) *angle = attflt;
}
}
}
else
Warning("rotated_latitude_longitude mapping parameter missing!");
}
/*
@Function
@Title
......
......@@ -1612,8 +1612,9 @@ void cgribexDefTime(int *isec1, int vdate, int vtime, int tsteptype, int numavg,
static
void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridID)
{
bool lcurvi = false;
static bool lwarning = true;
bool lrotatedll = false;
bool lcurvi = false;
memset(isec2, 0, 16*sizeof(int));
......@@ -1660,6 +1661,17 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI
lcurvi = true;
}
if ( gridtype == GRID_PROJECTION )
{
char mapping[CDI_MAX_NAME]; mapping[0] = 0;
cdiGridInqKeyStr(gridID, CDI_KEY_MAPPING, CDI_MAX_NAME, mapping);
if ( mapping[0] && strcmp(mapping, "rotated_latitude_longitude") == 0 )
{
gridtype = GRID_LONLAT;
lrotatedll = true;
}
}
ISEC2_Reduced = FALSE;
ISEC2_ScanFlag = 0;
......@@ -1676,7 +1688,7 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI
if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_GAUSSIAN_REDUCED )
ISEC2_GridType = GRIB1_GTYPE_GAUSSIAN;
else if ( gridtype == GRID_LONLAT && isRotated )
else if ( gridtype == GRID_LONLAT && (isRotated || lrotatedll) )
ISEC2_GridType = GRIB1_GTYPE_LATLON_ROT;
else
ISEC2_GridType = GRIB1_GTYPE_LATLON;
......@@ -1752,14 +1764,23 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI
ISEC2_ResFlag = ( ISEC2_LatIncr == 0 || ISEC2_LonIncr == 0 ) ? 0 : 128;
if ( isRotated )
{
ISEC2_LatSP = - (int)lround(gridInqYpole(gridID) * 1000);
ISEC2_LonSP = (int)lround((gridInqXpole(gridID) + 180) * 1000);
double angle = gridInqAngle(gridID);
if ( isRotated || lrotatedll )
{
double xpole = 0, ypole = 0, angle = 0;
if ( isRotated )
{
xpole = gridInqXpole(gridID);
ypole = gridInqYpole(gridID);
angle = gridInqAngle(gridID);
}
else
gridInqProjParamRLL(gridID, &xpole, &ypole, &angle);
ISEC2_LatSP = - (int)lround(ypole * 1000);
ISEC2_LonSP = (int)lround((xpole + 180) * 1000);
if ( fabs(angle) > 0 ) angle = -angle;
FSEC2_RotAngle = angle;
}
}
/* East -> West */
if ( ISEC2_LastLon < ISEC2_FirstLon ) ISEC2_ScanFlag += 128;
......
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