Added support for polar stereographic projection.

ChangeLog

@@ -4,6 +4,7 @@
 
 2018-09-13  Uwe Schulzweida
 
+	* Added support for polar stereographic projection
 	* grib2ScaleFactor: added support for negative scale factor (bug fix)
 
 2018-09-12  Uwe Schulzweida

src/grid.c

@@ -3932,6 +3932,34 @@ int gridVerifyGribParamLCC(double missval, double *lon_0, double *lat_0, double
   return 0;
 }
 
+
+int gridVerifyGribParamSTERE(double missval, double *lon_0, double *lat_ts, double *lat_0,
+                             double *a, double *xval_0, double *yval_0, double *x_0, double *y_0)
+{
+  static bool lwarn = true;
+
+  if ( lwarn )
+    {
+      // lwarn = false;
+      const char *projection = "lambert_conformal_conic";
+      if ( IS_EQUAL(*lon_0, missval) ) { Warning("%s mapping parameter %s missing!", projection, "straight_vertical_longitude_from_pole"); }
+      if ( IS_EQUAL(*lat_0, missval) ) { Warning("%s mapping parameter %s missing!", projection, "latitude_of_projection_origin"); }
+      if ( IS_EQUAL(*lat_ts, missval) ) { Warning("%s mapping parameter %s missing!", projection, "standard_parallel"); }
+      if ( IS_NOT_EQUAL(*x_0, missval) && IS_NOT_EQUAL(*y_0, missval) && (IS_EQUAL(*xval_0, missval) || IS_EQUAL(*yval_0, missval)) )
+        {
+          if ( proj_stere_to_lonlat_func )
+            {
+              *xval_0 = -(*x_0); *yval_0 = -(*y_0);
+              proj_stere_to_lonlat_func(missval, *lon_0, *lat_ts, *lat_0, *a, 0.0, 0.0, (size_t)1, xval_0, yval_0);
+            }
+          if ( IS_EQUAL(*xval_0, missval) || IS_EQUAL(*yval_0, missval) )
+            Warning("%s mapping parameter %s missing!", projection, "longitudeOfFirstGridPointInDegrees and latitudeOfFirstGridPointInDegrees");
+        }
+    }
+
+  return 0;
+}
+
 /*
 @Function  gridDefParamSTERE
 @Title     Define the parameter of a Polar stereographic grid

src/grid.h

@@ -166,6 +166,8 @@ struct addIfNewRes cdiVlistAddGridIfNew(int vlistID, grid_t *grid, int mode);
 
 int gridVerifyGribParamLCC(double missval, double *lon_0, double *lat_0, double *lat_1, double *lat_2,
                            double *a, double *rf, double *xval_0, double *yval_0, double *x_0, double *y_0);
+int gridVerifyGribParamSTERE(double missval, double *lon_0, double *lat_ts, double *lat_0,
+                             double *a, double *xval_0, double *yval_0, double *x_0, double *y_0);
 
 #endif
 /*

src/stream_gribapi.c

@@ -279,16 +279,25 @@ double grib2ScaleFactor(long factor)
   switch(factor)
     {
       case GRIB_MISSING_LONG: return 1;
-      case 0: return 1;
-      case 1: return 0.1;
-      case 2: return 0.01;
-      case 3: return 0.001;
-      case 4: return 0.0001;
-      case 5: return 0.00001;
-      case 6: return 0.000001;
-      case 7: return 0.0000001;
-      case 8: return 0.00000001;
-      case 9: return 0.000000001;
+      case -9: return 1000000000;
+      case -8: return 100000000;
+      case -7: return 10000000;
+      case -6: return 1000000;
+      case -5: return 100000;
+      case -4: return 10000;
+      case -3: return 1000;
+      case -2: return 100;
+      case -1: return 10;
+      case  0: return 1;
+      case  1: return 0.1;
+      case  2: return 0.01;
+      case  3: return 0.001;
+      case  4: return 0.0001;
+      case  5: return 0.00001;
+      case  6: return 0.000001;
+      case  7: return 0.0000001;
+      case  8: return 0.00000001;
+      case  9: return 0.000000001;
       default: return 0;
     }
 }
@@ -318,8 +327,6 @@ void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lboun
   int status = grib_get_long(gh, "typeOfFirstFixedSurface", &lpar); //1 byte
   if ( status == 0 )
     {
-      long llevel;
-
       *leveltype1 = (int) lpar;
 
       status = grib_get_long(gh, "typeOfSecondFixedSurface", &lpar); //1 byte
@@ -349,7 +356,7 @@ void grib2GetLevel(grib_handle *gh, int *leveltype1, int *leveltype2, int *lboun
             break;
         }
 
-      long factor;
+      long factor, llevel;
       GRIB_CHECK(grib_get_long(gh, "scaleFactorOfFirstFixedSurface", &factor), 0);      //1 byte
       GRIB_CHECK(grib_get_long(gh, "scaledValueOfFirstFixedSurface", &llevel), 0);      //4 byte
       *level1 = calcLevel(*level_sf, factor, llevel);
@@ -560,7 +567,7 @@ void gribapiDefProjSTERE(grib_handle *gh, int gridID)
   if ( proj_lonlat_to_stere_func )
     {
       x_0 = xval_0; y_0 = yval_0;
-      // proj_lonlat_to_stere_func(grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, (size_t)1, &x_0, &y_0);
+      proj_lonlat_to_stere_func(grid_missval, lon_0, lat_ts, lat_0, a, (size_t)1, &x_0, &y_0);
       if ( IS_NOT_EQUAL(x_0, grid_missval) && IS_NOT_EQUAL(y_0, grid_missval) )
         { x_0 = -x_0; y_0 = -y_0; }
     }
@@ -2123,7 +2130,7 @@ void gribapiDefGridSTERE(grib_handle *gh, int editionNumber, int gridID)
 
   double lon_0, lat_ts, lat_0, a, xval_0, yval_0, x_0, y_0;
   gridInqParamSTERE(gridID, grid_missval, &lon_0, &lat_ts, &lat_0, &a, &xval_0, &yval_0, &x_0, &y_0);
-  // gridVerifyGribParamSTERE(grid_missval, &lon_0, &lat_ts, &lat_0, &a, &xval_0, &yval_0, &x_0, &y_0);
+  gridVerifyGribParamSTERE(grid_missval, &lon_0, &lat_ts, &lat_0, &a, &xval_0, &yval_0, &x_0, &y_0);
   if ( xval_0 < 0 ) xval_0 += 360;
   int projflag = 0;