Skip to content

GitLab

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

Merged branch Hirlam_extensions.

parent 7ae445c7
Hide whitespace changes
Inline Side-by-side
src/cdi.h
View file @ cd2c0bac
......@@ -967,6 +967,9 @@ int gridInqComplexPacking(int gridID);
void gridDefUvRelativeToGrid(int gridID, int uvRelativeToGrid);
int gridInqUvRelativeToGrid(int gridID);
void gridDefScanningMode(int gridID, int mode);
int gridInqScanningMode(int gridID);
/* ZAXIS routines */
void zaxisName(int zaxistype, char *zaxisname);
......
src/cdi_int.h
View file @ cd2c0bac
......@@ -152,6 +152,14 @@ typedef struct
int ilevel2;
int ltype;
short tsteptype;
#ifdef HIRLAM_EXTENSIONS
// NOTE: tsteptype MUST be part of attributes used to compare variables!
// Modern NWP models (HARMONIE, HIRLAM) use timeRangeIndicator to specify
// if the field is instantanous or accumulated.
// Both types are typically in the same GRIB-file.
// (181; 105, 0, timeRangeIndicator=0) .. instantanous rain
// (181; 105, 0, timeRangeIndicator=4) .. accumulated rain .. both can be in the same grib file
#endif // HIRLAM_EXTENSIONS
short used;
short varID;
short levelID;
......@@ -303,6 +311,7 @@ typedef enum {
#define CDI_FILETYPE_UNDEF -1 /* Unknown/not yet defined file type */
extern int cdiDebugExt;
extern int CDI_Debug; /* If set to 1, debuggig (default 0) */
extern int CDI_Recopt;
extern int cdiGribApiDebug;
......
src/grb_write.c
View file @ cd2c0bac
......@@ -133,6 +133,71 @@ void grbCopyRecord(stream_t * streamptr2, stream_t * streamptr1)
}
}
#ifdef HIRLAM_EXTENSIONS
// Even if you are stream-copy records you might need to change a bit of grib-header !
if ( cdiGribDataScanningMode.active )
// allowed modes: <0, 64, 96>; Default is 64
// This will overrule the old scanning mode of the given grid
{
grib_handle *gh = NULL;
gh = grib_handle_new_from_message(NULL, (void *) gribbuffer, recsize);
int scanModeIN = gribapiGetScanningMode(gh);
grib_handle_delete(gh);
if (cdiDebugExt>=20) Message("Change GribDataScanningMode => %d (scanModeIN = %d)", cdiGribDataScanningMode.value, scanModeIN);
if (scanModeIN != cdiGribDataScanningMode.value)
{
int gridID;
int varID, levelID;
int vlistID;
//int zip;
int gridsize, nmiss = 0;
vlistID = streamptr1->vlistID;
varID = streamptr1->tsteps[tsID].records[recID].varID;
levelID = streamptr1->tsteps[tsID].records[recID].levelID;
//gribbuffer = (unsigned char *) streamptr->record->buffer;
// allocate above ..
gridID = vlistInqVarGrid(vlistID, varID);
gridsize = gridInqSize(gridID);
gridsize = vlistGridsizeMax(vlistID);
if ( vlistNumber(vlistID) != CDI_REAL ) gridsize *= 2;
double * data = (double *) malloc(gridsize*sizeof(double));
//int missval = vlistInqVarMissval(vlistID, varID);
//streamptr->numvals += gridsize;
// memtype: MEMTYPE_FLOAT or MEMTYPE_DOUBLE
//int statusDC = grbDecode(filetype, MEMTYPE_DOUBLE, gribbuffer, recsize, data, gridsize, streamptr1->unreduced, &nmiss, missval, vlistID, varID);
//int grbDecode(int filetype, int memtype, void *gribbuffer, int gribsize, void *data, size_t datasize,
// int unreduced, int *nmiss, double missval, int vlistID, int varID);
//streamptr1->tsteps[tsID].records[recID].zip = zip;
//gribapiSetScanningMode(gh, cdoGribDataScanningMode); // T.B.D. this will be done by grbDecode..
//varID = streamptr1->record->varID;
//levelID = streamptr1->record->levelID;
if (cdiDebugExt>=20) Message(" processing varID %d; levelID %d",varID,levelID);
grb_write_var_slice(streamptr2, varID, levelID, MEMTYPE_DOUBLE, (const void *) data, nmiss);
//grb_write_var_slice(streamptr, varID, levelID, memtype, ((double*)data)+levelID*gridsize, nmiss);
//grb_write_var(streamptr2, varID, MEMTYPE_DOUBLE, data, nmiss);
//grb_write_var(stream_t *streamptr, int varID, int memtype, const void *data, int nmiss)
//grb_write_var_slice(streamptr2, varID, levelID, MEMTYPE_DOUBLE, (const void *) data, nmiss);
free(data);
free(gribbuffer);
}
}
#endif // HIRLAM_EXTENSIONS
if ( filetype == CDI_FILETYPE_GRB )
{
size_t unzipsize;
......
src/gribapi_utilities.c
View file @ cd2c0bac
......@@ -382,6 +382,48 @@ int gribapiGetTsteptype(grib_handle *gh)
// printf("stepType: %s %ld %d\n", stepType, len, tsteptype);
}
#ifdef HIRLAM_EXTENSIONS
{
// Normaly cdo looks in grib for attribute called "stepType", see above.
// BUT NWP models such as Hirlam and Harmonie 37h1.2, use "timeRangeIndicator" instead!
// Where for example: 0: for instanteneous fields; 4: for accumulated fields
// 0: Forecast product valid at reference time + P1
// 2: Product with a valid time ranging between reference time + P1 and reference time + P2
// 4: Accumulation (reference time + P1 to reference time + P2)
// 5: Difference(reference time + P2 minus reference time + P1) product considered valid at reference time + P2
// More details on WMO standards:
// http://www.wmo.int/pages/prog/www/WDM/Guides/Guide-binary-2.html
//tsteptype = TSTEP_INSTANT; // default value for any case
long timeRangeIND = 0; // typically 0: for instanteneous fields; 4: for accumulated fields
int rc = grib_get_long(gh, "timeRangeIndicator", &timeRangeIND);
if (rc != 0) {
//if ( lprint )
Warning("Could not get 'stepType' either 'timeRangeIndicator'. Using defualt!");
return (tsteptype);
}
extern int cdiGribUseTimeRangeIndicator;
cdiGribUseTimeRangeIndicator = 1;
switch ( timeRangeIND )
{
case 0: tsteptype = TSTEP_INSTANT; break;
case 2: tsteptype = TSTEP_INSTANT2;
strcpy(stepType, "instant2"); // was incorrectly set before into accum
break;
case 4: tsteptype = TSTEP_ACCUM; break;
case 5: tsteptype = TSTEP_DIFF; break;
default:
if ( lprint )
{
if (CDI_Debug)
Warning("timeRangeIND = %d; stepType= %s; tsteptype=%d unsupported timeRangeIND at the moment, set to instant!", timeRangeIND, stepType, tsteptype);
lprint = FALSE;
}
break;
}
if (CDI_Debug)
Warning("timeRangeIND = %d; stepType= %s; tsteptype=%d", timeRangeIND, stepType, tsteptype);
}
#endif // HIRLAM_EXTENSIONS
}
return tsteptype;
......@@ -723,6 +765,34 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
grid->uvRelativeToGrid = (bool)temp;
}
if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_LONLAT || gridtype == GRID_PROJECTION || gridtype == GRID_LCC )
{
GRIB_CHECK(grib_get_long(gh, "iScansNegatively", &grid->iScansNegatively), 0);
GRIB_CHECK(grib_get_long(gh, "jScansPositively", &grid->jScansPositively), 0);
GRIB_CHECK(grib_get_long(gh, "jPointsAreConsecutive", &grid->jPointsAreConsecutive), 0);
grid->scanningMode = 128*grid->iScansNegatively + 64*grid->jScansPositively + 32*grid->jPointsAreConsecutive;
/* scanningMode = 128 * iScansNegatively + 64 * jScansPositively + 32 * jPointsAreConsecutive;
64 = 128 * 0 + 64 * 1 + 32 * 0
00 = 128 * 0 + 64 * 0 + 32 * 0
96 = 128 * 0 + 64 * 1 + 32 * 1
Default / implicit scanning mode is 64:
i and j scan positively, i points are consecutive (row-major) */
#ifdef HIRLAM_EXTENSIONS
if (cdiDebugExt>=30)
{
// indicatorOfParameter=33,indicatorOfTypeOfLevel=105,level
long paramId, levelTypeId, levelId;
GRIB_CHECK(grib_get_long(gh, "indicatorOfParameter", &paramId), 0);
GRIB_CHECK(grib_get_long(gh, "indicatorOfTypeOfLevel", &levelTypeId), 0);
GRIB_CHECK(grib_get_long(gh, "level", &levelId), 0);
Message("(param,ltype,level) = (%3d,%3d,%4d); Scanning mode = %02d -> bits:(%1d.%1d.%1d)*32; uvRelativeToGrid = %02d",\
(int)paramId, (int)levelTypeId, (int)levelId,
grid->scanningMode,(int)grid->jPointsAreConsecutive,(int)grid->jScansPositively,(int)grid->iScansNegatively,
(int)grid->uvRelativeToGrid);
}
#endif //HIRLAM_EXTENSIONS
}
grid->type = gridtype;
grid->projtype = projtype;
}
......
src/gribapi_utilities.h
View file @ cd2c0bac
......@@ -33,6 +33,11 @@ int gribapiGetParam(grib_handle *gh);
int gribapiGetGridType(grib_handle *gh);
void gribapiGetGrid(grib_handle *gh, grid_t *grid);
#ifdef HIRLAM_EXTENSIONS
void gribapiSetDataTimeRangeIndicator(grib_handle *gh, int timeRangeIndicator);
void gribapiGetDataTimeRangeIndicator(grib_handle *gh, int *timeRangeIndicator);
#endif // #ifdef HIRLAM_EXTENSIONS
#endif
#endif
......
src/grid.c
View file @ cd2c0bac
......@@ -158,6 +158,16 @@ void grid_init(grid_t *gridptr)
gridptr->atts.nalloc = MAX_ATTRIBUTES;
gridptr->atts.nelems = 0;
gridptr->uvRelativeToGrid = 0; // Some models deliver wind U,V relative to the grid-cell
gridptr->iScansNegatively = 0;
gridptr->jScansPositively = 1;
gridptr->jPointsAreConsecutive = 0;
gridptr->scanningMode = 128*gridptr->iScansNegatively + 64*gridptr->jScansPositively + 32*gridptr->jPointsAreConsecutive;
/* scanningMode = 128 * iScansNegatively + 64 * jScansPositively + 32 * jPointsAreConsecutive;
64 = 128 * 0 + 64 * 1 + 32 * 0
00 = 128 * 0 + 64 * 0 + 32 * 0
96 = 128 * 0 + 64 * 1 + 32 * 1
Default / implicit scanning mode is 64:
i and j scan positively, i points are consecutive (row-major) */
}
......@@ -2227,10 +2237,15 @@ bool gridCompare(int gridID, const grid_t *grid, bool coord_compare)
}
}
if ( (grid->uvRelativeToGrid != gridInqUvRelativeToGrid(gridID)) )
if ( (grid->scanningMode != gridInqScanningMode(gridID)) || (grid->uvRelativeToGrid != gridInqUvRelativeToGrid(gridID)) )
{
// often grid definition may differ in UV-relativeToGrid
differ = 1;
#ifdef HIRLAM_EXTENSIONS
if ( cdiDebugExt>=200 )
printf("gridCompare(gridID=%d): Differs: grid.scanningMode [%d] != gridInqScanningMode(gridID) [%d] or grid.uvRelativeToGrid [%ld] != gridInqUvRelativeToGrid(gridID) [%d]\n",
gridID, grid->scanningMode, gridInqScanningMode(gridID), grid->uvRelativeToGrid, gridInqUvRelativeToGrid(gridID) );
#endif // HIRLAM_EXTENSIONS
}
return differ;
}
......@@ -2293,6 +2308,7 @@ int gridCompareP(void *gridptr1, void *gridptr2)
if ( IS_NOT_EQUAL(g1->lcc.xinc , g2->lcc.xinc) ) return differ;
if ( IS_NOT_EQUAL(g1->lcc.yinc , g2->lcc.yinc) ) return differ;
if ( IS_NOT_EQUAL(g1->uvRelativeToGrid , g2->uvRelativeToGrid) ) return differ;
if ( IS_NOT_EQUAL(g1->scanningMode , g2->scanningMode) ) return differ;
const double *restrict g1_xvals = g1->vtable->inqXValsPtr(g1),
*restrict g2_xvals = g2->vtable->inqXValsPtr(g2);
......@@ -2632,6 +2648,10 @@ int gridGenerate(const grid_t *grid)
gridptr->number = grid->number;
gridptr->position = grid->position;
gridptr->uvRelativeToGrid = grid->uvRelativeToGrid;
gridptr->scanningMode = grid->scanningMode;
gridptr->iScansNegatively = grid->iScansNegatively;
gridptr->jScansPositively = grid->jScansPositively;
gridptr->jPointsAreConsecutive = grid->jPointsAreConsecutive;
memcpy(gridptr->uuid, grid->uuid, CDI_UUID_SIZE);
if ( gridtype == GRID_UNSTRUCTURED && grid->reference )
gridDefReference(gridID, grid->reference);
......@@ -3665,7 +3685,7 @@ void gridDefParamLCC(int gridID, double originLon, double originLat, double lonP
@Item yinc Y-direction grid lenght in meter.
@Item projflag Projection centre flag.
@Item scanflag Scanning mode flag.
@Description
The function @func{gridInqParamLCC} returns the parameter of a Lambert Conformal Conic grid.
......@@ -3962,6 +3982,31 @@ int gridInqUvRelativeToGrid(int gridID)
}
void gridDefScanningMode(int gridID, int mode)
{
grid_t *gridptr = grid_to_pointer(gridID);
if ( gridptr->scanningMode != mode )
{
gridMark4Update(gridID);
gridptr->scanningMode = mode;
}
}
int gridInqScanningMode(int gridID)
{
grid_t *gridptr = grid_to_pointer(gridID);
int scanningModeTMP = 128 * gridptr->iScansNegatively + 64 * gridptr->jScansPositively + 32 * gridptr->jPointsAreConsecutive;
if ( scanningModeTMP != gridptr->scanningMode )
Message("WARNING: scanningMode (%d) ! = (%d) 128 * iScansNegatively(%d) + 64 * jScansPositively(%d) + 32 * jPointsAreConsecutive(%d) ",
gridptr->scanningMode, scanningModeTMP, gridptr->iScansNegatively,gridptr->jScansPositively,gridptr->jPointsAreConsecutive );
return gridptr->scanningMode;
}
void cdiGridGetIndexList(unsigned ngrids, int * gridIndexList)
{
reshGetResHListOfType(ngrids, gridIndexList, &gridOps);
......@@ -4000,6 +4045,11 @@ enum {
GRID_PACK_INT_IDX_MEMBERMASK,
GRID_PACK_INT_IDX_XTSTDNNAME,
GRID_PACK_INT_IDX_YTSTDNNAME,
GRID_PACK_INT_IDX_UVRELATIVETOGRID,
GRID_PACK_INT_IDX_ISCANSNEGATIVELY,
GRID_PACK_INT_IDX_JSCANSPOSITIVELY,
GRID_PACK_INT_IDX_JPOINTSARECONSECUTIVE,
GRID_PACK_INT_IDX_SCANNINGMODE,
gridNint
};
......@@ -4213,6 +4263,11 @@ gridUnpack(char * unpackBuffer, int unpackBufferSize,
xystdname_tab[intBuffer[GRID_PACK_INT_IDX_XTSTDNNAME]][0];
gridP->y.stdname =
xystdname_tab[intBuffer[GRID_PACK_INT_IDX_YTSTDNNAME]][1];
gridP->uvRelativeToGrid = intBuffer[GRID_PACK_INT_IDX_UVRELATIVETOGRID];
gridP->iScansNegatively = intBuffer[GRID_PACK_INT_IDX_ISCANSNEGATIVELY];
gridP->jScansPositively = intBuffer[GRID_PACK_INT_IDX_JSCANSPOSITIVELY];
gridP->jPointsAreConsecutive = intBuffer[GRID_PACK_INT_IDX_JPOINTSARECONSECUTIVE];
gridP->scanningMode = intBuffer[GRID_PACK_INT_IDX_SCANNINGMODE];
}
if (memberMask & gridHasRowLonFlag)
......@@ -4407,6 +4462,12 @@ gridPack(void * voidP, void * packBuffer, int packBufferSize,
(int)((const char (*)[2][24])gridP->y.stdname
- (const char (*)[2][24])xystdname_tab[0][1]);
intBuffer[GRID_PACK_INT_IDX_UVRELATIVETOGRID] = gridP->uvRelativeToGrid;
intBuffer[GRID_PACK_INT_IDX_ISCANSNEGATIVELY] = gridP->iScansNegatively;
intBuffer[GRID_PACK_INT_IDX_JSCANSPOSITIVELY] = gridP->jScansPositively;
intBuffer[GRID_PACK_INT_IDX_JPOINTSARECONSECUTIVE] = gridP->jPointsAreConsecutive;
intBuffer[GRID_PACK_INT_IDX_SCANNINGMODE] = gridP->scanningMode;
serializePack(intBuffer, gridNint, CDI_DATATYPE_INT,
packBuffer, packBufferSize, packBufferPos, context);
d = cdiCheckSum(CDI_DATATYPE_INT, gridNint, intBuffer);
......
src/grid.h
View file @ cd2c0bac
......@@ -108,12 +108,21 @@ struct grid_t {
*/
bool lcomplex;
bool hasdims;
/* Some models deliver wind U,V relative to the grid-cell */
bool uvRelativeToGrid;
bool uvRelativeToGrid; /* Some models deliver wind U,V relative to the grid-cell */
struct gridaxis_t x;
struct gridaxis_t y;
const struct gridVirtTable *vtable;
cdi_atts_t atts;
long iScansNegatively;
long jScansPositively;
long jPointsAreConsecutive;
int scanningMode;
/* scanningMode = 128 * iScansNegatively + 64 * jScansPositively + 32 * jPointsAreConsecutive;
64 = 128 * 0 + 64 * 1 + 32 * 0
00 = 128 * 0 + 64 * 0 + 32 * 0
96 = 128 * 0 + 64 * 1 + 32 * 1
Default / implicit scanning mode is 64:
i and j scan positively, i points are consecutive (row-major) */
};
......
src/stream_grb.c
View file @ cd2c0bac
......@@ -10,6 +10,14 @@
#include "file.h"
#include "cgribex.h" /* gribZip gribGetZip gribGinfo */
int cdiDebugExt = 0; // Debug level for the KNMI extensions
#ifdef HIRLAM_EXTENSIONS
// *** RELATED to GRIB only ***
int cdiGribUseTimeRangeIndicator = 0; // normaly cdo looks in grib for attribute called "stepType"
// but NWP models such as Harmonie 37h1.2, use "timeRangeIndicator"
// where: 0: for instanteneous fields; 4: for accumulated fields
#endif // HIRLAM_EXTENSIONS
// Regarding operation to change parameter identification:
// change if cdiGribChangeParameterID.active
......@@ -43,12 +51,12 @@ struct cdiGribScanModeChange cdiGribDataScanningMode;
void streamGrbDefDataScanningMode(int scanmode)
{
cdiGribDataScanningMode.active = true;
cdiGribDataScanningMode.scanmode = scanmode;
cdiGribDataScanningMode.value = scanmode;
}
int streamGrbInqDataScanningMode(void)
{
return cdiGribDataScanningMode.scanmode;
return cdiGribDataScanningMode.value;
}
......
src/stream_grb.h
View file @ cd2c0bac
......@@ -38,7 +38,7 @@ struct cdiGribModeChange
struct cdiGribScanModeChange
{
int scanmode;
int value;
bool active;
};
......
src/stream_gribapi.c
View file @ cd2c0bac
......@@ -36,6 +36,14 @@ typedef struct {
int level2;
int ltype;
int tsteptype;
#ifdef HIRLAM_EXTENSIONS
// NOTE: tsteptype MUST be part of attributes used to compare variables!
// Modern NWP models (HARMONIE, HIRLAM) use timeRangeIndicator to specify
// if the field is instantanous or accumulated.
// Both types are typically in the same GRIB-file.
// (181; 105, 0, timeRangeIndicator=0) .. instantanous rain
// (181; 105, 0, timeRangeIndicator=4) .. accumulated rain .. both can be in the same grib file
#endif // HIRLAM_EXTENSIONS
char name[32];
var_tile_t tiles;
......@@ -43,6 +51,9 @@ typedef struct {
} compvar2_t;
int scanModeIN; // global variable
static
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
{
......@@ -1670,6 +1681,8 @@ void gribapiDefStepUnits(int editionNumber, grib_handle *gh, int timeunit, int p
}
else
{
// NOTE KNMI: HIRLAM model files LAMH_D11 are in grib1 and do NOT have key indicatorOfUnitForTimeRange
// Watch out for compatibility issues.
GRIB_CHECK(my_grib_set_long(gh, "indicatorOfUnitOfTimeRange", unitsOfTime), 0);
}
}
......@@ -2558,6 +2571,358 @@ void gribapiDefLevel(int editionNumber, grib_handle *gh, int zaxisID, int levelI
}
}
int gribapiGetScanningMode(grib_handle *gh)
{
long iScansNegatively;
long jScansPositively;
long jPointsAreConsecutive;
GRIB_CHECK(grib_get_long(gh, "iScansNegatively", &iScansNegatively), 0);
GRIB_CHECK(grib_get_long(gh, "jScansPositively", &jScansPositively), 0);
GRIB_CHECK(grib_get_long(gh, "jPointsAreConsecutive", &jPointsAreConsecutive), 0);
int scanningMode = 0;
scanningMode = 128*iScansNegatively + 64*jScansPositively + 32*jPointsAreConsecutive;
if (cdiDebugExt>=30)
printf("gribapiGetScanningMode(): Scanning mode = %02d (%1d%1d%1d)*32; \n",\
scanningMode,(int)jPointsAreConsecutive,(int)jScansPositively,(int)iScansNegatively);
return scanningMode;
}
void gribapiSetScanningMode(grib_handle *gh, int scanningMode)
{
long iScansNegatively;
long jScansPositively;
long jPointsAreConsecutive;
// 127: reserved for testing; generated test data will be in 64 scanning mode
//if (scanningMode== 127) scanningMode = 64;
iScansNegatively = (scanningMode & 128)/128;
jScansPositively = (scanningMode & 64)/64;
jPointsAreConsecutive = (scanningMode & 32)/32;
if (cdiDebugExt>=30)
{
long paramId, levelTypeId, levelId, uvRelativeToGrid;
GRIB_CHECK(grib_get_long(gh, "uvRelativeToGrid", &uvRelativeToGrid), 0);
GRIB_CHECK(grib_get_long(gh, "indicatorOfParameter", &paramId), 0);
GRIB_CHECK(grib_get_long(gh, "indicatorOfTypeOfLevel", &levelTypeId), 0);
GRIB_CHECK(grib_get_long(gh, "level", &levelId), 0);
printf("gribapiSetScanningMode(): (param,ltype,level) = (%3d,%3d,%4d); Scanning mode = %02d (%1d%1d%1d)*32; uvRelativeToGrid = %02d\n",\
(int)paramId, (int)levelTypeId, (int)levelId,
scanningMode,(int)jPointsAreConsecutive,(int)jScansPositively,(int)iScansNegatively,
(int)uvRelativeToGrid);
}
GRIB_CHECK(my_grib_set_long(gh, "iScansNegatively", iScansNegatively), 0);
GRIB_CHECK(my_grib_set_long(gh, "jScansPositively", jScansPositively), 0);
GRIB_CHECK(my_grib_set_long(gh, "jPointsAreConsecutive", jPointsAreConsecutive), 0);
}
void gribapiSetUvRelativeToGrid(grib_handle *gh, int mode)
{
long uvRelativeToGridMode = mode;
long uvRelativeToGridModeOld;
GRIB_CHECK(grib_get_long(gh, "uvRelativeToGrid", &uvRelativeToGridModeOld), 0);
if (cdiDebugExt>=30)
printf("gribapiSetUvRelativeToGrid(): uvRelativeToGrid: %02d (old) => %02d (new); \n",(int)uvRelativeToGridModeOld,(int)uvRelativeToGridMode);
GRIB_CHECK(my_grib_set_long(gh, "uvRelativeToGrid", uvRelativeToGridMode), 0);
}
/*
TABLE 8. SCANNING MODE FLAG
(GDS Octet 28)
BIT VALUE MEANING
1 0 Points scan in +i direction
1 Points scan in -i direction
2 0 Points scan in -j direction
1 Points scan in +j direction
3 0 Adjacent points in i direction are consecutive
(FORTRAN: (I,J))
1 Adjacent points in j direction are consecutive
(FORTRAN: (J,I))
=> Scanning Mode 0 0 0 0 0 0 0 0 (00 dec) +i, -j; i direction consecutive (row-major order West->East & North->South)
=> Scanning Mode 0 1 0 0 0 0 0 0 (64 dec) +i, +j; i direction consecutive (row-major order West->East & South->North )
=> Scanning Mode 1 1 0 0 0 0 0 0 (96 dec) +i, +j; j direction consecutive (column-major order South->North & West->East )
NOTE: South->North - As if you would plot the data as image on the screen
where [0,0] of the data is the top-left pixel.
grib2ppm LAMH_D11_201302150000_00000_oro | display ppm:-
ImageMagick (display): [0,0] of an image belongs to the top-left pixel
[DEFAULT] : 64 dec
iScansNegatively = 0;
jScansPositively = 1;
jPointsAreConsecutive = 0; => Scanning Mode 64
cdo selindexbox,1,726,100,550 LAMH_D11_201302150000_00000_oro LAMH_D11_201302150000_00000_oro_cropped
grib2ppm LAMH_D11_201302150000_00000_oro_cropped | /usr/bin/display ppm:- &
# ^^^ this image will be missing the souther parts of data
grib2ppm LAMH_D11_201302150000_00000_oro | /usr/bin/display ppm:- &
# ^ full domain data
*/
void verticallyFlipGridDefinitionWhenScanningModeChanged(grib_handle *gh, double yfirst, double ylast, double yinc )
{
/*
Nj = 550;
latitudeOfFirstGridPointInDegrees = -30.8;
latitudeOfLastGridPointInDegrees = 24.1;
iScansNegatively = 0;
jScansPositively = 0;
jPointsAreConsecutive = 0;
jDirectionIncrementInDegrees = 0.1;
When switching from scanning mode 0 <=> 64
yfirst = -30.8 + (550-1)*0.1
yfirst = yfirst + (ysize-1) * yinc
yinc = -1.0*yinc
*/
//long jDim=0;
//GRIB_CHECK(grib_get_long(gh, "Nj", &jDim), 0);
double latitudeOfFirstGridPointInDegrees;
double latitudeOfLastGridPointInDegrees;
double jDirectionIncrementInDegrees;
//GRIB_CHECK(grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &latitudeOfFirstGridPointInDegrees), 0); // yfirst
//GRIB_CHECK(grib_get_double(gh, "latitudeOfLastGridPointInDegrees", &latitudeOfLastGridPointInDegrees), 0); // ylast
//GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &jDirectionIncrementInDegrees), 0); // yinc
if (cdiDebugExt>=10)
{
Message(" BEFORE: yfirst = %f; ylast = %f; yinc = %f; ", yfirst,ylast, yinc);
}
GRIB_CHECK(my_grib_set_double(gh, "latitudeOfFirstGridPointInDegrees", ylast), 0);
GRIB_CHECK(my_grib_set_double(gh, "latitudeOfLastGridPointInDegrees", yfirst), 0);
//yinc *= -1.0; // don't set yinc here ...
//GRIB_CHECK(my_grib_set_double(gh, "jDirectionIncrementInDegrees", yinc), 0);
if (cdiDebugExt>=10)
{
GRIB_CHECK(grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &latitudeOfFirstGridPointInDegrees), 0); // yfirst
GRIB_CHECK(grib_get_double(gh, "latitudeOfLastGridPointInDegrees", &latitudeOfLastGridPointInDegrees), 0); // ylast
GRIB_CHECK(grib_get_double(gh, "jDirectionIncrementInDegrees", &jDirectionIncrementInDegrees), 0); // yinc
Message("CHANGED INTO: yfirst = %f, ylast = %f, yinc = %f",latitudeOfFirstGridPointInDegrees,latitudeOfLastGridPointInDegrees, jDirectionIncrementInDegrees);
}
}
#ifdef HIRLAM_EXTENSIONS
void convertDataScanningMode(int scanModeOUT, double *data, int gridsize, int iDim, int jDim)
{
int i,j;
int idxIN, idxOUT;
// 127: reserved for testing; it will generate test data in 64 scanning mode
if (scanModeOUT== 127) // fill with testdata ...
{
scanModeOUT = 64;
if (cdiDebugExt>=30) printf("convertDataScanningMode(): Generating test data in 64 scanning mode..\n");
for (j=0; j<jDim; j++)
{
int jXiDim = j*iDim;
for (i=0; i<iDim; i++)
{
idxIN = i + jXiDim;
data[idxIN] = (double) (100.0*j +i);
}
}
}
if ( (iDim*jDim)!= gridsize)
{
if (cdiDebugExt>=30) printf("convertDataScanningMode(): ERROR: (iDim*jDim)!= gridsize; (%d * %d) != %d\n", iDim,jDim, gridsize);
return;
}
if (cdiDebugExt>=30) printf("convertDataScanningMode(): scanModeIN=%02d => scanModeOUT=%02d ; where: (iDim * jDim == gridsize) (%d*%d == %d)\n",scanModeIN, scanModeOUT, iDim,jDim, gridsize);
if (cdiDebugExt>=100)
{
printf("convertDataScanningMode(): data IN:\n");
for (j=0; j<jDim; j++)
{
int jXiDim = j*iDim;
for (i=0; i<iDim;