Prepare gridDefParamLCC() for refactoing.

src/grid.c

@@ -114,6 +114,8 @@ void grid_init(grid_t *gridptr)
   gridptr->lcc.lat2      = 0.0;
   gridptr->lcc.xinc      = 0.0;
   gridptr->lcc.yinc      = 0.0;
+  gridptr->lcc.a         = 0.0;
+  gridptr->lcc.rf        = 0.0;
   gridptr->lcc.projflag  = 0;
   gridptr->lcc.scanflag  = 0;
   gridptr->lcc.defined   = FALSE;
@@ -1009,6 +1011,13 @@ void gridVerifyProj(int gridID)
       gridSetName(gridptr->x.units, "degrees");
       gridSetName(gridptr->y.units, "degrees");
     }
+  else if ( projtype == CDI_PROJ_LCC )
+    {
+      gridptr->x.stdname = xystdname_tab[grid_xystdname_projection][0];
+      gridptr->y.stdname = xystdname_tab[grid_xystdname_projection][1];
+      gridSetName(gridptr->x.units, "m");
+      gridSetName(gridptr->y.units, "m");
+    }
 }
 
 /*
@@ -2474,9 +2483,7 @@ void gridComplete(grid_t *grid)
 	    double *xvals = (double *) Malloc((size_t)grid->x.size * sizeof (double));
 	    gridGenXvals(grid->x.size, grid->x.first, grid->x.last, grid->x.inc, xvals);
 	    grid->x.vals = xvals;
-	    /*
-	    gridDefXinc(gridID, grid->x.inc);
-	    */
+	    // gridDefXinc(gridID, grid->x.inc);
 	  }
 
 	if ( grid->y.flag == 2 )
@@ -2485,9 +2492,7 @@ void gridComplete(grid_t *grid)
 	    double *yvals = (double *) Malloc((size_t)grid->y.size * sizeof (double));
 	    gridGenYvals(gridtype, grid->y.size, grid->y.first, grid->y.last, grid->y.inc, yvals);
 	    grid->y.vals = yvals;
-	    /*
-	    gridDefYinc(gridID, grid->y.inc);
-	    */
+	    // gridDefYinc(gridID, grid->y.inc);
 	  }
 
 	if ( grid->projtype == CDI_PROJ_RLL )
@@ -3667,6 +3672,32 @@ void gridDefParamLCC(int gridID, double originLon, double originLat, double lonP
       gridptr->lcc.defined   = TRUE;
       gridMark4Update(gridID);
     }
+
+#ifdef TEST_LCC
+  cdiGridDefKeyStr(gridID, CDI_KEY_MAPPING, CDI_MAX_NAME, "Lambert_Conformal");
+
+  const char *mapping = "lambert_conformal_conic";
+  cdiGridDefKeyStr(gridID, CDI_KEY_MAPNAME, CDI_MAX_NAME, mapping);
+  cdiDefAttTxt(gridID, CDI_GLOBAL, "grid_mapping_name", strlen(mapping), mapping);
+  int nlats = 0;
+  double lats[2];
+  lats[nlats++] = lat1;
+  if ( IS_NOT_EQUAL(lat1, lat2) ) lats[nlats++] = lat2;
+  cdiDefAttFlt(gridID, CDI_GLOBAL, "standard_parallel", CDI_DATATYPE_FLT64, nlats, lats);
+  cdiDefAttFlt(gridID, CDI_GLOBAL, "longitude_of_central_meridian", CDI_DATATYPE_FLT64, 1, &lonParY);
+  cdiDefAttFlt(gridID, CDI_GLOBAL, "latitude_of_projection_origin", CDI_DATATYPE_FLT64, 1, &lat2);
+  double a = 0;
+  if ( a > 0 ) cdiDefAttFlt(gridID, CDI_GLOBAL, "earth_radius", CDI_DATATYPE_FLT64, 1, &a);
+  double rf = 0;
+  if ( rf > 0 ) cdiDefAttFlt(gridID, CDI_GLOBAL, "inverse_flattening", CDI_DATATYPE_FLT64, 1, &rf);
+  cdiDefAttFlt(gridID, CDI_GLOBAL, "longitudeOfFirstGridPointInDegrees", CDI_DATATYPE_FLT64, 1, &originLon);
+  cdiDefAttFlt(gridID, CDI_GLOBAL, "latitudeOfFirstGridPointInDegrees", CDI_DATATYPE_FLT64, 1, &originLat);
+
+  //grid_t *gridptr = grid_to_pointer(gridID);
+  gridptr->projtype = CDI_PROJ_LCC;
+
+  gridVerifyProj(gridID);
+#endif
 }
 
 /*

src/grid.h

@@ -67,6 +67,8 @@ struct grid_lcc_t {
   double  lat2;
   double  xinc;
   double  yinc;
+  double  a;
+  double  rf;
   int     projflag;
   short   scanflag;             /* must be combination of 128, 64, 32 */
   short   defined;

src/stream_cgribex.c

@@ -41,12 +41,7 @@ int cgribexGetGridType(int *isec2)
     case  GRIB1_GTYPE_LATLON:     { gridtype = GRID_LONLAT;     break; }
     case  GRIB1_GTYPE_LATLON_ROT: { gridtype = GRID_PROJECTION; break; }
     case  GRIB1_GTYPE_LCC:        { gridtype = GRID_LCC;        break; }
-    case  GRIB1_GTYPE_GAUSSIAN:   { if ( ISEC2_Reduced )
-	                              gridtype = GRID_GAUSSIAN_REDUCED;
-                         	    else
-				      gridtype = GRID_GAUSSIAN;
-          	                    break;
-                                  }
+    case  GRIB1_GTYPE_GAUSSIAN:   { gridtype = ISEC2_Reduced ? GRID_GAUSSIAN_REDUCED : GRID_GAUSSIAN; break; }
     case  GRIB1_GTYPE_SPECTRAL:   { gridtype = GRID_SPECTRAL;   break; }
     case  GRIB1_GTYPE_GME:        { gridtype = GRID_GME;        break; }
     }
@@ -287,7 +282,7 @@ void cgribexGetGrid(stream_t *streamptr, int *isec2, int *isec4, grid_t *grid, i
         Error("numberOfPoints (%d) and gridSize (%d) differ!",
               ISEC4_NumValues, ISEC2_NumLon*ISEC2_NumLat);
 
-      grid->size  = ISEC4_NumValues;
+      grid->size   = ISEC4_NumValues;
       grid->x.size = ISEC2_NumLon;
       grid->y.size = ISEC2_NumLat;
 
@@ -395,6 +390,21 @@ void cgribexAddRecord(stream_t * streamptr, int param, int *isec1, int *isec2, d
           double angle = - FSEC2_RotAngle;
           gridDefParamRLL(gridID, xpole, ypole, angle);
         }
+      else if ( gridptr->projtype == CDI_PROJ_LCC )
+        {
+          double lcc_xinc      = ISEC2_Lambert_dx;
+          double lcc_yinc      = ISEC2_Lambert_dy;
+          double lcc_originLon = ISEC2_FirstLon * 0.001;
+          double lcc_originLat = ISEC2_FirstLat * 0.001;
+          double lcc_lonParY   = ISEC2_Lambert_Lov * 0.001;
+          double lcc_lat1      = ISEC2_Lambert_LatS1 * 0.001;
+          double lcc_lat2      = ISEC2_Lambert_LatS2 * 0.001;
+          double lcc_projflag  = ISEC2_Lambert_ProjFlag;
+          double lcc_scanflag  = (short)ISEC2_ScanFlag;
+          gridDefParamLCC(gridID, lcc_originLon, lcc_originLat, lcc_lonParY,
+                          lcc_lat1, lcc_lat2, lcc_xinc, lcc_yinc,
+                          lcc_projflag, lcc_scanflag);
+        }
     }
   else
     Free(gridptr);
@@ -677,10 +687,7 @@ int cgribexScanTimestep1(stream_t *streamptr)
 
       gribDateTime(isec1, &vdate, &vtime);
 
-      if ( ISEC4_NumBits > 0 && ISEC4_NumBits <= 32 )
-	datatype = ISEC4_NumBits;
-      else
-        datatype = CDI_DATATYPE_PACK;
+      datatype = (ISEC4_NumBits > 0 && ISEC4_NumBits <= 32) ? ISEC4_NumBits : CDI_DATATYPE_PACK;
 
       if ( nrecs == 0 )
 	{