Skip to content

GitLab

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

Refactor GRID_LCC to GRID_PROJECTION.

parent 439cbce6
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 277 additions and 339 deletions
ChangeLog
View file @ f999f454
......@@ -2,6 +2,10 @@
* Version 1.8.1 released
2017-04-09 Uwe Schulzweida
* Refactor GRID_LCC to GRID_PROJECTION
2017-04-05 Uwe Schulzweida
* Scalar Z-Coordinate: added support for bounds (bug fix)
......
app/printinfo.h
View file @ f999f454
......@@ -279,6 +279,10 @@ void printGridInfoKernel(int gridID, int index, bool lproj)
// int prec = gridInqPrec(gridID);
// int dig = (prec == CDI_DATATYPE_FLT64) ? 15 : 7;
int dig = 7;
#ifdef CDO
extern int CDO_flt_digits;
dig = CDO_flt_digits;
#endif
if ( !lproj )
{
......@@ -406,27 +410,6 @@ void printGridInfoKernel(int gridID, int index, bool lproj)
print_xyvals2D(gridID, dig);
}
else if ( gridtype == GRID_LCC )
{
double originLon, originLat, lonParY, lat1, lat2, xincm, yincm;
int projflag, scanflag;
gridInqParamLCC(gridID, &originLon, &originLat, &lonParY, &lat1, &lat2, &xincm, &yincm,
&projflag, &scanflag);
#ifdef CDO
set_text_color(stdout, RESET, GREEN);
#endif
fprintf(stdout, "points=%d (%dx%d) ", gridsize, xsize, ysize);
if ( (projflag&128) == 0 )
fprintf(stdout, "North Pole\n");
else
fprintf(stdout, "South Pole\n");
my_reset_text_color(stdout);
fprintf(stdout, "%33s : originLon=%g originLat=%g lonParY=%g\n", " ", originLon, originLat, lonParY);
fprintf(stdout, "%33s : lat1=%g lat2=%g xinc=%g m yinc=%g m\n", " ", lat1, lat2, xincm, yincm);
}
else /* if ( gridtype == GRID_GENERIC ) */
{
#ifdef CDO
......@@ -495,6 +478,10 @@ void printZaxisInfo(int vlistID)
// int prec = zaxisInqPrec(zaxisID);
// int dig = (prec == CDI_DATATYPE_FLT64) ? 15 : 7;
int dig = 7;
#ifdef CDO
extern int CDO_flt_digits;
dig = CDO_flt_digits;
#endif
zaxisName(zaxistype, zaxisname);
zaxisInqName(zaxisID, zname);
......
src/cdi.h
View file @ f999f454
......@@ -932,8 +932,8 @@ void gridDefParamGME(int gridID, int nd, int ni, int ni2, int ni3);
void gridInqParamGME(int gridID, int *nd, int *ni, int *ni2, int *ni3);
/* Lambert Conformal Conic grid (GRIB version) */
void gridDefParamLCC(int gridID, double originLon, double originLat, double lonParY, double lat1, double lat2, double xinc, double yinc, int projflag, int scanflag);
void gridInqParamLCC(int gridID, double *originLon, double *originLat, double *lonParY, double *lat1, double *lat2, double *xinc, double *yinc, int *projflag, int *scanflag);
void gridDefParamLCC(int gridID, double a, double rf, double xval_0, double yval_0, double lon_0, double lat_1, double lat_2);
void gridInqParamLCC(int gridID, double *a, double *rf, double *xval_0, double *yval_0, double *lon_0, double *lat_1, double *lat_2);
void gridDefArea(int gridID, const double area[]);
void gridInqArea(int gridID, double area[]);
......
src/gribapi_utilities.c
View file @ f999f454
......@@ -498,6 +498,11 @@ 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 ( gridtype == GRID_LCC )
{
gridtype = GRID_PROJECTION;
projtype = CDI_PROJ_LCC;
}
/*
if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED )
{
......@@ -531,7 +536,7 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
if ( numberOfPoints != nlon*nlat )
Error("numberOfPoints (%ld) and gridSize (%ld) differ!", numberOfPoints, nlon*nlat);
grid->size = (int)numberOfPoints;
grid->size = (int)numberOfPoints;
grid->x.size = (int)nlon;
grid->y.size = (int)nlat;
grid->x.inc = 0;
......@@ -649,13 +654,12 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
grid->y.flag = 2;
}
}
else if ( gridtype == GRID_LCC )
else if ( projtype == CDI_PROJ_LCC )
{
int nlon, nlat;
FAIL_ON_GRIB_ERROR(grib_get_long, gh, "Nx", &lpar);
nlon = (int)lpar;
int nlon = (int)lpar;
FAIL_ON_GRIB_ERROR(grib_get_long, gh, "Ny", &lpar);
nlat = (int)lpar;
int nlat = (int)lpar;
if ( numberOfPoints != nlon*nlat )
Error("numberOfPoints (%d) and gridSize (%d) differ!", (int)numberOfPoints, nlon*nlat);
......@@ -664,24 +668,18 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
grid->x.size = nlon;
grid->y.size = nlat;
FAIL_ON_GRIB_ERROR(grib_get_double, gh, "DxInMetres", &grid->lcc.xinc);
FAIL_ON_GRIB_ERROR(grib_get_double, gh, "DyInMetres", &grid->lcc.yinc);
FAIL_ON_GRIB_ERROR(grib_get_double, gh, "longitudeOfFirstGridPointInDegrees", &grid->lcc.originLon);
FAIL_ON_GRIB_ERROR(grib_get_double, gh, "latitudeOfFirstGridPointInDegrees", &grid->lcc.originLat);
FAIL_ON_GRIB_ERROR(grib_get_double, gh, "LoVInDegrees", &grid->lcc.lonParY);
FAIL_ON_GRIB_ERROR(grib_get_double, gh, "Latin1InDegrees", &grid->lcc.lat1);
FAIL_ON_GRIB_ERROR(grib_get_double, gh, "Latin2InDegrees", &grid->lcc.lat2);
if ( editionNumber <= 1 )
{
FAIL_ON_GRIB_ERROR(grib_get_long, gh, "projectionCenterFlag", &lpar);
grid->lcc.projflag = (int) lpar;
FAIL_ON_GRIB_ERROR(grib_get_long, gh, "scanningMode", &lpar);
grid->lcc.scanflag = (short) lpar;
}
grid->x.flag = 0;
grid->y.flag = 0;
double xinc, yinc;
FAIL_ON_GRIB_ERROR(grib_get_double, gh, "DxInMetres", &xinc);
FAIL_ON_GRIB_ERROR(grib_get_double, gh, "DyInMetres", &yinc);
grid->x.first = 0;
grid->x.last = 0;
grid->x.inc = xinc;
grid->y.first = 0;
grid->y.last = 0;
grid->y.inc = yinc;
grid->x.flag = 2;
grid->y.flag = 2;
}
else if ( gridtype == GRID_SPECTRAL )
{
......
src/grid.c
View file @ f999f454
This diff is collapsed.
src/grid.h
View file @ f999f454
......@@ -58,22 +58,6 @@ struct gridaxis_t {
double *bounds;
};
// Lambert Conformal Conic
struct grid_lcc_t {
double originLon;
double originLat;
double lonParY;
double lat1;
double lat2;
double xinc;
double yinc;
double a;
double rf;
int projflag;
short scanflag; /* must be combination of 128, 64, 32 */
short defined;
};
// GME Grid
struct grid_gme_t {
int nd, ni, ni2, ni3; /* parameter for GRID_GME */
......@@ -93,7 +77,6 @@ struct grid_t {
mask_t *mask;
mask_t *mask_gme;
double *area;
struct grid_lcc_t lcc;
struct grid_gme_t gme;
int number, position; /* parameter for GRID_REFERENCE */
int trunc; /* parameter for GRID_SPECTEAL */
......
src/stream_cgribex.c
View file @ f999f454
......@@ -29,8 +29,6 @@ typedef struct {
#if defined (HAVE_LIBCGRIBEX)
static bool bit_get(int number, int bit) { return (bool)((number >> bit) & 1); }
static
int cgribexGetGridType(int *isec2)
{
......@@ -143,11 +141,16 @@ void cgribexGetGrid(stream_t *streamptr, int *isec2, int *isec4, grid_t *grid, i
bool compyinc = true;
int gridtype = cgribexGetGridType(isec2);
int projtype = (gridtype == GRID_PROJECTION && cgribexGetIsRotated(isec2)) ? CDI_PROJ_RLL : CDI_UNDEFID;
if ( gridtype == GRID_LCC )
{
gridtype = GRID_PROJECTION;
projtype = CDI_PROJ_LCC;
}
if ( streamptr->unreduced && gridtype == GRID_GAUSSIAN_REDUCED && iret != -801 )
{
int ilat, nlon = 0;
for ( ilat = 0; ilat < ISEC2_NumLat; ++ilat )
int nlon = 0;
for ( int ilat = 0; ilat < ISEC2_NumLat; ++ilat )
if ( ISEC2_RowLon(ilat) > nlon ) nlon = ISEC2_RowLon(ilat);
gridtype = GRID_GAUSSIAN;
ISEC2_NumLon = nlon;
......@@ -160,12 +163,14 @@ void cgribexGetGrid(stream_t *streamptr, int *isec2, int *isec4, grid_t *grid, i
if ( gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN || projtype == CDI_PROJ_RLL )
{
bool ijDirectionIncrementGiven = bit_get(ISEC2_ResFlag, 7);
bool uvRelativeToGrid = bit_get(ISEC2_ResFlag, 3);
bool ijDirectionIncrementGiven = gribbyte_get_bit(ISEC2_ResFlag, 1);
bool uvRelativeToGrid = gribbyte_get_bit(ISEC2_ResFlag, 5);
if ( uvRelativeToGrid ) grid->uvRelativeToGrid = 1;
if ( ISEC4_NumValues != ISEC2_NumLon*ISEC2_NumLat )
Error("numberOfPoints (%d) and gridSize (%d) differ!", ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);
grid->size = ISEC4_NumValues;
grid->size = ISEC4_NumValues;
grid->x.size = ISEC2_NumLon;
grid->y.size = ISEC2_NumLat;
if ( gridtype == GRID_GAUSSIAN ) grid->np = ISEC2_NumPar;
......@@ -232,8 +237,8 @@ void cgribexGetGrid(stream_t *streamptr, int *isec2, int *isec4, grid_t *grid, i
}
else if ( gridtype == GRID_GAUSSIAN_REDUCED )
{
bool ijDirectionIncrementGiven = bit_get(ISEC2_ResFlag, 7);
bool uvRelativeToGrid = bit_get(ISEC2_ResFlag, 3);
bool ijDirectionIncrementGiven = gribbyte_get_bit(ISEC2_ResFlag, 1);
bool uvRelativeToGrid = gribbyte_get_bit(ISEC2_ResFlag, 5);
if ( uvRelativeToGrid ) grid->uvRelativeToGrid = 1;
grid->np = ISEC2_NumPar;
grid->size = ISEC4_NumValues;
......@@ -273,31 +278,26 @@ void cgribexGetGrid(stream_t *streamptr, int *isec2, int *isec4, grid_t *grid, i
grid->y.flag = 2;
}
}
else if ( gridtype == GRID_LCC )
else if ( projtype == CDI_PROJ_LCC )
{
bool uvRelativeToGrid = bit_get(ISEC2_ResFlag, 3);
bool uvRelativeToGrid = gribbyte_get_bit(ISEC2_ResFlag, 5);
if ( uvRelativeToGrid ) grid->uvRelativeToGrid = 1;
if ( ISEC4_NumValues != ISEC2_NumLon*ISEC2_NumLat )
Error("numberOfPoints (%d) and gridSize (%d) differ!",
ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);
Error("numberOfPoints (%d) and gridSize (%d) differ!", ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);
grid->size = ISEC4_NumValues;
grid->x.size = ISEC2_NumLon;
grid->y.size = 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 = (short)ISEC2_ScanFlag;
grid->x.flag = 0;
grid->y.flag = 0;
grid->x.first = 0;
grid->x.last = 0;
grid->x.inc = ISEC2_Lambert_dx;
grid->y.first = 0;
grid->y.last = 0;
grid->y.inc = ISEC2_Lambert_dy;
grid->x.flag = 2;
grid->y.flag = 2;
}
else if ( gridtype == GRID_SPECTRAL )
{
......@@ -392,18 +392,17 @@ void cgribexAddRecord(stream_t * streamptr, int param, int *isec1, int *isec2, d
}
else if ( gridptr->projtype == CDI_PROJ_LCC )
{
double lcc_xinc = ISEC2_Lambert_dx;
double lcc_yinc = ISEC2_Lambert_dy;
double lcc_originLon = ISEC2_FirstLon * 0.001;
double lcc_originLat = ISEC2_FirstLat * 0.001;
double lcc_lonParY = ISEC2_Lambert_Lov * 0.001;
double lcc_lat1 = ISEC2_Lambert_LatS1 * 0.001;
double lcc_lat2 = ISEC2_Lambert_LatS2 * 0.001;
double lcc_projflag = ISEC2_Lambert_ProjFlag;
double lcc_scanflag = (short)ISEC2_ScanFlag;
gridDefParamLCC(gridID, lcc_originLon, lcc_originLat, lcc_lonParY,
lcc_lat1, lcc_lat2, lcc_xinc, lcc_yinc,
lcc_projflag, lcc_scanflag);
double a = 6367470., rf = 0;
bool earthIsOblate = gribbyte_get_bit(ISEC2_ResFlag, 2);
if ( earthIsOblate ) { a = 6378160.; rf = 297.0; }
double xval_0 = ISEC2_FirstLon * 0.001;
double yval_0 = ISEC2_FirstLat * 0.001;
double lon_0 = ISEC2_Lambert_Lov * 0.001;
double lat_1 = ISEC2_Lambert_LatS1 * 0.001;
double lat_2 = ISEC2_Lambert_LatS2 * 0.001;
bool lsouth = gribbyte_get_bit(ISEC2_Lambert_ProjFlag, 1);
if ( lsouth ) { lat_1 = -lat_1; lat_2 = -lat_2; }
gridDefParamLCC(gridID, a, rf, xval_0, yval_0, lon_0, lat_1, lat_2);
}
}
else
......@@ -1648,10 +1647,17 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI
}
}
if ( gridtype == GRID_PROJECTION && gridInqProjType(gridID) == CDI_PROJ_RLL )
if ( gridtype == GRID_PROJECTION )
{
gridtype = GRID_LONLAT;
lrotated = true;
if ( gridInqProjType(gridID) == CDI_PROJ_RLL )
{
gridtype = GRID_LONLAT;
lrotated = true;
}
else if ( gridInqProjType(gridID) == CDI_PROJ_LCC )
{
gridtype = GRID_LCC;
}
}
ISEC2_Reduced = FALSE;
......@@ -1737,8 +1743,9 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI
if ( ISEC2_NumLon == 1 && ISEC2_NumLat > 1 )
if ( ISEC2_LonIncr == 0 && ISEC2_LatIncr != 0 ) ISEC2_LonIncr = ISEC2_LatIncr;
ISEC2_ResFlag = ( ISEC2_LatIncr == 0 || ISEC2_LonIncr == 0 ) ? 0 : 128; // Set bit 7
if ( gridInqUvRelativeToGrid(gridID) ) ISEC2_ResFlag += 8; // Set bit 3
ISEC2_ResFlag = 0;
if ( ISEC2_LatIncr && ISEC2_LonIncr ) gribbyte_set_bit(&ISEC2_ResFlag, 1);
if ( gridInqUvRelativeToGrid(gridID) ) gribbyte_set_bit(&ISEC2_ResFlag, 5);
if ( lrotated )
{
......@@ -1751,43 +1758,48 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI
FSEC2_RotAngle = angle;
}
/* East -> West */
if ( ISEC2_LastLon < ISEC2_FirstLon ) ISEC2_ScanFlag += 128;
/* South -> North */
if ( ISEC2_LastLat > ISEC2_FirstLat ) ISEC2_ScanFlag += 64;
ISEC2_ScanFlag = 0;
if ( ISEC2_LastLon < ISEC2_FirstLon ) gribbyte_set_bit(&ISEC2_ScanFlag, 1); // East -> West
if ( ISEC2_LastLat > ISEC2_FirstLat ) gribbyte_set_bit(&ISEC2_ScanFlag, 2); // South -> North
break;
}
case GRID_LCC:
{
double originLon = 0.0, originLat = 0.0, lonParY = 0.0,
lat1 = 0.0, lat2 = 0.0, xincm = 0.0, yincm = 0.0;
int projflag = 0, scanflag = 0;
int xsize = gridInqXsize(gridID);
int ysize = gridInqYsize(gridID);
int xsize = gridInqXsize(gridID),
ysize = gridInqYsize(gridID);
double a = 0, rf = 0, xval_0 = 0, yval_0 = 0, lon_0 = 0, lat_1 = 0, lat_2 = 0;
gridInqParamLCC(gridID, &a, &rf, &xval_0, &yval_0, &lon_0, &lat_1, &lat_2);
bool lsouth = (lat_1 < 0);
if ( lsouth ) { lat_1 = -lat_2; lat_2 = -lat_2; }
gridInqParamLCC(gridID, &originLon, &originLat, &lonParY, &lat1, &lat2, &xincm, &yincm,
&projflag, &scanflag);
double xinc = gridInqXinc(gridID);
double yinc = gridInqYinc(gridID);
ISEC2_GridType = GRIB1_GTYPE_LCC;
ISEC2_NumLon = xsize;
ISEC2_NumLat = ysize;
ISEC2_FirstLon = (int)lround(originLon * 1000);
ISEC2_FirstLat = (int)lround(originLat * 1000);
ISEC2_Lambert_Lov = (int)lround(lonParY * 1000);
ISEC2_Lambert_LatS1 = (int)lround(lat1 * 1000);
ISEC2_Lambert_LatS2 = (int)lround(lat2 * 1000);
ISEC2_Lambert_dx = (int)lround(xincm);
ISEC2_Lambert_dy = (int)lround(yincm);
ISEC2_FirstLon = (int)lround(xval_0 * 1000);
ISEC2_FirstLat = (int)lround(yval_0 * 1000);
ISEC2_Lambert_Lov = (int)lround(lon_0 * 1000);
ISEC2_Lambert_LatS1 = (int)lround(lat_1 * 1000);
ISEC2_Lambert_LatS2 = (int)lround(lat_2 * 1000);
ISEC2_Lambert_dx = (int)lround(xinc);
ISEC2_Lambert_dy = (int)lround(yinc);
ISEC2_Lambert_LatSP = 0;
ISEC2_Lambert_LonSP = 0;
ISEC2_Lambert_ProjFlag = projflag;
ISEC2_ScanFlag = scanflag;
ISEC2_Lambert_ProjFlag = 0;
if ( lsouth ) gribbyte_set_bit(&ISEC2_Lambert_ProjFlag, 1);
bool earthIsOblate = (IS_EQUAL(a, 6378160.) && IS_EQUAL(rf, 297.));
ISEC2_ResFlag = 0;
if ( ISEC2_Lambert_dx && ISEC2_Lambert_dy ) gribbyte_set_bit(&ISEC2_ResFlag, 1);
if ( earthIsOblate ) gribbyte_set_bit(&ISEC2_ResFlag, 2);
if ( gridInqUvRelativeToGrid(gridID) ) gribbyte_set_bit(&ISEC2_ResFlag, 5);
ISEC2_ResFlag = ( ISEC2_Lambert_dx == 0 || ISEC2_Lambert_dy == 0 ) ? 0 : 128; // Set bit 7
if ( gridInqUvRelativeToGrid(gridID) ) ISEC2_ResFlag += 8; // Set bit 3
ISEC2_ScanFlag = 0;
gribbyte_set_bit(&ISEC2_ScanFlag, 2); // South -> North
break;
}
......@@ -1844,11 +1856,11 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI
{
// this will overrule/change the UvRelativeToGrid flag;
// typically when the wind is rotated with respect to north pole
bool uvRelativeToGrid = bit_get(ISEC2_ResFlag, 3);
bool uvRelativeToGrid = gribbyte_get_bit(ISEC2_ResFlag, 5);
if ( uvRelativeToGrid && !cdiGribChangeModeUvRelativeToGrid.mode )
ISEC2_ResFlag -= 8;
gribbyte_clear_bit(&ISEC2_ResFlag, 5);
else if ( !uvRelativeToGrid && cdiGribChangeModeUvRelativeToGrid.mode )
ISEC2_ResFlag += 8;
gribbyte_set_bit(&ISEC2_ResFlag, 5);
}
}
......
src/stream_grb.h
View file @ f999f454
#ifndef _STREAM_GRB_H
#define _STREAM_GRB_H
static inline bool gribbyte_get_bit(int number, int bit) { return (bool)((number >> (8-bit)) & 1); }
static inline void gribbyte_set_bit(int *number, int bit) { *number |= 1 << (8-bit); }
static inline void gribbyte_clear_bit(int *number, int bit) { *number &= ~(1 << (8-bit)); }
int grbBitsPerValue(int datatype);
int grbInqContents(stream_t *streamptr);
......
src/stream_gribapi.c
View file @ f999f454
......@@ -480,12 +480,31 @@ void gribapiAddRecord(stream_t * streamptr, int param, grib_handle *gh,
grib_get_double(gh, "latitudeOfSouthernPoleInDegrees", &ypole);
grib_get_double(gh, "longitudeOfSouthernPoleInDegrees", &xpole);
grib_get_double(gh, "angleOfRotation", &angle);
xpole = xpole - 180;
xpole -= 180;
if ( fabs(ypole) > 0 ) ypole = -ypole; // change from south to north pole
if ( fabs(angle) > 0 ) angle = -angle;
gridDefParamRLL(gridID, xpole, ypole, angle);
}
else if ( grid->projtype == CDI_PROJ_LCC )
{
double a = 6367470., rf = 0;
long earthIsOblate;
grib_get_long(gh, "earthIsOblate", &earthIsOblate);
if ( earthIsOblate ) { a = 6378160.; rf = 297.0; }
double xval_0, yval_0, lon_0, lat_1, lat_2;
long projflag = 0;
grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &xval_0);
grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &yval_0);
grib_get_double(gh, "LoVInDegrees", &lon_0);
grib_get_double(gh, "Latin1InDegrees", &lat_1);
grib_get_double(gh, "Latin2InDegrees", &lat_2);
grib_get_long(gh, "projectionCenterFlag", &projflag);
bool lsouth = gribbyte_get_bit(projflag, 1);
if ( lsouth ) { lat_1 = -lat_1; lat_2 = -lat_2; }
gridDefParamLCC(gridID, a, rf, xval_0, yval_0, lon_0, lat_1, lat_2);
}
int zaxistype = gribapiGetZaxisType(gribEditionNumber(gh), leveltype1);
......@@ -1938,10 +1957,17 @@ void gribapiDefGrid(int editionNumber, grib_handle *gh, int gridID, int comptype
}
}
if ( gridtype == GRID_PROJECTION && gridInqProjType(gridID) == CDI_PROJ_RLL )
if ( gridtype == GRID_PROJECTION )
{
gridtype = GRID_LONLAT;
lrotated = true;
if ( gridInqProjType(gridID) == CDI_PROJ_RLL )
{
gridtype = GRID_LONLAT;
lrotated = true;
}
else if ( gridInqProjType(gridID) == CDI_PROJ_LCC )
{
gridtype = GRID_LCC;
}
}
if ( gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN )
......@@ -2070,7 +2096,7 @@ void gribapiDefGrid(int editionNumber, grib_handle *gh, int gridID, int comptype
double xpole = 0, ypole = 0, angle = 0;
gridInqParamRLL(gridID, &xpole, &ypole, &angle);
xpole = xpole + 180;
xpole += 180;
if ( fabs(ypole) > 0 ) ypole = -ypole; // change from north to south pole
if ( fabs(angle) > 0 ) angle = -angle;
GRIB_CHECK(my_grib_set_double(gh, "latitudeOfSouthernPoleInDegrees", ypole), 0);
......@@ -2103,14 +2129,19 @@ void gribapiDefGrid(int editionNumber, grib_handle *gh, int gridID, int comptype
}
case GRID_LCC:
{
double originLon, originLat, lonParY, lat1, lat2, xincm, yincm;
int projflag, scanflag;
int xsize = gridInqXsize(gridID);
int ysize = gridInqYsize(gridID);
gridInqParamLCC(gridID, &originLon, &originLat, &lonParY, &lat1, &lat2, &xincm, &yincm,
&projflag, &scanflag);
double a, rf, xval_0, yval_0, lon_0, lat_1, lat_2;
gridInqParamLCC(gridID, &a, &rf, &xval_0, &yval_0, &lon_0, &lat_1, &lat_2);
if ( xval_0 < 0 ) xval_0 += 360;
bool lsouth = (lat_1 < 0);
if ( lsouth ) { lat_1 = -lat_2; lat_2 = -lat_2; }
int projflag = 0;
if ( lsouth ) gribbyte_set_bit(&projflag, 1);
double xinc = gridInqXinc(gridID);
double yinc = gridInqYinc(gridID);
static const char mesg[] = "lambert";
size_t len = sizeof(mesg) -1;
......@@ -2118,22 +2149,24 @@ void gribapiDefGrid(int editionNumber, grib_handle *gh, int gridID, int comptype
GRIB_CHECK(my_grib_set_long(gh, "Nx", xsize), 0);
GRIB_CHECK(my_grib_set_long(gh, "Ny", ysize), 0);
GRIB_CHECK(my_grib_set_long(gh, "DxInMetres", lround(xincm)), 0);
GRIB_CHECK(my_grib_set_long(gh, "DyInMetres", lround(yincm)), 0);
GRIB_CHECK(my_grib_set_double(gh, "longitudeOfFirstGridPointInDegrees", originLon), 0);
GRIB_CHECK(my_grib_set_double(gh, "latitudeOfFirstGridPointInDegrees", originLat), 0);
GRIB_CHECK(my_grib_set_double(gh, "LoVInDegrees", lonParY), 0);
GRIB_CHECK(my_grib_set_double(gh, "Latin1InDegrees", lat1), 0);
GRIB_CHECK(my_grib_set_double(gh, "Latin2InDegrees", lat2), 0);
if ( editionNumber <= 1 )
{
GRIB_CHECK(my_grib_set_long(gh, "projectionCenterFlag", projflag), 0);
GRIB_CHECK(my_grib_set_long(gh, "scanningMode", scanflag), 0);
}
GRIB_CHECK(my_grib_set_long(gh, "DxInMetres", lround(xinc)), 0);
GRIB_CHECK(my_grib_set_long(gh, "DyInMetres", lround(yinc)), 0);
GRIB_CHECK(my_grib_set_double(gh, "longitudeOfFirstGridPointInDegrees", xval_0), 0);
GRIB_CHECK(my_grib_set_double(gh, "latitudeOfFirstGridPointInDegrees", yval_0), 0);
GRIB_CHECK(my_grib_set_double(gh, "LoVInDegrees", lon_0), 0);
GRIB_CHECK(my_grib_set_double(gh, "Latin1InDegrees", lat_1), 0);
GRIB_CHECK(my_grib_set_double(gh, "Latin2InDegrees", lat_2), 0);
GRIB_CHECK(my_grib_set_long(gh, "projectionCenterFlag", projflag), 0);
long uvRelativeToGrid = gridInqUvRelativeToGrid(gridID);
if ( uvRelativeToGrid ) GRIB_CHECK(my_grib_set_long(gh, "uvRelativeToGrid", uvRelativeToGrid), 0);
long earthIsOblate = (IS_EQUAL(a, 6378160.) && IS_EQUAL(rf, 297.));
if ( earthIsOblate ) GRIB_CHECK(my_grib_set_long(gh, "earthIsOblate", earthIsOblate), 0);
int scanflag = 0;
gribbyte_set_bit(&scanflag, 2);
if ( editionNumber <= 1 )
GRIB_CHECK(my_grib_set_long(gh, "scanningMode", (long)scanflag), 0);
break;
}
......
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