core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

GRIB-574 Polar Stereographic. Use -/+90 for std parallel not LaD

This commit is contained in:
Shahram Najm 2014-10-03 16:45:17 +01:00
parent 45daa319ef
commit dad43b45f5
3 changed files with 22 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
definitions/grib1/grid_definition_5.def
View File

@ -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,

11
definitions/grib2/template.3.20.def
View File

@ -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,

12
src/grib_iterator_class_polar_stereographic.c
View File

@ -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,&centralLongitude)) != 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;