GRIB rotated grids: invert angle of rotation (north to south pole) (bug fix)

ChangeLog

@@ -2,6 +2,10 @@
 
 	* Version 1.7.0 released
 
+2015-09-10  Uwe Schulzweida
+
+	* GRIB rotated grids: invert angle of rotation (north to south pole) (bug fix)
+
 2015-08-18  Uwe Schulzweida
 
 	* added support for netCDF Scalar Coordinate Variables

src/gribapi_utilities.c

@@ -35,7 +35,7 @@ char* gribCopyString(grib_handle* gribHandle, const char* key)
     char* result = (char *)xmalloc(length);
     if(!grib_get_string(gribHandle, key, result, &length))
     result = (char *)xrealloc(result, length);
-  
+
     else
       {
           free(result);
@@ -762,8 +762,9 @@ void gribapiGetGrid(grib_handle *gh, grid_t *grid)
       FAIL_ON_GRIB_ERROR(grib_get_double, gh, "longitudeOfSouthernPoleInDegrees", &grid->xpole);
       FAIL_ON_GRIB_ERROR(grib_get_double, gh, "angleOfRotation", &grid->angle);
       /* change from south to north pole */
-      grid->ypole = -grid->ypole;
+      if ( fabs(grid->ypole) > 0 ) grid->ypole = -grid->ypole;
       grid->xpole =  grid->xpole - 180;
+      if ( fabs(grid->angle) > 0 ) grid->angle = -grid->angle;
     }
 
   grid->xvals = NULL;

src/grid.c

@@ -1597,6 +1597,7 @@ double gridInqYinc(int gridID)
 */
 double gridInqXpole(int gridID)
 {
+  // Xpole -> grid_north_pole_longitude
   grid_t *gridptr = gridID2Ptr(gridID);
 
   return gridptr->xpole;
@@ -1614,6 +1615,7 @@ double gridInqXpole(int gridID)
 */
 void gridDefXpole(int gridID, double xpole)
 {
+  // Xpole -> grid_north_pole_longitude
   grid_t *gridptr = gridID2Ptr(gridID);
 
   if ( memcmp(gridptr->xstdname, "grid", 4) != 0 )
@@ -1639,6 +1641,7 @@ void gridDefXpole(int gridID, double xpole)
 */
 double gridInqYpole(int gridID)
 {
+  // Ypole -> grid_north_pole_latitude
   grid_t *gridptr = gridID2Ptr(gridID);
 
   return (gridptr->ypole);
@@ -1656,6 +1659,7 @@ double gridInqYpole(int gridID)
 */
 void gridDefYpole(int gridID, double ypole)
 {
+  // Ypole -> grid_north_pole_latitude
   grid_t *gridptr = gridID2Ptr(gridID);
 
   if ( memcmp(gridptr->ystdname, "grid", 4) != 0 )
@@ -1681,6 +1685,7 @@ void gridDefYpole(int gridID, double ypole)
 */
 double gridInqAngle(int gridID)
 {
+  // Angle -> north_pole_grid_longitude
   grid_t *gridptr = gridID2Ptr(gridID);
 
   return (gridptr->angle);
@@ -1698,6 +1703,7 @@ double gridInqAngle(int gridID)
 */
 void gridDefAngle(int gridID, double angle)
 {
+  // Angle -> north_pole_grid_longitude
   grid_t *gridptr = gridID2Ptr(gridID);
 
   if ( gridptr->isRotated != TRUE || IS_NOT_EQUAL(gridptr->angle, angle) )

src/stream_cgribex.c

@@ -349,7 +349,7 @@ void cgribexGetGrid(stream_t *streamptr, int *isec2, double *fsec2, int *isec4,
       grid->isRotated = TRUE;
       grid->ypole     = - ISEC2_LatSP*0.001;
       grid->xpole     =   ISEC2_LonSP*0.001 - 180;
-      grid->angle     = FSEC2_RotAngle;
+      grid->angle     = - FSEC2_RotAngle;
     }
 
   grid->xvals = NULL;
@@ -1762,7 +1762,9 @@ void cgribexDefGrid(int *isec1, int *isec2, double *fsec2, int *isec4, int gridI
 	  {
 	    ISEC2_LatSP = - (int)lround(gridInqYpole(gridID) * 1000);
 	    ISEC2_LonSP =   (int)lround((gridInqXpole(gridID) + 180) * 1000);
-            FSEC2_RotAngle = gridInqAngle(gridID);
+            double angle = gridInqAngle(gridID);
+            if ( fabs(angle) > 0 ) angle = -angle;
+            FSEC2_RotAngle = angle;
 	  }
 
 	/* East -> West */

src/stream_gribapi.c

@@ -2037,8 +2037,9 @@ void gribapiDefGrid(int editionNumber, grib_handle *gh, int gridID, int comptype
 	    ypole = gridInqYpole(gridID);
 	    angle = gridInqAngle(gridID);
 	    /* change from north to south pole */
-	    ypole = -ypole;
+	    if ( fabs(ypole) > 0 ) ypole = -ypole;
 	    xpole =  xpole + 180;
+            if ( fabs(angle) > 0 ) angle = -angle;
 	    GRIB_CHECK(my_grib_set_double(gh, "latitudeOfSouthernPoleInDegrees",  ypole), 0);
 	    GRIB_CHECK(my_grib_set_double(gh, "longitudeOfSouthernPoleInDegrees", xpole), 0);
 	    GRIB_CHECK(my_grib_set_double(gh, "angleOfRotation", angle), 0);