From d2126f055b800cab3dc2a5f529a9e2038eb996ab Mon Sep 17 00:00:00 2001
From: Michal Koutek <koutek@knmi.nl>
Date: Wed, 18 Jan 2017 11:11:23 +0100
Subject: [PATCH] Adding Hirlam_extension branch for new CDO/LIBCDI functions
for HIRLAM NWP model processing
---
.gitignore | 40 +++
configure | 24 ++
src/basetime.c | 11 +-
src/binary.h | 2 +
src/cdi.h | 64 ++++
src/cdi_int.h | 8 +
src/config.h.in | 3 +
src/grb_write.c | 93 ++++++
src/gribapi_utilities.c | 86 ++++++
src/gribapi_utilities.h | 5 +
src/grid.c | 656 ++++++++++++++++++++++++++++++++++++++++
src/grid.h | 17 ++
src/stream_grb.c | 36 ++-
src/stream_gribapi.c | 502 +++++++++++++++++++++++++++++-
src/varscan.c | 6 +
src/vlist.c | 1 +
16 files changed, 1546 insertions(+), 8 deletions(-)
create mode 100644 .gitignore
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 000000000..ad0116e55
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,40 @@
+*.test
+*~
+*.la
+*.lo
+*.o
+*.so
+*.log
+*.trs
+*.dSYM
+*.pyd
+*.pyc
+*.swp
+*.plo
+*.prof
+**/.deps
+**/.libs
+.dirstamp
+stamp-h1
+config.log
+config.status
+config.h
+build
+Makefile
+**/Makefile
+**/config.h
+**/config.lt
+**/config.log
+**/config.status
+*.dirstamp
+cdo.settings
+config.status
+autom4te.cache
+libtool
+src/cdo
+src/stamp-h1
+contrib/bindings/python
+cdi.settings
+.hgignore
+.hg/
+
diff --git a/configure b/configure
index 747cdfcff..dfa81850c 100755
--- a/configure
+++ b/configure
@@ -1567,6 +1567,7 @@ Optional Features:
--enable-ieg Use the ieg library [default=yes]
--enable-all-static build a completely statically linked CDO binary
[default=no]
+ --enable-HIRLAM-extensions This will ensure that HIRLAM_EXTENSIONS is defined.
--enable-mpi Compile with MPI compiler [default=no]
--enable-iso-c-interface
Create Fortran Interface via iso_c_bindings facility
@@ -28263,6 +28264,29 @@ fi
$as_echo "$enable_cgribex" >&6; }
ENABLE_CGRIBEX=$enable_cgribex
+# ----------------------------------------------------------------------
+# Compile interface with HIRLAM extensions
+{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for HIRLAM extensions" >&5
+$as_echo_n "checking for HIRLAM extensions... " >&6; }
+# Check whether --enable-HIRLAM-extensions was given.
+if test "${enable-HIRLAM-extensions+set}" = set; then :
+ enableval=$enable_HIRLAMext; if test "x$enable_HIRLAMext" != 'xno'; then :
+
+$as_echo "#define HIRLAM_EXTENSIONS 1" >>confdefs.h
+
+ enable_HIRLAMext=yes
+fi
+else
+
+$as_echo "#define HIRLAM_EXTENSIONS 1" >>confdefs.h
+
+ enable_HIRLAMext=yes
+fi
+
+{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $enable_HIRLAMext" >&5
+$as_echo "$enable_HIRLAMext" >&6; }
+ENABLE_HIRLAMEXT=$enable_HIRLAMext
+
# ----------------------------------------------------------------------
# Compile interface with internal SERVICE library
{ $as_echo "$as_me:${as_lineno-$LINENO}: checking for SERVICE support" >&5
diff --git a/src/basetime.c b/src/basetime.c
index 3d4d5e366..e4d4b3888 100644
--- a/src/basetime.c
+++ b/src/basetime.c
@@ -3,6 +3,7 @@
#endif
#include <stdio.h>
+#include <stdbool.h>
#include "error.h"
#include "cdi.h"
@@ -16,11 +17,11 @@ void basetimeInit(basetime_t *basetime)
if ( basetime == NULL )
Error("Internal problem! Basetime not allocated.");
- basetime->ncvarid = UNDEFID;
- basetime->ncdimid = UNDEFID;
- basetime->ncvarboundsid = UNDEFID;
- basetime->leadtimeid = UNDEFID;
- basetime->lwrf = 0;
+ basetime->ncvarid = CDI_UNDEFID;
+ basetime->ncdimid = CDI_UNDEFID;
+ basetime->ncvarboundsid = CDI_UNDEFID;
+ basetime->leadtimeid = CDI_UNDEFID;
+ basetime->lwrf = false;
basetime->timevar_cache = NULL;
}
/*
diff --git a/src/binary.h b/src/binary.h
index 26a0a6da3..1bd560f9e 100644
--- a/src/binary.h
+++ b/src/binary.h
@@ -8,7 +8,9 @@
#include <inttypes.h>
#include <stdio.h>
+#ifndef _DTYPES_H
#include "dtypes.h"
+#endif
#ifndef HOST_ENDIANNESS
#ifdef __cplusplus
diff --git a/src/cdi.h b/src/cdi.h
index 69b253d8a..6237c524f 100644
--- a/src/cdi.h
+++ b/src/cdi.h
@@ -9,6 +9,7 @@
#include <stdio.h>
#include <sys/types.h>
+#include <stdbool.h>
#ifdef __cplusplus
extern "C" {
@@ -1230,6 +1231,69 @@ int vlistInqVarSubtype(int vlistID, int varID);
void gribapiLibraryVersion(int *major_version, int *minor_version, int *revision_version);
+#ifdef HIRLAM_EXTENSIONS
+// Make sure that CDO is configured with: --enable-HIRLAM-extensions
+// ./configure --enable-HIRLAM-extensions --prefix=${InstPathCDO} --with-hdf5=${InstPathHDF} --with-netcdf=${InstPathNetCDF} --enable-cdi-lib --enable-python --enable-swig --disable-cgribex --with-grib_api=${InstPathGribAPI} --with-proj=${InstPathProj4}
+
+// WITHOUT extensions:
+// ./configure --prefix=${InstPathCDO} --with-hdf5=${InstPathHDF} --with-netcdf=${InstPathNetCDF} --enable-cdi-lib --enable-python --enable-swig --disable-cgribex --with-grib_api=${InstPathGribAPI} --with-proj=${InstPathProj4}
+
+
+// The following global are defined in libcdi/src/stream_grb.c
+
+ extern int cdoDebugExt; // debug level for the KNMI extensions
+// *** RELATED to GRIB only ***
+extern int cdoGribDataScanningMode; // -1: not used; allowed modes: <0, 64, 96>; Default is 64
+extern int cdoGribChangeModeUvRelativeToGrid; // -1: don't set; 0: set to '0'; 1: set to '1'
+extern int cdoGribUseTimeRangeIndicator; // 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
+// Regarding operation to change parameter identification:
+// -1: don't change; 1: change (cdoGribChangeParID_code, cdoGribChangeParID_ltype, cdoGribChangeParID_lev)
+extern int cdoGribChangeParameterIdentification;
+extern int cdoGribChangeParID_code;
+extern int cdoGribChangeParID_ltype;
+extern int cdoGribChangeParID_lev;
+
+// *** RELATED to scanning mode and uv-grid-relative; defined in grid.c ***
+int gridInqScanningMode(int gridID);
+void gridInqScanningModeBits(int gridID, int *iScansNegatively, int *jScansPositively, int *jPointsAreConsecutive);
+int gridInqUvRelativeToGrid(int gridID);
+void gridSetScanningMode(int gridID, int mode);
+void gridSetScanningModeBits(int gridID, int iScansNegatively, int jScansPositively, int jPointsAreConsecutive);
+void gridSetUvRelativeToGrid(int gridID,int flag);
+
+void convertDataScanningMode(int scanModeOUT, double *data, int gridsize, int iDim, int jDim);
+void streamGrbChangeParameterIdentification(int code, int ltype, int lev);
+
+// Forward declarations:
+/*
+int gribapiGetScanningMode(grib_handle *gh);
+void gribapiSetScanningMode(grib_handle *gh, int scanningMode);
+int gribapiGetUvRelativeToGrid(grib_handle *gh);
+void gribapiSetUvRelativeToGrid(grib_handle *gh, int mode);
+void gribapiChangeParameterIdentification(grib_handle *gh, int code, int ltype, int lev);
+*/
+
+// *** RELATED to de-staggering, grid-sampling and sub-gridding; defined in grid.c ***
+int define_destagered_grid(int gridID_u_stag, int gridID_v_stag, double *destagGridOffsets);
+int define_sample_grid(int gridSrcID, int sampleFactor);
+int define_subgrid_grid(int gridSrcID, int gridIDcurvl, int subI0, int subI1, int subJ0, int subJ1);
+
+const double *gridInqXvalsPtr(int gridID);
+const double *gridInqYvalsPtr(int gridID);
+
+
+// The following global are defined in libcdi/src/stream_cdf.c
+// THIS NetCDF FEATURE is currently NOT included in this CDO branch...
+
+//extern bool CDI_netcdf_MAKE_LONG_AXISNAMES;
+ // Some applications may prefer to have long names for the dimension variables
+ // instead of lon, lat => longitude, latitude will be used..
+//extern bool CDI_netcdf_INCLUDE_NON_STANDARD_DIM_ATTRS;
+ // Some applications (Visualiasation Toolkit VTK, ParaView, IDV) cannot handle these non-standard attributes in NetCDF files.
+ // Set this to false NOT to includes these attributes.
+#endif // HIRLAM_EXTENSIONS
#if defined (__cplusplus)
}
diff --git a/src/cdi_int.h b/src/cdi_int.h
index f558ac682..96190e756 100644
--- a/src/cdi_int.h
+++ b/src/cdi_int.h
@@ -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;
diff --git a/src/config.h.in b/src/config.h.in
index 8f6932de4..79642a169 100644
--- a/src/config.h.in
+++ b/src/config.h.in
@@ -100,6 +100,9 @@
/* GRIB_API library is present if defined to 1 */
#undef HAVE_LIBGRIB_API
+/* Define to 1 to use KNMI extensions */
+#undef HIRLAM_EXTENSIONS
+
/* Define to 1 for IEG interface */
#undef HAVE_LIBIEG
diff --git a/src/grb_write.c b/src/grb_write.c
index dd8cc9f63..d63ca2892 100644
--- a/src/grb_write.c
+++ b/src/grb_write.c
@@ -15,6 +15,20 @@
#include "gribapi.h"
#include "namespace.h"
+#ifdef HIRLAM_EXTENSIONS
+// Forward declarations:
+int gribapiGetScanningMode(grib_handle *gh);
+void gribapiSetScanningMode(grib_handle *gh, int scanningMode);
+int gribapiGetUvRelativeToGrid(grib_handle *gh);
+void gribapiSetUvRelativeToGrid(grib_handle *gh, int mode);
+void gribapiChangeParameterIdentification(grib_handle *gh, int code, int ltype, int lev);
+
+// This function is defined in grb_read.c
+// But we need it here as well.
+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);
+
+#endif // HIRLAM_EXTENSIONS
static
size_t grbEncode(int filetype, int memtype, int varID, int levelID, int vlistID, int gridID, int zaxisID,
@@ -116,6 +130,85 @@ void grbCopyRecord(stream_t * streamptr2, stream_t * streamptr1)
size_t nbytes = recsize;
+#ifdef HIRLAM_EXTENSIONS
+ // Even if you are stream-copy records you might need to change a bit of grib-header !
+ grib_handle *gh = NULL;
+ gh = grib_handle_new_from_message(NULL, (void *) gribbuffer, recsize);
+
+ if (cdoGribChangeModeUvRelativeToGrid>=0)
+ {
+ if (cdoDebugExt>=20) Message("Change ModeUvRelativeToGrid =>%d", cdoGribChangeModeUvRelativeToGrid);
+ gribapiSetUvRelativeToGrid(gh,cdoGribChangeModeUvRelativeToGrid);
+ }
+ if (cdoGribChangeParameterIdentification==1)
+ {
+ if (cdoDebugExt>=20) Message("gribapiChangeParameterIdentification => (%d, %d, %d)", cdoGribChangeParID_code, cdoGribChangeParID_ltype, cdoGribChangeParID_lev);
+ gribapiChangeParameterIdentification(gh, cdoGribChangeParID_code, cdoGribChangeParID_ltype, cdoGribChangeParID_lev);
+ cdoGribChangeParameterIdentification = -1; // after grib attributes have been changed turn it off again
+ }
+
+ if (cdoGribDataScanningMode != -1) // -1: not used; allowed modes: <0, 64, 96>; Default is 64
+ // This will overrule the old scanning mode of the given grid
+ {
+ int scanModeIN;
+ scanModeIN = gribapiGetScanningMode(gh);
+
+ if (scanModeIN != cdoGribDataScanningMode)
+ {
+ if (cdoDebugExt>=20) Message("Change GribDataScanningMode => %d (scanModeIN = %d)", cdoGribDataScanningMode, scanModeIN);
+
+ int status = 0;
+ int gridID;
+ int varID, levelID;
+ int vlistID;
+ int zip;
+ int gridsize, nmiss;
+
+ 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 (cdoDebugExt>=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);
+ return status;
+ }
+ else
+ if (cdoDebugExt>=20) Message("Change GribDataScanningMode => %d (scanModeIN = %d)", cdoGribDataScanningMode, scanModeIN);
+
+ }
+#endif // HIRLAM_EXTENSIONS
+
if ( filetype == FILETYPE_GRB )
{
long unzipsize;
diff --git a/src/gribapi_utilities.c b/src/gribapi_utilities.c
index eabd1d1f6..85d27c1f4 100644
--- a/src/gribapi_utilities.c
+++ b/src/gribapi_utilities.c
@@ -378,11 +378,67 @@ 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);
+ }
+ cdoGribUseTimeRangeIndicator = 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);
}
+#ifdef HIRLAM_EXTENSIONS
+void gribapiGetDataTimeRangeIndicator(grib_handle *gh, int *timeRangeIndicator)
+{
+ long timeRangeInd = 0;
+ // typically 0: for instanteneous fields; 4: for accumulated fields
+ GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &timeRangeInd), 0);
+ *timeRangeIndicator = (int) timeRangeInd;
+}
+
+void gribapiSetDataTimeRangeIndicator(grib_handle *gh, int timeRangeIndicator)
+{
+ GRIB_CHECK(grib_set_long(gh, "timeRangeIndicator", timeRangeIndicator), 0);
+}
+#endif // #ifdef HIRLAM_EXTENSIONS
+
int gribGetDatatype(grib_handle* gribHandle)
{
int datatype;
@@ -752,6 +808,36 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
}
}
+#ifdef HIRLAM_EXTENSIONS
+ if ( gridtype == GRID_LONLAT || gridtype == GRID_CURVILINEAR || 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) */
+
+ GRIB_CHECK(grib_get_long(gh, "uvRelativeToGrid", &grid->uvRelativeToGrid), 0);
+ if (cdoDebugExt>=30)
+ {
+ // indicatorOfParameter=33,indicatorOfTypeOfLevel=105,level
+ long paramId, levelTypeId, levelId;
+ GRIB_CHECK(grib_get_long(gh, "indicatorOfParameter", ¶mId), 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->isRotated = FALSE;
if ( gribapiGetIsRotated(gh) )
{
diff --git a/src/gribapi_utilities.h b/src/gribapi_utilities.h
index b9f20840f..2495f8315 100644
--- a/src/gribapi_utilities.h
+++ b/src/gribapi_utilities.h
@@ -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
diff --git a/src/grid.c b/src/grid.c
index f963581c1..f42ce16f5 100644
--- a/src/grid.c
+++ b/src/grid.c
@@ -167,6 +167,19 @@ void grid_init(grid_t *gridptr)
gridptr->name = NULL;
gridptr->vtable = &cdiGridVtable;
gridptr->extraData = NULL;
+#ifdef HIRLAM_EXTENSIONS
+ 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) */
+ gridptr->uvRelativeToGrid = 0; /* Some models deliver wind U,V relative to the grid-cell */
+#endif // HIRLAM_EXTENSIONS
}
@@ -1793,6 +1806,40 @@ double gridInqYinc(int gridID)
}
}
+ #ifdef HIRLAM_EXTENSIONS__NOT_WORKING_OK
+ /* This is needed to correctly aquire scanning mode when working with NetCDF files.
+ Grib files contain the scanning mode implicitely.
+ Without the following would the conversion NetCDF => grib go wrong for the scanning mode.
+ */
+ double xinc = gridInqXinc(gridID); gridptr->xinc = xinc;
+
+ if (xinc > 0)
+ gridptr->iScansNegatively = 0;
+ else
+ gridptr->iScansNegatively = 1; // this is rare but possible ...
+
+ /* .. TO BE DONE: Detection of array ordening with NetCDF ...
+ gridptr->jPointsAreConsecutive = 0; ~ C-style: Default
+
+ C-style ordering (column-major) ~ the raw elements are stored first.
+ FORTRAN-style ordering (row-major) ~ the column elements are stored first
+ */
+ if (yinc > 0)
+ gridptr->jScansPositively = 1;
+ else
+ gridptr->jScansPositively = 0;
+ gridptr->scanningMode = 128*gridptr->iScansNegatively + 64*gridptr->jScansPositively + 32*gridptr->jPointsAreConsecutive;
+ if ( cdoDebugExt>=10 )
+ printf("gridInqYinc(): xinc=%f; yinc=%f; scanningMode=%d\n", xinc, yinc, gridptr->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) */
+ #endif // HIRLAM_EXTENSIONS
+
return yinc;
}
@@ -2404,6 +2451,17 @@ int gridCompare(int gridID, const grid_t *grid)
}
}
+#ifdef HIRLAM_EXTENSIONS
+ if ( (grid->scanningMode != gridInqScanningMode(gridID)) || (grid->uvRelativeToGrid != gridInqUvRelativeToGrid(gridID)) )
+ {
+ // often grid definition may differ in UV-relativeToGrid, simply new or other GRID will be used in CDO
+ if ( cdoDebugExt>=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) );
+ differ = 1;
+ }
+#endif // HIRLAM_EXTENSIONS
+
return differ;
}
@@ -2476,6 +2534,10 @@ int gridCompareP ( void * gridptr1, void * gridptr2 )
if ( IS_NOT_EQUAL(g1->xpole , g2->xpole) ) return differ;
if ( IS_NOT_EQUAL(g1->ypole , g2->ypole) ) return differ;
if ( IS_NOT_EQUAL(g1->angle , g2->angle) ) return differ;
+#ifdef HIRLAM_EXTENSIONS
+ if ( IS_NOT_EQUAL(g1->scanningMode , g2->scanningMode) ) return differ;
+ if ( IS_NOT_EQUAL(g1->uvRelativeToGrid , g2->uvRelativeToGrid) ) return differ;
+#endif // HIRLAM_EXTENSIONS
const double *restrict g1_xvals = g1->vtable->inqXValsPtr(g1),
*restrict g2_xvals = g2->vtable->inqXValsPtr(g2);
@@ -2850,6 +2912,13 @@ int gridGenerate(const grid_t *grid)
gridptr->lcc_scanflag = grid->lcc_scanflag;
gridptr->number = grid->number;
gridptr->position = grid->position;
+ #ifdef HIRLAM_EXTENSIONS
+ gridptr->scanningMode = grid->scanningMode;
+ gridptr->iScansNegatively = grid->iScansNegatively;
+ gridptr->jScansPositively = grid->jScansPositively;
+ gridptr->jPointsAreConsecutive = grid->jPointsAreConsecutive;
+ gridptr->uvRelativeToGrid = grid->uvRelativeToGrid;
+#endif // HIRLAM_EXTENSIONS
memcpy(gridptr->uuid, grid->uuid, CDI_UUID_SIZE);
if ( gridtype == GRID_UNSTRUCTURED && grid->reference )
gridDefReference(gridID, grid->reference);
@@ -3670,6 +3739,12 @@ static void gridPrintKernel(grid_t * gridptr, int index, int opt, FILE *fp)
if ( uuidOfHGridStr[0] != 0 && strlen(uuidOfHGridStr) == 36 )
fprintf(fp, "uuid = %s\n", uuidOfHGridStr);
}
+#ifdef HIRLAM_EXTENSIONS
+ {
+ int scanningMode = gridInqScanningMode(gridID);
+ fprintf(fp, "scanningMode = %d\n", scanningMode);
+ }
+#endif // HIRLAM_EXTENSIONS
if ( gridptr->mask )
{
@@ -4184,8 +4259,13 @@ gridTxCode ()
return GRID;
}
+#ifdef HIRLAM_EXTENSIONS
+enum { gridNint = 33,
+ gridNdouble = 29,
+#else
enum { gridNint = 28,
gridNdouble = 24,
+#endif // HIRLAM_EXTENSIONS
gridHasMaskFlag = 1 << 0,
gridHasGMEMaskFlag = 1 << 1,
gridHasXValsFlag = 1 << 2,
@@ -4379,6 +4459,14 @@ gridUnpack(char * unpackBuffer, int unpackBufferSize,
memberMask = intBuffer[25];
gridP->xstdname = xystdname_tab[intBuffer[26]][0];
gridP->ystdname = xystdname_tab[intBuffer[27]][1];
+#ifdef HIRLAM_EXTENSIONS
+ gridP->iScansNegatively = intBuffer[28];
+ gridP->jScansPositively = intBuffer[29];
+ gridP->jPointsAreConsecutive = intBuffer[30];
+ gridP->scanningMode = intBuffer[31];
+ gridP->uvRelativeToGrid = intBuffer[32];
+#endif // HIRLAM_EXTENSIONS
+
}
if (memberMask & gridHasRowLonFlag)
@@ -4424,6 +4512,14 @@ gridUnpack(char * unpackBuffer, int unpackBufferSize,
gridP->xpole = doubleBuffer[21];
gridP->ypole = doubleBuffer[22];
gridP->angle = doubleBuffer[23];
+#ifdef HIRLAM_EXTENSIONS
+ gridP->iScansNegatively = doubleBuffer[24];
+ gridP->jScansPositively = doubleBuffer[25];
+ gridP->jPointsAreConsecutive = doubleBuffer[26];
+ gridP->scanningMode = doubleBuffer[27];
+ gridP->uvRelativeToGrid = doubleBuffer[28];
+#endif // HIRLAM_EXTENSIONS
+
}
int irregular = gridP->type == GRID_UNSTRUCTURED
@@ -4584,6 +4680,13 @@ gridPack(void * voidP, void * packBuffer, int packBufferSize,
- xystdname_tab);
intBuffer[27] = (int)((const char (*)[2][24])gridP->ystdname
- (const char (*)[2][24])xystdname_tab[0][1]);
+#ifdef HIRLAM_EXTENSIONS
+ intBuffer[28] = gridP->iScansNegatively;
+ intBuffer[29] = gridP->jScansPositively;
+ intBuffer[30] = gridP->jPointsAreConsecutive;
+ intBuffer[31] = gridP->scanningMode;
+ intBuffer[32] = gridP->uvRelativeToGrid;
+#endif // HIRLAM_EXTENSIONS
serializePack(intBuffer, gridNint, DATATYPE_INT,
packBuffer, packBufferSize, packBufferPos, context);
@@ -4630,6 +4733,14 @@ gridPack(void * voidP, void * packBuffer, int packBufferSize,
doubleBuffer[21] = gridP->xpole;
doubleBuffer[22] = gridP->ypole;
doubleBuffer[23] = gridP->angle;
+#ifdef HIRLAM_EXTENSIONS
+ doubleBuffer[24] = gridP->iScansNegatively;
+ doubleBuffer[25] = gridP->jScansPositively;
+ doubleBuffer[26] = gridP->jPointsAreConsecutive;
+ doubleBuffer[27] = gridP->scanningMode;
+ doubleBuffer[28] = gridP->uvRelativeToGrid;
+#endif // HIRLAM_EXTENSIONS
+
serializePack(doubleBuffer, gridNdouble, DATATYPE_FLT64,
packBuffer, packBufferSize, packBufferPos, context);
@@ -4876,6 +4987,551 @@ const struct gridVirtTable cdiGridVtable
.inqYBoundsPtr = gridInqYBoundsPtrSerial,
};
+#ifdef HIRLAM_EXTENSIONS
+
+extern int cdoDebugExt; // defined in cdo.c
+
+grid_t *grid_to_pointer(int gridID)
+{
+ grid_t *gridptr = NULL;
+
+ gridptr = ( grid_t *) reshGetVal ( gridID, &gridOps );
+ return gridptr;
+}
+
+static
+void gridCheckPtr(const char *caller, int gridID, grid_t *gridptr)
+{
+ if ( gridptr == NULL )
+ Errorc("grid %d undefined!", gridID);
+}
+
+#define grid_check_ptr(gridID, gridptr) gridCheckPtr(__func__, gridID, gridptr)
+
+
+
+/*
+@Function define_destagered_grid
+@Title Define a de-staggered grid for U and V
+
+@Prototype int define_destagered_grid(int gridID_u_stag, int gridID_v_stag, double *destagGridOffsets)
+@Parameter
+ @Item grid_u_stag Staggered grid of u-wind component
+ @Item grid_v_stag Staggered grid of v-wind component
+ @Item grid_uv_destag Destaggered grid of uv-wind
+
+@Description
+The function @func{define_destagered_grid} defines a de-staggered grid for U and V
+
+@EndFunction
+*/
+int define_destagered_grid(int gridID_u_stag, int gridID_v_stag, double *destagGridOffsets)
+{
+/* Example of horizontal grids (Hirlam LAMH_D11):
+ U : lonlat > size : dim = 399300 nlon = 726 nlat = 550
+ rlon : first = -30.15 last = 42.35 inc = 0.1 degrees
+ rlat : first = -30.8 last = 24.1 inc = 0.1 degrees
+ northpole : lon = -195 lat = 30
+ V : lonlat > size : dim = 399300 nlon = 726 nlat = 550
+ rlon : first = -30.2 last = 42.3 inc = 0.1 degrees
+ rlat : first = -30.75 last = 24.15 inc = 0.1 degrees
+ northpole : lon = -195 lat = 30
+=> RESULT:
+ R : lonlat > size : dim = 399300 nlon = 726 nlat = 550
+ rlon : first = -30.2 last = 42.3 inc = 0.1 degrees
+ rlat : first = -30.8 last = 24.1 inc = 0.1 degrees
+ northpole : lon = -195 lat = 30
+*/
+ if ( cdoDebugExt ) printf("cdo uvDestag: define_destagered_grid(gridID_u=%d,gridID_v=%d,destagGridOffsets(%02.1f,%02.1f)) ...\n",gridID_u_stag, gridID_v_stag, destagGridOffsets[0],destagGridOffsets[1]);
+
+ if ( cdoDebugExt>1 ) {
+ gridPrint(gridID_u_stag,1,0);
+ gridPrint(gridID_v_stag,1,0);
+ }
+
+ grid_t *grid_u_stag, *grid_v_stag;
+ grid_t *grid_uv_destag;
+ int gridID_uv_destag;
+
+ grid_u_stag = grid_to_pointer(gridID_u_stag);
+ grid_check_ptr(gridID_u_stag, grid_u_stag);
+
+ grid_v_stag = grid_to_pointer(gridID_v_stag);
+ grid_check_ptr(gridID_v_stag, grid_v_stag);
+
+
+ const double *xvals_U = gridInqXvalsPtr(gridID_u_stag);
+ const double *yvals_U = gridInqYvalsPtr(gridID_u_stag);
+ const double *xvals_V = gridInqXvalsPtr(gridID_v_stag);
+ const double *yvals_V = gridInqYvalsPtr(gridID_v_stag);
+
+ int gridXsize = gridInqXsize(gridID_u_stag);
+ int gridYsize = gridInqYsize(gridID_u_stag);
+
+ double xfirst_U = gridInqXval(gridID_u_stag,0); // staggered grid of u-wind
+ double yfirst_U = gridInqYval(gridID_u_stag,0);
+ double xlast_U = gridInqXval(gridID_u_stag, gridXsize-1);
+ double ylast_U = gridInqYval(gridID_u_stag, gridYsize-1);
+ double xfirst_V = gridInqXval(gridID_v_stag,0); // staggered grid of v-wind
+ double yfirst_V = gridInqYval(gridID_v_stag,0);
+ double xlast_V = gridInqXval(gridID_v_stag, gridXsize-1);
+ double ylast_V = gridInqYval(gridID_v_stag, gridYsize-1);
+ double xinc = gridInqXinc(gridID_u_stag);
+ double yinc = gridInqYinc(gridID_u_stag);
+
+
+ gridID_uv_destag = gridDuplicate(gridID_u_stag);
+
+ grid_uv_destag = grid_to_pointer(gridID_u_stag);
+ grid_check_ptr(gridID_uv_destag, grid_uv_destag);
+
+ grid_uv_destag->scanningMode = grid_u_stag->scanningMode;
+ grid_uv_destag->iScansNegatively = grid_u_stag->iScansNegatively;
+ grid_uv_destag->jScansPositively = grid_u_stag->jScansPositively;
+ grid_uv_destag->jPointsAreConsecutive = grid_u_stag->jPointsAreConsecutive;
+ grid_uv_destag->uvRelativeToGrid = grid_u_stag->uvRelativeToGrid;
+
+ gridPrint(gridID_uv_destag, 1,0);
+
+ if ( cdoDebugExt )
+ {
+ printf("cdo uvDestag: define_destagered_grid(): (grid_u_stag->xdef=%d, grid_u_stag->ydef=%d); (gridXsize=%d, gridYsize=%d)\n", grid_u_stag->xdef, grid_u_stag->ydef, gridXsize, gridYsize);
+ printf("cdo uvDestag: define_destagered_grid(): (xfirst_U = %3.2f; yfirst_U = %3.2f); (xfirst_V = %3.2f; yfirst_V = %3.2f);\n",xfirst_U,yfirst_U,xfirst_V,yfirst_V);
+ printf("define_cdo uvDestag: destagered_grid(): (xlast_U = %3.2f; ylast_U = %3.2f); (xlast_V = %3.2f; ylast_V = %3.2f);\n",xlast_U,ylast_U,xlast_V,ylast_V);
+ // NOTE: xdef, ydef, xfirst, yfirst, xlast, ylast == 0 ! Not defined after reading the GRIB file with GRIBAPI
+ //printf("define_destagered_grid(): (xfirst_U=%3.2f;yfirst_U=%3.2f); (xfirst_V=%3.2f;yfirst_V=%3.2f);\n",grid_u_stag->xfirst, grid_u_stag->xfirst, grid_v_stag->xfirst, grid_v_stag->yfirst);
+ //printf("define_destagered_grid(): (xlast_U =%3.2f;ylast_U =%3.2f); (xlast_V =%3.2f;ylast_V =%3.2f);\n",grid_u_stag->xlast, grid_u_stag->ylast, grid_v_stag->xlast, grid_v_stag->ylast);
+ }
+ grid_uv_destag->xinc = xinc;
+ grid_uv_destag->yinc = yinc;
+
+ if ( (destagGridOffsets[0]==-0.5) && (destagGridOffsets[1]==-0.5) )
+ {
+ grid_uv_destag->xfirst = xfirst_V;
+ grid_uv_destag->xlast = xlast_V;
+
+ grid_uv_destag->yfirst = yfirst_U;
+ grid_uv_destag->ylast = ylast_U;
+ }
+ else
+ if ( (destagGridOffsets[0]==0.5) && (destagGridOffsets[1]==0.5) )
+ {
+ grid_uv_destag->xfirst = xfirst_V + xinc*destagGridOffsets[0];
+ grid_uv_destag->xlast = xlast_V + xinc*destagGridOffsets[0];
+
+ grid_uv_destag->yfirst = yfirst_U + yinc*destagGridOffsets[1];
+ grid_uv_destag->ylast = ylast_U + yinc*destagGridOffsets[1];
+ }
+ else
+ Error("cdo uvDestag: define_destagered_grid() Unsupported destaggered grid offsets! We support only: (-0.5,-0.5) or (0.5,0.5)");
+
+ double *xvals = (double *) malloc(grid_uv_destag->xsize*sizeof(double));
+ gridGenXvals(grid_uv_destag->xsize, grid_uv_destag->xfirst, grid_uv_destag->xlast, grid_uv_destag->xinc, xvals);
+ gridDefXvals(gridID_uv_destag, xvals); free(xvals);
+ //gridDefXinc(gridID_uv_destag, grid_uv_destag->xinc); was not needed ..
+
+ double *yvals = (double *) malloc(grid_uv_destag->ysize*sizeof(double));
+ gridGenYvals(grid_uv_destag->type, grid_uv_destag->ysize, grid_uv_destag->yfirst, grid_uv_destag->ylast, grid_uv_destag->yinc, yvals);
+ gridDefYvals(gridID_uv_destag, yvals); free(yvals);
+ //gridDefYinc(gridID_uv_destag, grid_uv_destag->yinc); was not needed ..
+
+ if ( cdoDebugExt )
+ {
+ printf("cdo uvDestag: define_destagered_grid(): \n");
+ gridPrint(gridID_uv_destag, 1,0);
+ }
+
+ return gridID_uv_destag;
+}
+
+/*
+@Function define_sample_grid
+@Title Define a sampled grid of another grid
+
+@Prototype int define_sample_grid(int gridSrcID, int sampleFactor)
+@Parameter
+ @Item gridSrcID Source grid
+ @Item sampleFactor sampleFactor; typically 2,3,4 ...
+
+@Description
+The function @func{define_sample_grid} defines a sampled grid of another grid
+
+@EndFunction
+*/
+int define_sample_grid(int gridSrcID, int sampleFactor)
+{
+/* Example of horizontal grids (Harmonie HARM36_L25):
+ #
+ # gridID 2
+ #
+ gridtype = lcc
+ gridsize = 622521
+ xsize = 789
+ ysize = 789
+ originLon = -7.89
+ originLat = 42.935
+ lonParY = 0
+ lat1 = 52.5
+ lat2 = 52.5
+ xinc = 2500
+ yinc = 2500
+ projection = northpole
+=> RESULT:
+ #
+ # gridID 2
+ #
+ gridtype = lcc
+ gridsize = 156025
+ xsize = 395
+ ysize = 395
+ originLon = -7.89
+ originLat = 42.935
+ lonParY = 0
+ lat1 = 52.5
+ lat2 = 52.5
+ xinc = 5000
+ yinc = 5000
+ projection = northpole
+*/
+ if ( cdoDebugExt )
+ printf("cdo SampleGrid: define_sample_grid(gridSrcID=%d, sampleFactor=%d) ...\n",gridSrcID, sampleFactor);
+
+ int gridXsize = gridInqXsize(gridSrcID);
+ int gridYsize = gridInqYsize(gridSrcID);
+
+ if ( (sampleFactor<1) || (gridXsize<1) || (gridYsize<1) || (sampleFactor > (gridXsize/4) ) || (sampleFactor > (gridYsize/4)) )
+ Error("cdo SampleGrid: define_sample_grid() Unsupported sampleFactor (%d)! Note that: gridXsize = %d, gridYsize = %d", sampleFactor, gridXsize, gridYsize);
+
+ if ( cdoDebugExt>20 ) {
+ gridPrint(gridSrcID,1,0);
+ }
+
+ grid_t *grid_src;
+ grid_t *grid_sampled;
+ int gridID_sampled;
+
+ grid_src = grid_to_pointer(gridSrcID);
+ grid_check_ptr(gridSrcID, grid_src);
+
+ //const double *xvals = gridInqXvalsPtr(gridSrcID);
+ //const double *yvals = gridInqYvalsPtr(gridSrcID);
+ /*
+ double xfirst = gridInqXval(gridSrcID,0); // staggered grid of u-wind
+ double yfirst = gridInqYval(gridSrcID,0);
+ double xlast = gridInqXval(gridSrcID, gridXsize-1);
+ double ylast = gridInqYval(gridSrcID, gridYsize-1);
+ double xinc = gridInqXinc(gridSrcID);
+ double yinc = gridInqYinc(gridSrcID);
+ */
+
+ gridID_sampled = gridDuplicate(gridSrcID);
+ grid_sampled = grid_to_pointer(gridID_sampled);
+ grid_check_ptr(gridID_sampled, grid_sampled);
+
+ grid_sampled->scanningMode = grid_src->scanningMode;
+ grid_sampled->iScansNegatively = grid_src->iScansNegatively;
+ grid_sampled->jScansPositively = grid_src->jScansPositively;
+ grid_sampled->jPointsAreConsecutive = grid_src->jPointsAreConsecutive;
+ grid_sampled->uvRelativeToGrid = grid_src->uvRelativeToGrid;
+
+ grid_sampled->xsize = (grid_src->xsize + (sampleFactor-1)) / sampleFactor; // HARM36_L25: (789 + 2-1) / 2 = 395
+ grid_sampled->ysize = (grid_src->ysize + (sampleFactor-1)) / sampleFactor;
+ grid_sampled->size = grid_sampled->xsize * grid_sampled->ysize;
+
+ // for the case of Lambert projection ...
+ grid_sampled->lcc_xinc = grid_src->lcc_xinc * sampleFactor;
+ grid_sampled->lcc_yinc = grid_src->lcc_yinc * sampleFactor;
+
+ grid_sampled->xinc = grid_src->xinc * sampleFactor;
+ grid_sampled->yinc = grid_src->yinc * sampleFactor;
+
+ double *xvals = (double *) malloc(grid_sampled->xsize*sizeof(double));
+ //void gridGenXvals(int xsize, double xfirst, double xlast, double xinc, double *restrict xvals)
+ gridGenXvals(grid_sampled->xsize, grid_sampled->xfirst, grid_sampled->xlast, grid_sampled->xinc, xvals);
+ gridDefXvals(gridID_sampled, xvals); free(xvals);
+ //gridDefXinc(gridID_sampled, grid_sampled->xinc); was not needed ..
+
+ double *yvals = (double *) malloc(grid_sampled->ysize*sizeof(double));
+ //void gridGenYvals(int gridtype, int ysize, double yfirst, double ylast, double yinc, double *restrict yvals)
+ gridGenYvals(grid_sampled->type, grid_sampled->ysize, grid_sampled->yfirst, grid_sampled->ylast, grid_sampled->yinc, yvals);
+ gridDefYvals(gridID_sampled, yvals); free(yvals);
+ //gridDefYinc(gridID_sampled, grid_sampled->yinc); was not needed ..
+
+ if ( grid_sampled->type == GRID_LCC )
+ gridDefLCC( gridID_sampled, grid_sampled->lcc_originLon, grid_sampled->lcc_originLat, grid_sampled->lcc_lonParY,
+ grid_sampled->lcc_lat1, grid_sampled->lcc_lat2, grid_sampled->lcc_xinc, grid_sampled->lcc_yinc,
+ grid_sampled->lcc_projflag, grid_sampled->lcc_scanflag);
+
+ if ( cdoDebugExt>20 )
+ {
+ printf("cdo SampleGrid: define_sample_grid(): \n");
+ gridPrint(gridID_sampled, 1,0);
+ }
+
+ return gridID_sampled;
+}
+
+
+/*
+@Function define_subgrid_grid
+@Title Define a sub-grid of another grid (LCC)
+
+@Prototype int define_subgrid_grid(int gridIDsrc, int subI0, int subI1, int subJ0, int subJ1)
+@Parameter
+ @Item gridSrcID Source grid
+ @Item subI0,subI1, subJ0, subJ1 Sub-grid indices
+
+@Description
+The function @func{define_subgrid_grid} defines a sub-grid of another grid (LCC)
+
+@EndFunction
+*/
+int define_subgrid_grid(int gridSrcID, int gridIDcurvl, int subI0, int subI1, int subJ0, int subJ1)
+{
+/* Example of horizontal grids (Harmonie HARM36_L25):
+ #
+ # gridID 2
+ #
+ gridtype = lcc
+ gridsize = 622521
+ xsize = 789
+ ysize = 789
+ originLon = -7.89
+ originLat = 42.935
+ lonParY = 0
+ lat1 = 52.5
+ lat2 = 52.5
+ xinc = 2500
+ yinc = 2500
+ projection = northpole
+=> RESULT:
+ #
+ # gridID 2
+ #
+ gridtype = lcc
+ gridsize = 156025
+ xsize = 350
+ ysize = 350
+ originLon = ...
+ originLat = ...
+ lonParY = 0
+ lat1 = 52.5
+ lat2 = 52.5
+ xinc = 2500
+ yinc = 2500
+ projection = northpole
+*/
+ if ( cdoDebugExt )
+ printf("cdo SampleGrid: define_subgrid_grid(gridSrcID=%d, (subI0,subI1,subJ0,subJ1) =(%d,%d,%d,%d) ...\n",gridSrcID, subI0,subI1, subJ0, subJ1 );
+
+
+ int gridXsize = gridInqXsize(gridSrcID);
+ int gridYsize = gridInqYsize(gridSrcID);
+ int maxIndexI = gridXsize-1;
+ int maxIndexJ = gridYsize-1;
+
+
+ if ( \
+ (subI0<0) || (subI0>maxIndexI) || \
+ (subI1<=subI0) || (subI1>maxIndexI) ||
+ (subJ0<0) || (subJ0>maxIndexJ) || \
+ (subJ1<=subJ0) || (subJ1>maxIndexJ) )
+ Error("cdo SampleGrid: define_subgrid_grid() Incorrect subgrid specified! (subI0,subI1,subJ0,subJ1) =(%d,%d,%d,%d) Note that: gridXsize = %d, gridYsize = %d", subI0,subI1, subJ0, subJ1, gridXsize, gridYsize);
+
+
+ grid_t *grid_src;
+ grid_t *grid_sampled;
+ int gridID_sampled;
+
+ double originLon, originLat;
+
+ grid_src = grid_to_pointer(gridSrcID);
+ grid_check_ptr(gridSrcID, grid_src);
+
+ if ( grid_src->type != GRID_LCC )
+ Error("cdo SampleGrid: define_subgrid_grid() Error; Only LCC grid is supported; use selindexbox");
+
+
+ double lonParY; double lat1; double lat2; double xinc; double yinc; int projflag; int scanflag;
+ gridInqLCC(gridSrcID, &originLon, &originLat, &lonParY, &lat1, &lat2, &xinc, &yinc, &projflag, &scanflag);
+
+ if ( cdoDebugExt>20 ) {
+ gridPrint(gridSrcID,1,0);
+ }
+
+ if (cdoDebugExt)
+ {
+ Message("cdo SampleGrid: define_subgrid_grid() Original LCC grid:");
+ Message("grid Xsize %d, grid Ysize %d", gridInqXsize(gridSrcID), gridInqYsize(gridSrcID));
+ Message("originLon %4.3f, originLon %4.3f", originLon, originLat);
+ Message("grid Xinc %4.3f, grid Yinc %4.3f", xinc, yinc);
+ }
+ originLon = gridInqXval(gridIDcurvl, 0);
+ originLat = gridInqYval(gridIDcurvl, 0);
+
+ if (cdoDebugExt)
+ {
+ Message("\ncdo SampleGrid: define_subgrid_grid() Original LCC grid as curvilinear (with lats-lons computed):");
+ Message("grid Xsize %d, grid Ysize %d", gridInqXsize(gridIDcurvl), gridInqYsize(gridIDcurvl));
+ Message("grid Xfirst %4.3f, grid Yfirst %4.3f", gridInqXval(gridIDcurvl, 0), gridInqYval(gridIDcurvl, 0));
+ Message("grid Xlast %4.3f, grid Ylast %4.3f", gridInqXval(gridIDcurvl, gridInqSize(gridIDcurvl) -1), gridInqYval(gridIDcurvl, gridInqSize(gridIDcurvl) -1));
+ Message("originLon %4.3f, originLon %4.3f", originLon, originLat);
+ }
+
+
+
+ //const double *xvals = gridInqXvalsPtr(gridSrcID);
+ //const double *yvals = gridInqYvalsPtr(gridSrcID);
+ /*
+ double xfirst = gridInqXval(gridSrcID,0); // staggered grid of u-wind
+ double yfirst = gridInqYval(gridSrcID,0);
+ double xlast = gridInqXval(gridSrcID, gridXsize-1);
+ double ylast = gridInqYval(gridSrcID, gridYsize-1);
+ double xinc = gridInqXinc(gridSrcID);
+ double yinc = gridInqYinc(gridSrcID);
+ */
+
+ gridID_sampled = gridDuplicate(gridSrcID);
+ grid_sampled = grid_to_pointer(gridID_sampled);
+ grid_check_ptr(gridID_sampled, grid_sampled);
+
+ grid_sampled->scanningMode = grid_src->scanningMode;
+ grid_sampled->iScansNegatively = grid_src->iScansNegatively;
+ grid_sampled->jScansPositively = grid_src->jScansPositively;
+ grid_sampled->jPointsAreConsecutive = grid_src->jPointsAreConsecutive;
+ grid_sampled->uvRelativeToGrid = grid_src->uvRelativeToGrid;
+
+ grid_sampled->xsize = subI1 - subI0 + 1;
+ grid_sampled->ysize = subJ1 - subJ0 + 1;
+ grid_sampled->size = grid_sampled->xsize * grid_sampled->ysize;
+ grid_sampled->xinc = grid_src->xinc;
+ grid_sampled->yinc = grid_src->yinc;
+
+ /* LCC: Lambert Conformal Conic
+ double lcc_originLon;
+ double lcc_originLat;
+ double lcc_lonParY;
+ double lcc_lat1;
+ double lcc_lat2;
+ double lcc_xinc;
+ double lcc_yinc;
+ int lcc_projflag;
+ int lcc_scanflag;
+ int lcc_defined; */
+
+ originLon = gridInqXval(gridIDcurvl, subJ0*gridInqXsize(gridIDcurvl) + subI0 );
+ originLat = gridInqYval(gridIDcurvl, subJ0*gridInqXsize(gridIDcurvl) + subI0 );
+
+ if (cdoDebugExt)
+ {
+ Message("\ncdo SampleGrid: define_subgrid_grid() Sub-grid:");
+ Message("grid Xsize %d, grid Ysize %d", gridInqXsize(gridID_sampled), gridInqYsize(gridID_sampled));
+ Message("originLon %4.3f, originLon %4.3f", originLon, originLat);
+ }
+
+ grid_sampled->lcc_originLon = originLon;
+ grid_sampled->lcc_originLat = originLat;
+ grid_sampled->lcc_lat1 = grid_src->lcc_lat1;
+ grid_sampled->lcc_lat2 = grid_src->lcc_lat2;
+ grid_sampled->lcc_xinc = grid_src->lcc_xinc;
+ grid_sampled->lcc_yinc = grid_src->lcc_yinc;
+ grid_sampled->lcc_projflag = grid_src->lcc_projflag;
+ grid_sampled->lcc_scanflag = grid_src->lcc_scanflag;
+ grid_sampled->lcc_defined = grid_src->lcc_defined;
+
+ if ( grid_sampled->type == GRID_LCC )
+ gridDefLCC(gridID_sampled, grid_sampled->lcc_originLon, grid_sampled->lcc_originLat, grid_sampled->lcc_lonParY,
+ grid_sampled->lcc_lat1, grid_sampled->lcc_lat2, grid_sampled->lcc_xinc, grid_sampled->lcc_yinc,
+ grid_sampled->lcc_projflag, grid_sampled->lcc_scanflag);
+
+
+ if ( cdoDebugExt>20 )
+ {
+ printf("cdo SampleGrid: define_subgrid_grid(): \n");
+ gridPrint(gridID_sampled, 1,0);
+ }
+
+ return gridID_sampled;
+}
+
+
+
+int gridInqScanningMode(int gridID)
+{
+ grid_t *gridptr;
+
+ gridptr = grid_to_pointer(gridID);
+
+ grid_check_ptr(gridID, gridptr);
+ 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 gridInqScanningModeBits(int gridID, int *iScansNegatively, int *jScansPositively, int *jPointsAreConsecutive)
+{
+ grid_t *gridptr;
+
+ gridptr = grid_to_pointer(gridID);
+
+ grid_check_ptr(gridID, gridptr);
+
+ *iScansNegatively = gridptr->iScansNegatively;
+ *jScansPositively = gridptr->jScansPositively;
+ *jPointsAreConsecutive = 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) */
+ return;
+}
+
+int gridInqUvRelativeToGrid(int gridID)
+{
+ grid_t *gridptr;
+
+ gridptr = grid_to_pointer(gridID);
+
+ grid_check_ptr(gridID, gridptr);
+
+ return ( gridptr->uvRelativeToGrid );
+}
+
+void gridSetScanningMode(int gridID, int mode)
+{
+ grid_t *gridptr;
+ gridptr = grid_to_pointer(gridID);
+ grid_check_ptr(gridID, gridptr);
+
+ gridptr->scanningMode = mode;
+}
+
+void gridSetScanningModeBits(int gridID, int iScansNegatively, int jScansPositively, int jPointsAreConsecutive)
+{
+ grid_t *gridptr;
+ gridptr = grid_to_pointer(gridID);
+ grid_check_ptr(gridID, gridptr);
+
+ gridptr->iScansNegatively = iScansNegatively;
+ gridptr->jScansPositively = jScansPositively;
+ gridptr->jPointsAreConsecutive = jPointsAreConsecutive;
+}
+
+void gridSetUvRelativeToGrid(int gridID,int flag)
+{
+ grid_t *gridptr;
+ gridptr = grid_to_pointer(gridID);
+ grid_check_ptr(gridID, gridptr);
+
+ gridptr->uvRelativeToGrid = flag;
+}
+#endif // HIRLAM_EXTENSIONS
+
/*
* Local Variables:
* c-file-style: "Java"
diff --git a/src/grid.h b/src/grid.h
index 67e4b73b4..5e8058899 100644
--- a/src/grid.h
+++ b/src/grid.h
@@ -107,6 +107,19 @@ struct grid_t {
char ylongname[CDI_MAX_NAME];
char xunits[CDI_MAX_NAME];
char yunits[CDI_MAX_NAME];
+#ifdef HIRLAM_EXTENSIONS
+ 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) */
+ long uvRelativeToGrid; /* Some models deliver wind U,V relative to the grid-cell */
+#endif // HIRLAM_EXTENSIONS
char *name;
const struct gridVirtTable *vtable;
void *extraData;
@@ -125,6 +138,10 @@ unsigned cdiGridCount(void);
const double *gridInqXvalsPtr(int gridID);
const double *gridInqYvalsPtr(int gridID);
+#ifdef HIRLAM_EXTENSIONS
+grid_t *grid_to_pointer(int gridID);
+#endif // HIRLAM_EXTENSIONS
+
const double *gridInqXboundsPtr(int gridID);
const double *gridInqYboundsPtr(int gridID);
const double *gridInqAreaPtr(int gridID);
diff --git a/src/stream_grb.c b/src/stream_grb.c
index a56313a9c..295c2a4e1 100644
--- a/src/stream_grb.c
+++ b/src/stream_grb.c
@@ -10,7 +10,41 @@
#include "file.h"
#include "cgribex.h" /* gribZip gribGetZip gribGinfo */
#include "gribapi.h"
-
+#ifdef HIRLAM_EXTENSIONS
+#include "grib_api.h"
+
+int cdoDebugExt = 0; // Debug level for the KNMI extensions
+// *** RELATED to GRIB only ***
+int cdoGribDataScanningMode = -1; // -1: not used; allowed modes: <0, 64, 96>; Default is 64
+int cdoGribChangeModeUvRelativeToGrid = -1; // -1: don't set; 0: set to '0'; 1: set to '1'
+int cdoGribUseTimeRangeIndicator = 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
+extern void gribapiSetUvRelativeToGrid(grib_handle *gh, int mode);
+extern int gribapiGetScanningMode(grib_handle *gh);
+extern void gribapiSetScanningMode(grib_handle *gh, int scanningMode);
+
+// Regarding operation to change parameter identification:
+// -1: don't change; 1: change (cdoGribChangeParID_code, cdoGribChangeParID_ltype, cdoGribChangeParID_lev)
+int cdoGribChangeParameterIdentification = -1;
+int cdoGribChangeParID_code;
+int cdoGribChangeParID_ltype;
+int cdoGribChangeParID_lev;
+
+extern void gribapiChangeParameterIdentification(grib_handle *gh, int code, int ltype, int lev); // Function defined in stream_gribapi.c
+
+void streamGrbChangeParameterIdentification(int code, int ltype, int lev)
+{ // NOTE this is a "PROXY" function for gribapiChangeParameterIdentification();
+ // This just sets the globals. There are probably better solutions to this.
+ // The parameter change is done by function gribapiChangeParameterIdentification() in stream_gribapi.c
+ cdoGribChangeParameterIdentification = 1;
+ // Setting this control variable to 1 will cause calling fnc. gribapiChangeParameterIdentification later.
+ // After grib attributes have been changed this variable goes to -1.
+ cdoGribChangeParID_code = code;
+ cdoGribChangeParID_ltype = ltype;
+ cdoGribChangeParID_lev = lev;
+}
+#endif // HIRLAM_EXTENSIONS
int grib1ltypeToZaxisType(int grib_ltype)
{
diff --git a/src/stream_gribapi.c b/src/stream_gribapi.c
index bb7292479..e136d77ed 100644
--- a/src/stream_gribapi.c
+++ b/src/stream_gribapi.c
@@ -35,12 +35,31 @@ 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;
} compvar2_t;
+#ifdef HIRLAM_EXTENSIONS
+// Forward declarations:
+int gribapiGetScanningMode(grib_handle *gh);
+void gribapiSetScanningMode(grib_handle *gh, int scanningMode);
+void convertDataScanningMode(int scanModeOUT, double *data, int gridsize, int iDim, int jDim);
+int gribapiGetUvRelativeToGrid(grib_handle *gh);
+void gribapiSetUvRelativeToGrid(grib_handle *gh, int mode);
+void gribapiChangeParameterIdentification(grib_handle *gh, int code, int ltype, int lev);
+
+int scanModeIN; // global variable
+#endif // HIRLAM_EXTENSIONS
static
int gribapiGetZaxisType(long editionNumber, int grib_ltype)
@@ -107,7 +126,15 @@ int gribapiGetTimeUnits(grib_handle *gh)
int timeunits = -1;
long unitsOfTime = -1;
+#ifndef HIRLAM_EXTENSIONS
grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
+#else
+ int status = grib_get_long(gh, "indicatorOfUnitOfTimeRange", &unitsOfTime);
+ if ( status != 0 )
+ { // We should look for "stepUnits" in GRIB 1 according WMO !
+ status = grib_get_long(gh, "stepUnits", &unitsOfTime);
+ }
+#endif // HIRLAM_EXTENSIONS
GRIB_CHECK(my_grib_set_long(gh, "stepUnits", unitsOfTime), 0);
@@ -130,6 +157,7 @@ void gribapiGetSteps(grib_handle *gh, int timeunits, int *startStep, int *endSte
if ( status == 0 ) *startStep = (int) lpar;
else
{
+ // KNMI-NOTE: We should look for "startStep" in GRIB 1 according WMO !
status = grib_get_long(gh, "startStep", &lpar);
if ( status == 0 )
*startStep = (int) (((double)lpar * timeunit_factor(timeunits, timeunits2)) + 0.5);
@@ -167,7 +195,7 @@ int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
int timeUnits, startStep = 0, endStep;
int tstepRange = 0;
int range;
- long sigofrtime = 3;
+ long sigofrtime = -1;// KNMI-FIX: put it to an undefined state; otherwise for GRIB1 with relative forecast time we always get TRUE for if ( sigofrtime == 3 ) !
if ( gribEditionNumber(gh) > 1 )
{
@@ -181,6 +209,7 @@ int gribapiGetValidityDateTime(grib_handle *gh, int *vdate, int *vtime)
if ( sigofrtime == 3 ) //XXX: This looks like a bug to me, because timeRangeIndicator == 3 does not seem to have the same meaning as significanceOfReferenceTime == 3. I would recommend replacing this condition with `if(!gribapiTimeIsFC())`.
{
gribapiGetDataDateTime(gh, vdate, vtime);
+ // printf("gribapiGetValidityDateTime(): sigofrtime == 3\n");
}
else
{
@@ -1667,6 +1696,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);
}
}
@@ -1762,6 +1793,17 @@ int gribapiDefDateTimeRel(int editionNumber, grib_handle *gh, int rdate, int rti
{
int proDefTempNum = gribapiDefSteptype(editionNumber, gh, productDefinitionTemplate, typeOfGeneratingProcess, tsteptype, gcinit);
+ if ( cdoDebugExt >=100 ) {
+ long productDefinitionTemplateNumber;
+ GRIB_CHECK(grib_get_long(gh, "productDefinitionTemplateNumber", &productDefinitionTemplateNumber), 0);
+ long timeRangeIndicator;
+ GRIB_CHECK(grib_get_long(gh, "timeRangeIndicator", &timeRangeIndicator), 0);
+ long indicatorOfParameter;
+ GRIB_CHECK(grib_get_long(gh, "indicatorOfParameter", &indicatorOfParameter), 0);
+ printf("gribapiDefDateTimeRel()::productDefinitionTemplateNumber = %ld; indicatorOfParameter = %ld; timeRangeIndicator =%ld\n",productDefinitionTemplateNumber, indicatorOfParameter, timeRangeIndicator);
+ printf("gribapiDefDateTimeRel():: proDefTempNum %d,productDefinitionTemplate %d, typeOfGeneratingProcess %d, tsteptype %d)\n",proDefTempNum,productDefinitionTemplate, typeOfGeneratingProcess, tsteptype);
+ }
+
gribapiDefStepUnits(editionNumber, gh, timeunit, proDefTempNum, gcinit);
endStep = (int) ((days*86400.0 + secs)/factor);
@@ -2138,7 +2180,7 @@ void gribapiDefGrid(int editionNumber, grib_handle *gh, int gridID, int comptype
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, "scanningMode", scanflag), 0); // Does this field really exists? NOT according WMO; 128*iScansNegatively + 64*jScansPositively + 32*jPointsAreConsecutive
}
break;
@@ -2718,6 +2760,81 @@ size_t gribapiEncode(int varID, int levelID, int vlistID, int gridID, int zaxisI
GRIB_CHECK(my_grib_set_double(gh, "missingValue", vlistInqVarMissval(vlistID, varID)), 0);
}
+#ifdef HIRLAM_EXTENSIONS
+
+ /*
+ Nj = 550;
+ latitudeOfFirstGridPointInDegrees = -30.8;
+ latitudeOfLastGridPointInDegrees = 24.1;
+ iScansNegatively = 0;
+ jScansPositively = 0;
+ jPointsAreConsecutive = 0;
+ jDirectionIncrementInDegrees = 0.1; */
+
+
+ int iDim=0;
+ int jDim=0;
+ long gridsize;
+ double yfirst, ylast, yinc;
+
+ gridsize = gridInqSize(gridID);
+ iDim = gridInqXsize(gridID);
+ jDim = gridInqYsize(gridID);
+
+ //printf("gribapiEncode(): \n");
+
+ yfirst = gridInqYval(gridID, 0);
+ ylast = gridInqYval(gridID, jDim-1);
+ yinc = gridInqYinc(gridID);
+
+ int scanModeOUT;
+ //scanModeIN = gribapiGetScanningMode(gh);
+ scanModeIN = gridInqScanningMode(gridID); // use global variable scanModeIN
+
+ if (cdoDebugExt>=100)
+ Message("(iDim,jDim) = (%d,%d); scanModeIN=%d; gridsize=%ld", iDim,jDim, scanModeIN, gridsize);
+
+
+ if (cdoGribDataScanningMode != -1) // -1: not used; allowed modes: <0, 64, 96>; Default is 64
+ {
+ scanModeOUT = cdoGribDataScanningMode;
+ convertDataScanningMode(scanModeOUT , data, datasize, iDim, jDim);
+ // This will overrule the old scanning mode of the given grid
+ if (cdoDebugExt>=10) Message("Set GribDataScanningMode (%d) => (%d)", scanModeIN, cdoGribDataScanningMode);
+ gribapiSetScanningMode(gh, cdoGribDataScanningMode);
+
+ if (((scanModeIN==00) && (cdoGribDataScanningMode==64)) ||
+ ((scanModeIN==64) && (cdoGribDataScanningMode==00)) )
+ verticallyFlipGridDefinitionWhenScanningModeChanged(gh, yfirst, ylast, yinc);
+
+ } else
+ {
+ if (cdoDebugExt>=100) Message("Set GribDataScanningMode => (%d) based on used grid", scanModeIN);
+ gribapiSetScanningMode(gh, scanModeIN);
+ }
+
+
+ //else
+ //{ // keep the same scanning mode ..
+ // cdoGribDataScanningMode = scanModeIN;
+ // convertDataScanningMode(scanModeIN , data, gridsize, iDim, jDim);
+ //}
+
+
+ if (cdoGribChangeModeUvRelativeToGrid>=0)
+ // this will overrule/change the UvRelativeToGrid flag;
+ // typically when the wind is rotated with respect to north pole
+ {
+ if (cdoDebugExt>=100) Message("Set ModeUvRelativeToGrid =>%d ( note grid has: %d)", cdoGribChangeModeUvRelativeToGrid, gridInqUvRelativeToGrid(gridID));
+ gribapiSetUvRelativeToGrid(gh,cdoGribChangeModeUvRelativeToGrid);
+ } else
+ {
+ if (cdoDebugExt>=100) Message("Set ModeUvRelativeToGrid =>%d based on used grid", gridInqUvRelativeToGrid(gridID));
+ gribapiSetUvRelativeToGrid(gh,gridInqUvRelativeToGrid(gridID));
+ }
+#endif //HIRLAM_EXTENSIONS
+
+ //Message("About to call grib_set_double_array( data=%x, datasize=%ld)..", data, datasize);
GRIB_CHECK(grib_set_double_array(gh, "values", data, (size_t)datasize), 0);
/* get the size of coded message */
@@ -2739,6 +2856,387 @@ size_t gribapiEncode(int varID, int levelID, int vlistID, int gridID, int zaxisI
}
#endif
+#ifdef HIRLAM_EXTENSIONS
+
+#if defined (HAVE_LIBGRIB_API)
+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 (cdoDebugExt>=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 (cdoDebugExt>=30)
+ {
+ long paramId, levelTypeId, levelId, uvRelativeToGrid;
+ GRIB_CHECK(grib_get_long(gh, "uvRelativeToGrid", &uvRelativeToGrid), 0);
+ GRIB_CHECK(grib_get_long(gh, "indicatorOfParameter", ¶mId), 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);
+
+}
+
+int gribapiGetUvRelativeToGrid(grib_handle *gh)
+{
+ long uvRelativeToGridMode;
+
+ GRIB_CHECK(my_grib_set_long(gh, "uvRelativeToGrid", &uvRelativeToGridMode), 0);
+ if (cdoDebugExt>=20)
+ printf("gribapiGetUvRelativeToGrid(): uvRelativeToGrid mode = %02d ; \n",(int)uvRelativeToGridMode);
+
+ return uvRelativeToGridMode;
+}
+
+void gribapiChangeParameterIdentification(grib_handle *gh, int code, int ltype, int lev)
+{
+ long indicatorOfParameter, indicatorOfTypeOfLevel, level; // timeRangeIndicator: could be included later
+ indicatorOfParameter = code;
+ indicatorOfTypeOfLevel = ltype;
+ level = lev;
+
+ if (indicatorOfParameter!=-1)
+ GRIB_CHECK(my_grib_set_long(gh, "indicatorOfParameter", indicatorOfParameter), 0);
+ if (indicatorOfTypeOfLevel!=-1)
+ GRIB_CHECK(my_grib_set_long(gh, "indicatorOfTypeOfLevel", indicatorOfTypeOfLevel), 0);
+ if (level!=-1)
+ GRIB_CHECK(my_grib_set_long(gh, "level", level), 0);
+}
+
+void gribapiSetUvRelativeToGrid(grib_handle *gh, int mode)
+{
+ long uvRelativeToGridMode = mode;
+ long uvRelativeToGridModeOld;
+
+ GRIB_CHECK(grib_get_long(gh, "uvRelativeToGrid", &uvRelativeToGridModeOld), 0);
+
+ if (cdoDebugExt>=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 (cdoDebugExt>=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 (cdoDebugExt>=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);
+ }
+}
+
+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 (cdoDebugExt>=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 (cdoDebugExt>=30) printf("convertDataScanningMode(): ERROR: (iDim*jDim)!= gridsize; (%d * %d) != %d\n", iDim,jDim, gridsize);
+ return;
+ }
+ if (cdoDebugExt>=30) printf("convertDataScanningMode(): scanModeIN=%02d => scanModeOUT=%02d ; where: (iDim * jDim == gridsize) (%d*%d == %d)\n",scanModeIN, scanModeOUT, iDim,jDim, gridsize);
+
+ if (cdoDebugExt>=100)
+ {
+ printf("convertDataScanningMode(): data IN:\n");
+ for (j=0; j<jDim; j++)
+ {
+ int jXiDim = j*iDim;
+ for (i=0; i<iDim; i++)
+ {
+ idxIN = i + jXiDim;
+ printf("%03.0f, ",data[idxIN]);
+ }
+ printf("\n");
+ }
+ }
+
+ if (scanModeIN==scanModeOUT)
+ {
+ if (cdoDebugExt>=30) printf("convertDataScanningMode(): INFO: Nothing to do; scanModeIN==scanModeOUT..\n");
+ return;
+ }
+
+ if (0)
+ {
+ return;
+ if (scanModeOUT==00)
+ {
+ if (cdoDebugExt>0) printf("convertDataScanningMode(): Leave data unchaged BUT set scanModeOUT=00.\n");
+ // CHECK: Looks like that GRIB-API provide (no matter what) data in the scannning mode 00, even it is store in the gribfile as 64 !!
+ return;
+ }
+ }
+ double *dataCopy = NULL;
+ dataCopy = (double *) malloc(gridsize*sizeof(double));
+
+ memcpy((void*)dataCopy,(void*) data, gridsize*sizeof(double));
+
+ if (scanModeIN==64) // Scanning Mode (00 dec) +i, -j; i direction consecutive (row-major order West->East & South->North )
+ { // Scanning Mode (64 dec) +i, +j; i direction consecutive (row-major order West->East & North->South )
+ // Scanning Mode (96 dec) +i, +j; j direction consecutive (column-major order North->South & West->East )
+ if (scanModeOUT==00)
+ // CHECK: Looks like that GRIB-API provide (no matter what) data in the scannning mode 00, even it is store in the gribfile as 64 !!
+#define VERTICAL_FLIP
+#ifdef VERTICAL_FLIP
+ { // flip the data vertically ..
+ idxIN= 0; idxOUT= (jDim-1)*iDim;
+ if (cdoDebugExt>=30) printf("convertDataScanningMode(): copying rows nr. (%04d : %04d)\n",0,jDim-1);
+ for (j=0; j<jDim; j++)
+ {
+ memcpy((void*)&data[idxOUT], (void*)&dataCopy[idxIN], iDim*sizeof(double));
+ idxIN += iDim; idxOUT -= iDim;
+ }
+ } // end if (scanModeOUT==00)*/
+#endif
+#ifdef HORIZONTAL_FLIP
+ { // flip data horizontally ...
+ if (1)
+ {
+ if (cdoDebugExt>=30) printf("convertDataScanningMode(): copying columns nr. (%04d : %04d);\n", 0, iDim-1);
+ for (i=0; i<iDim; i++)
+ {
+ for (j=0; j<jDim; j++)
+ {
+ int jXiDim = j*iDim;
+ idxIN = i + jXiDim;
+ //data[idxIN] = (double) (100.0*j +i); // just some testdata ..
+ idxOUT = iDim - i -1 + jXiDim;
+ //printf("[%03d=>%03d] = %f;",idxIN,idxOUT,dataCopy[idxIN]);
+ data[idxOUT] = dataCopy[idxIN];
+ }
+ }
+ }
+ } // end if (scanModeOUT==00)
+#endif
+
+ if (scanModeOUT==96)
+ { // transpose the data
+ if (cdoDebugExt>=30) printf("convertDataScanningMode(): transpose data rows=>columns nr. (%04d : %04d) => (%04d : %04d);\n", 0, iDim-1, 0, jDim-1);
+ for (j=0; j<jDim; j++)
+ {
+ int jXiDim = j*iDim;
+ for (i=0; i<iDim; i++)
+ {
+ idxIN = i + jXiDim;
+ idxOUT = j + i*jDim;
+ //printf("[%03d=>%03d] = %f;",idxIN,idxOUT,dataCopy[idxIN]);
+ data[idxOUT] = dataCopy[idxIN];
+ }
+ //printf(".\n");
+ }
+ } // end if (scanModeOUT==96)
+ } // end if (scanModeIN==64)
+
+ if (scanModeIN==00) // Scanning Mode (00 dec) +i, -j; i direction consecutive (row-major order West->East & South->North )
+ { // Scanning Mode (64 dec) +i, +j; i direction consecutive (row-major order West->East & North->South )
+ // Scanning Mode (96 dec) +i, +j; j direction consecutive (column-major order North->South & West->East )
+ if (scanModeOUT==64)
+ { // flip the data vertically ..
+ idxIN= 0; idxOUT= (jDim-1)*iDim;
+ for (j=0; j<jDim; j++)
+ {
+ if (cdoDebugExt>=25) printf("convertDataScanningMode(): copying row nr. %04d; [idxIN=%08d] => [idxOUT=%08d]\n",j, idxIN, idxOUT);
+ memcpy((void*)&data[idxOUT], (void*)&dataCopy[idxIN], iDim*sizeof(double));
+ idxIN += iDim; idxOUT -= iDim;
+ }
+ } // end if (scanModeOUT==64)
+
+ if (scanModeOUT==96)
+ { // transpose the data
+ int jInv;
+ for (j=0; j<jDim; j++)
+ {
+ if (cdoDebugExt>=30) printf("convertDataScanningMode(): processing row nr. %04d;\n", j);
+ jInv = (jDim-1) -j;
+ for (i=0; i<iDim; i++)
+ data[j + i*jDim] = dataCopy[i + jInv*iDim]; // source data has -j
+ }
+ } // end if (scanModeOUT==96)
+ } // end if (scanModeIN==00)
+
+ if (scanModeIN==96) // Scanning Mode (00 dec) +i, -j; i direction consecutive (row-major order West->East & South->North )
+ { // Scanning Mode (64 dec) +i, +j; i direction consecutive (row-major order West->East & North->South )
+ // Scanning Mode (96 dec) +i, +j; j direction consecutive (column-major order North->South & West->East )
+ if (scanModeOUT==64)
+ { // transpose the data
+ for (j=0; j<jDim; j++)
+ {
+ if (cdoDebugExt>=30) printf("convertDataScanningMode(): processing row nr. %04d;\n", j);
+ int jXiDim = j*iDim;
+ for (i=0; i<iDim; i++)
+ //data[j + i*jDim] = dataCopy[i + j*iDim];
+ data[i + jXiDim] = dataCopy[j + i*jDim];
+ }
+ } // end if (scanModeOUT==64)
+
+ if (scanModeOUT==00)
+ { // transpose the data
+ idxIN= 0; idxOUT= 0;
+ int jInv;
+ for (j=0; j<jDim; j++)
+ {
+ if (cdoDebugExt>=30) printf("convertDataScanningMode(): processing row nr. %04d;\n", j);
+ jInv = (jDim-1) -j;
+ int jXiDim = j*iDim;
+ for (i=0; i<iDim; i++)
+ //data[jInv + iXjDim] = dataCopy[i + jXiDim]; // target data has -j
+ data[i + jXiDim] = dataCopy[jInv + i*jDim]; // target data has -j
+ }
+ } // end if (scanModeOUT==00)
+ } // end if (scanModeIN==96)
+
+ if (cdoDebugExt>=100)
+ {
+ printf("convertDataScanningMode(): data OUT (new scanning mode):\n");
+ for (j=0; j<jDim; j++)
+ {
+ int jXiDim = j*iDim;
+ for (i=0; i<iDim; i++)
+ {
+ idxIN = i + jXiDim;
+ printf("%03.0f, ",data[idxIN]);
+ }
+ printf("\n");
+ }
+ }
+
+ free(dataCopy); return;
+}
+#endif // HAVE_LIBGRIB_API
+
+#endif // HIRLAM_EXTENSIONS
+
+
/*
* Local Variables:
* c-file-style: "Java"
diff --git a/src/varscan.c b/src/varscan.c
index d89fa3297..72d4ac354 100644
--- a/src/varscan.c
+++ b/src/varscan.c
@@ -579,6 +579,7 @@ void cdi_generate_vars(stream_t *streamptr)
double level_sf = 1;
int vlistID = streamptr->vlistID;
+// printf("cdi_generate_vars() ..\n"); // KNMI-debug-prints
int *varids = (int *) Malloc(nvars*sizeof(int));
for ( unsigned varID = 0; varID < nvars; varID++ ) varids[varID] = (int)varID;
@@ -626,6 +627,8 @@ void cdi_generate_vars(stream_t *streamptr)
zaxisID = UNDEFID;
+// printf("Processing (index=%d): param.nr. %d; gridID=%d, zaxistype=%d \n", index, param, gridID, zaxistype); // KNMI-debug-prints
+
/* consistency check: test if all subtypes have the same levels: */
unsigned nlevels = vartable[varid].recordTable[0].nlevels;
for (int isub=1; isub<vartable[varid].nsubtypes; isub++) {
@@ -727,6 +730,7 @@ void cdi_generate_vars(stream_t *streamptr)
}
const char *unitptr = cdiUnitNamePtr(vartable[varid].level_unit);
+ //Message("calling varDefZaxis(zaxistype: %d, nlevels: %d) ", zaxistype, nlevels); // KNMI-debug-prints
zaxisID = varDefZaxis(vlistID, zaxistype, (int)nlevels, dlevels, lbounds, dlevels1, dlevels2,
(int)Vctsize, Vct, NULL, NULL, unitptr, 0, 0, ltype1);
@@ -745,6 +749,7 @@ void cdi_generate_vars(stream_t *streamptr)
if ( lbounds ) Free(dlevels1);
if ( lbounds ) Free(dlevels2);
Free(dlevels);
+// printf("streamNewVar(streamID, gridID=%d, zaxisID=%d)\n",gridID, zaxisID); // KNMI-debug-prints
/* define new subtype for tile set */
int tilesetID = CDI_UNDEFID;
@@ -1019,6 +1024,7 @@ int varDefZaxis(int vlistID, int zaxistype, int nlevels, const double *levels, i
{
if ( ! zaxisglobdefined )
{
+ //Message("calling zaxisCreate(zaxistype: %d, nlevels: %d) ", zaxistype, nlevels); // KNMI-debug-prints
zaxisID = zaxisCreate(zaxistype, nlevels);
zaxisDefLevels(zaxisID, levels);
if ( lbounds )
diff --git a/src/vlist.c b/src/vlist.c
index f94141399..c4ab9e88d 100644
--- a/src/vlist.c
+++ b/src/vlist.c
@@ -489,6 +489,7 @@ int vlist_generate_zaxis(int vlistID, int zaxistype, int nlevels, const double *
{
if ( ! zaxisglobdefined )
{
+ //Message("calling zaxisCreate(zaxistype: %d, nlevels: %d) ", zaxistype, nlevels); // KNMI-debug-prints
zaxisID = zaxisCreate(zaxistype, nlevels);
zaxisDefLevels(zaxisID, levels);
if ( has_bounds )
--
GitLab