diff --git a/definitions/grib1/grid_definition_5.def b/definitions/grib1/grid_definition_5.def index ab8db174f..9ec34a8ee 100644 --- a/definitions/grib1/grid_definition_5.def +++ b/definitions/grib1/grid_definition_5.def @@ -2,7 +2,7 @@ # # This software is licensed under the terms of the Apache Licence Version 2.0 # which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. -# +# # In applying this licence, ECMWF does not waive the privileges and immunities granted to it by # virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. # @@ -57,6 +57,11 @@ alias geography.LaDInDegrees=LaDInDegrees; # Projection centre flag unsigned[1] projectionCentreFlag : dump ; alias projectionCenterFlag=projectionCentreFlag; +# Note our flagbit numbers go from 7 to 0, while WMO convention is from 1 to 8 +# If bit 1 is 0, then the North Pole is on the projection plane +# If bit 1 is 1, then the South Pole is on the projection plane +flagbit southPoleOnProjectionPlane(projectionCentreFlag,7) : dump; # WMO bit 1 + # for change_scanning_direction alias yFirst=latitudeOfFirstGridPointInDegrees; @@ -74,7 +79,8 @@ meta numberOfValues number_of_values(values,bitsPerValue,numberOfDataPoints,bitm iterator polar_stereographic(numberOfPoints,missingValue,values, radius,Nx,Ny, latitudeOfFirstGridPointInDegrees,longitudeOfFirstGridPointInDegrees, - LaDInDegrees,orientationOfTheGridInDegrees, + southPoleOnProjectionPlane, + orientationOfTheGridInDegrees, Dx,Dy, iScansNegatively, jScansPositively, diff --git a/definitions/grib2/template.3.20.def b/definitions/grib2/template.3.20.def index f04b848a7..6d79ee4d9 100644 --- a/definitions/grib2/template.3.20.def +++ b/definitions/grib2/template.3.20.def @@ -52,19 +52,23 @@ alias LoV = orientationOfTheGrid ; meta geography.orientationOfTheGridInDegrees scale(orientationOfTheGrid,oneConstant,grib2divider,truncateDegrees) : dump; # Dx - X-direction grid length -# NOTE 3: Grid length is in units of 10–3 m at the latitude specified by LaD +# NOTE 3: Grid length is in units of 10-3 m at the latitude specified by LaD unsigned[4] Dx : edition_specific; meta geography.DxInMetres scale(Dx,one,thousand,truncateDegrees) : dump; alias xDirectionGridLength=Dx; # Dy - Y-direction grid length -# NOTE 3: Grid length is in units of 10–3 m at the latitude specified by LaD +# NOTE 3: Grid length is in units of 10-3 m at the latitude specified by LaD unsigned[4] Dy : edition_specific; meta geography.DyInMetres scale(Dy,one,thousand,truncateDegrees) : dump; alias yDirectionGridLength=Dy; # Projection centre flag flags[1] projectionCentreFlag 'grib2/tables/[tablesVersion]/3.5.table' : dump; +# Note our flagbit numbers go from 7 to 0, while WMO convention is from 1 to 8 +# If bit 1 is 0, then the North Pole is on the projection plane +# If bit 1 is 1, then the South Pole is on the projection plane +flagbit southPoleOnProjectionPlane(projectionCentreFlag,7) : dump; # WMO bit 1 include "template.3.scanning_mode.def"; @@ -72,7 +76,8 @@ include "template.3.scanning_mode.def"; iterator polar_stereographic(numberOfPoints,missingValue,values, radius,Nx,Ny, latitudeOfFirstGridPointInDegrees,longitudeOfFirstGridPointInDegrees, - LaDInDegrees,orientationOfTheGridInDegrees, + southPoleOnProjectionPlane, + orientationOfTheGridInDegrees, DxInMetres,DyInMetres, iScansNegatively, jScansPositively, diff --git a/src/grib_iterator_class_polar_stereographic.c b/src/grib_iterator_class_polar_stereographic.c index d090acc4a..5ce36eedb 100644 --- a/src/grib_iterator_class_polar_stereographic.c +++ b/src/grib_iterator_class_polar_stereographic.c @@ -105,10 +105,10 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args) double *lats,*lons; double lonFirstInDegrees,latFirstInDegrees,lonFirst,latFirst,radius=0; long nx,ny,standardParallel,centralLongitude; - double phi,lambda0,xFirst,yFirst,x,y,Dx,Dy; + double lambda0,xFirst,yFirst,x,y,Dx,Dy; double k,sinphi1,cosphi1; long alternativeRowScanning,iScansNegatively; - long jScansPositively,jPointsAreConsecutive; + long jScansPositively,jPointsAreConsecutive, southPoleOnPlane; double sinphi,cosphi,cosdlambda,sindlambda; double cosc,sinc; long i,j; @@ -120,7 +120,7 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args) const char* sny = grib_arguments_get_name(h,args,self->carg++); const char* slatFirstInDegrees = grib_arguments_get_name(h,args,self->carg++); const char* slonFirstInDegrees = grib_arguments_get_name(h,args,self->carg++); - const char* sstandardParallel = grib_arguments_get_name(h,args,self->carg++); + const char* ssouthPoleOnPlane = grib_arguments_get_name(h,args,self->carg++); const char* scentralLongitude = grib_arguments_get_name(h,args,self->carg++); const char* sDx = grib_arguments_get_name(h,args,self->carg++); const char* sDy = grib_arguments_get_name(h,args,self->carg++); @@ -129,8 +129,7 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args) const char* sjPointsAreConsecutive = grib_arguments_get_name(h,args,self->carg++); const char* salternativeRowScanning = grib_arguments_get_name(h,args,self->carg++); double c,rho; - /*double pi4=acos(0.0)/2.0;*/ - sinphi1 = cosphi1 = phi = 0.0; + sinphi1 = cosphi1 = 0.0; if((ret = grib_get_double_internal(h, sradius,&radius)) != GRIB_SUCCESS) return ret; @@ -147,7 +146,7 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args) return ret; if((ret = grib_get_double_internal(h, slonFirstInDegrees,&lonFirstInDegrees)) != GRIB_SUCCESS) return ret; - if((ret = grib_get_long_internal(h, sstandardParallel,&standardParallel)) != GRIB_SUCCESS) + if((ret = grib_get_long_internal(h, ssouthPoleOnPlane,&southPoleOnPlane)) != GRIB_SUCCESS) return ret; if((ret = grib_get_long_internal(h, scentralLongitude,¢ralLongitude)) != GRIB_SUCCESS) return ret; @@ -164,6 +163,7 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args) if((ret = grib_get_long_internal(h, salternativeRowScanning,&alternativeRowScanning)) != GRIB_SUCCESS) return ret; + standardParallel = (southPoleOnPlane == 1) ? -90 : +90; sinphi1 = sin(standardParallel*DEG2RAD); cosphi1 = cos(standardParallel*DEG2RAD); lambda0 = centralLongitude*DEG2RAD;