Stereographic Projection

[Michael Fiorino]

dear reto,

i'm a bit confused about 'ecmwf-products' generated with grads... are
these products plots?

the key thing to remember is there are two (really three) projections in
grads.  for plots or what is displayed on the screen, it is the 'display
projection' which is controlled by 'set mproj'.  the display project can
be be (very) different than the "data projection", e.g., data on a
uniform lat/lon grid but displayed on an orthographic (satellite)
projection...

if your question is about the display projection, then you'd have to
consult the source code for precise details, specifically in gxwmap.c

the relevant code from gxwmap.c for north polar stereographic is below,
but also consult:

http://www.iges.org/grads/gadoc/map.html

best (re)grads, mike


-------------------- gxwmap.c ------------

/* Routine for north polar stereographic.  Projection scaling
   is set along with level 1 linear scaling.   The only difficult
   aspect to this is to set the level 1 linear scaling such that
   the proper aspect ratio is maintained.   */

static float londif;

int gxnste (struct mapprj *mpj) {
gadouble x1,x2,y1,y2,dum,lonave;
gadouble w1,xave,yave;
gadouble lonmn, lonmx, latmn, latmx, xmin, xmax, ymin, ymax;

  lonmn = mpj->lnmn; lonmx = mpj->lnmx;
  latmn = mpj->ltmn; latmx = mpj->ltmx;
  xmin = mpj->xmn; xmax = mpj->xmx;
  ymin = mpj->ymn; ymax = mpj->ymx;

  if ((lonmx-lonmn) > 360.0) {
    return (1);
  }
  if (lonmn<-360.0 || lonmx>360.0) {
    return (1);
  }
  if (latmn<-80.0 || latmx>90.0) {
    return (1);
  }
  if (latmn>=latmx||lonmn>=lonmx||xmin>=xmax||ymin>=ymax) {
    return(1);
  }

  lonave = (lonmx+lonmn)/2.0;            /* Longitude adjustment to put */
  londif = -90.0 - lonave;               /*  central meridian at bottom.*/
  lonref = lonave;

  /* Plotting limits depend on how much of the hemisphere we are
     actually plotting.  */

  if ( (lonmx-lonmn) < 180.0 ) {

     gxnpst ( lonmn, latmn, &x1, &dum );         /* Left side coord  */
     gxnpst ( lonmx, latmn, &x2, &dum );         /* Right side coord */
     gxnpst ( lonmn, latmx, &dum, &y2 );         /* Top coord        */
     gxnpst ( lonave, latmn, &dum, &y1 );        /* Bottom coord     */

  } else {

     gxnpst ( lonave-90.0, latmn, &x1, &dum );   /* Left side coord  */
     gxnpst ( lonave+90.0, latmn, &x2, &dum );   /* Right side coord */
     gxnpst ( lonmn, latmn, &dum, &y2 );         /* Top coord        */
     gxnpst ( lonave, latmn, &dum, &y1 );        /* Bottom coord     */

  }

  /* Set up linear level scaling while maintaining aspect ratio.   */

  if ( ((xmax-xmin)/(ymax-ymin)) > ((x2-x1)/(y2-y1)) ) {
    w1 = 0.5*(ymax-ymin)*(x2-x1)/(y2-y1);
    xave = (xmax+xmin)/2.0;
    gxscal ( xave-w1, xave+w1, ymin, ymax, x1, x2, y1, y2 );
    mpj->axmn = xave-w1;  mpj->axmx = xave+w1;
    mpj->aymn = ymin;  mpj->aymx = ymax;
  } else {
    w1 = 0.5*(xmax-xmin)*(y2-y1)/(x2-x1);
    yave = (ymax+ymin)/2.0;
    gxscal ( xmin, xmax, yave-w1, yave+w1, x1, x2, y1, y2 );
    mpj->axmn = xmin;  mpj->axmx = xmax;
    mpj->aymn = yave-w1;  mpj->aymx = yave+w1;
  }

  gxproj (gxnpst);
  gxback (gxnrev);
  adjtyp = 1;
  return (0);

}

void gxnpst (gadouble rlon, gadouble rlat, gadouble *x, gadouble *y) {
gadouble radius,theta;

  radius = tan (0.785315-(0.00872572*rlat));
  theta = (rlon+londif)*0.0174514;
  *x = radius * cos(theta);
  *y = radius * sin(theta);
}

/* Routine for back transform for npst */

void gxnrev (gadouble x, gadouble y, gadouble *rlon, gadouble *rlat) {
gadouble rad,alpha;

  rad = hypot(x,y);
  alpha = 180.0*atan(rad)/pi;
  *rlat = 90.0 - 2.0*alpha;

  if (x==0.0 && y==0.0) *rlon = 0.0;
  else {
    *rlon = (180.0*atan2(y,x)/pi)-londif;
    while (*rlon < lonref-180.0) *rlon += 360.0;
    while (*rlon > lonref+180.0) *rlon -= 360.0;
  }
}


Reto Stauffer (UNI) wrote:
> Dear GrADS Usergroup
>
> I am working on a MSG satellite visualisation project at our University. We
> already have a lot of ECMWF-Products which are generatet by a GrADS script.
>
> Now we would like to reproduce the same projection with Matlab (Matlab is
> our major script language). The problem is: there is just one line wich
> defines the current projection with mproj and lat/lon bordervalues.
>
> Now my problem: if i am right GrADS use a standard projection point. I
> searched a lot in the internet but cannot find any proper solvements. There
> is just this (http://www.iges.org/grads/gadoc/gradcomdsetmproj.html) =).
>
> Does anyone has some hints and tips about the (standard) used radius and
> projection-point in GrADS?
>
> Thanks a lot
> Best wishes
> Reto Stauffer
> Institute for Meteorology and Geophysics
> UNI Innsbruck Austria