From e84f98ed46e898eb4f8c97415b3bf7d5d2ae2feb Mon Sep 17 00:00:00 2001
From: Uwe Schulzweida <uwe.schulzweida@mpimet.mpg.de>
Date: Mon, 5 Jun 2023 21:10:42 +0200
Subject: [PATCH] gribapi LCC: add support for parameter xpole/ypole (bug fix)
---
ChangeLog | 1 +
src/cdi.h | 2 ++
src/grid.c | 9 +++++++++
src/stream_gribapi.c | 5 +++++
4 files changed, 17 insertions(+)
diff --git a/ChangeLog b/ChangeLog
index cf2b3830f..4255cf751 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,6 @@
2023-06-04 Uwe Schulzweida
+ * gribapi LCC: add support for parameter xpole/ypole (bug fix)
* gribapiDefGridLCC: add LaDInDegrees for GRIB2 (bug fix)
* gribapiDefGridLCC: store DxInMetres/DyInMetres as double (bug fix)
diff --git a/src/cdi.h b/src/cdi.h
index 4c48f4edd..1cf74eab1 100644
--- a/src/cdi.h
+++ b/src/cdi.h
@@ -1401,6 +1401,8 @@ struct CDI_GridProjParams
double yval_0; // Latitude of the first grid point in degree (optional)
double x_0; // False easting (optional)
double y_0; // False northing (optional)
+ double xpole; // Longitude of southern pole
+ double ypole; // Latitude of southern pole
};
void gridProjParamsInit(struct CDI_GridProjParams *gridProjParams);
diff --git a/src/grid.c b/src/grid.c
index 9ac8cf138..3d417e0de 100644
--- a/src/grid.c
+++ b/src/grid.c
@@ -3651,6 +3651,9 @@ gridProjParamsInit(struct CDI_GridProjParams *gpp)
gpp->yval_0 = CDI_Grid_Missval; // Latitude of the first grid point in degree (optional)
gpp->x_0 = CDI_Grid_Missval; // False easting (optional)
gpp->y_0 = CDI_Grid_Missval; // False northing (optional)
+ gpp->xpole = CDI_Grid_Missval; // Longitude of southern pole
+ gpp->ypole = CDI_Grid_Missval; // Latitude of southern pole
+
// clang-format on
}
@@ -3676,6 +3679,10 @@ gridDefParamsCommon(int gridID, struct CDI_GridProjParams gpp)
cdiDefAttFlt(gridID, CDI_GLOBAL, "longitudeOfFirstGridPointInDegrees", CDI_DATATYPE_FLT64, 1, &gpp.xval_0);
if (IS_NOT_EQUAL(gpp.yval_0, gpp.mv))
cdiDefAttFlt(gridID, CDI_GLOBAL, "latitudeOfFirstGridPointInDegrees", CDI_DATATYPE_FLT64, 1, &gpp.yval_0);
+ if (IS_NOT_EQUAL(gpp.xpole, gpp.mv))
+ cdiDefAttFlt(gridID, CDI_GLOBAL, "longitudeOfSouthernPoleInDegrees", CDI_DATATYPE_FLT64, 1, &gpp.xpole);
+ if (IS_NOT_EQUAL(gpp.ypole, gpp.mv))
+ cdiDefAttFlt(gridID, CDI_GLOBAL, "latitudeOfSouthernPoleInDegrees", CDI_DATATYPE_FLT64, 1, &gpp.ypole);
}
/*
@@ -3774,6 +3781,8 @@ gridInqParamsLCC(int gridID, struct CDI_GridProjParams *gpp)
else if (str_is_equal(attname, "false_northing")) gpp->y_0 = attflt[0];
else if (str_is_equal(attname, "longitudeOfFirstGridPointInDegrees")) gpp->xval_0 = attflt[0];
else if (str_is_equal(attname, "latitudeOfFirstGridPointInDegrees")) gpp->yval_0 = attflt[0];
+ else if (str_is_equal(attname, "longitudeOfSouthernPoleInDegrees")) gpp->xpole = attflt[0];
+ else if (str_is_equal(attname, "latitudeOfSouthernPoleInDegrees")) gpp->ypole = attflt[0];
else if (str_is_equal(attname, "standard_parallel"))
{
gpp->lat_1 = attflt[0];
diff --git a/src/stream_gribapi.c b/src/stream_gribapi.c
index ef41f4473..cdc4cc736 100644
--- a/src/stream_gribapi.c
+++ b/src/stream_gribapi.c
@@ -612,6 +612,8 @@ gribapiDefProjLCC(grib_handle *gh, int gridID)
long projflag = 0;
grib_get_double(gh, "longitudeOfFirstGridPointInDegrees", &gpp.xval_0);
grib_get_double(gh, "latitudeOfFirstGridPointInDegrees", &gpp.yval_0);
+ grib_get_double(gh, "longitudeOfSouthernPoleInDegrees", &gpp.xpole);
+ grib_get_double(gh, "latitudeOfSouthernPoleInDegrees", &gpp.ypole);
grib_get_double(gh, "LoVInDegrees", &gpp.lon_0);
grib_get_double(gh, "Latin1InDegrees", &gpp.lat_1);
grib_get_double(gh, "Latin2InDegrees", &gpp.lat_2);
@@ -2268,6 +2270,9 @@ gribapiDefGridLCC(grib_handle *gh, int editionNumber, int gridID, int uvRelative
GRIB_CHECK(my_grib_set_double(gh, "Latin2InDegrees", gpp.lat_2), 0);
GRIB_CHECK(my_grib_set_long(gh, "projectionCentreFlag", projflag), 0);
+ if (gpp.xpole >= -180 && gpp.xpole <= 360) GRIB_CHECK(my_grib_set_double(gh, "longitudeOfSouthernPoleInDegrees", gpp.xpole), 0);
+ if (gpp.ypole >= -90 && gpp.ypole <= 90) GRIB_CHECK(my_grib_set_double(gh, "latitudeOfSouthernPoleInDegrees", gpp.ypole), 0);
+
long shapeOfTheEarth = radiusToShapeOfTheEarth(gpp.a);
if (shapeOfTheEarth) GRIB_CHECK(my_grib_set_long(gh, "shapeOfTheEarth", shapeOfTheEarth), 0);
--
GitLab