GrADS源代码pdef插值部分

rookie

以下是GrADS中关于Grd数据中有pdef变量的插值，红色部分为pdef说明文档中的双线性插值（Bi-linear interpolation ，参考How         PDEF Grid Interpolation Works，http://grads.iges.org/grads/gadoc/pdef.html），以文档为例
经纬度转化为网格点坐标31.24 ，67.88 后（可使用gaddes.c文件中 ll2lc函数），原始网格中四个点行列号为：
点1：31,67  
点2 ：31,68
点3：32,67
点4：32,68
对应红色代码中P[0]为点1的变量值，P[1]为点2的变量值..........
dx=31.24-31； dy=67.88-67；



for (i=0; i<len; i++) {
      ig = (y-1) * pfi->dnum[0] + x + i - 1;
      ioff = *(pfi->ppi[0]+ig);              /* ioff index values start at 0 */
      if (ioff<0) {
             *gr = pgr->undef;
             *gru = 0;
      }
      else {
         dx = (gadouble)*(pfi->ppf[0]+ig);
         dy = (gadouble)*(pfi->ppf[1]+ig);
         pos = (e-1)*(pfi->dnum[3]*pfi->tsiz) + (tt-1)*(pfi->tsiz) + pvr->offset + (z-1)*(pfi->gsiz) + ioff;

            /* Get the first two pre-projected grid values */
        if (pfi->idxflg) {                     
               rc = gairow(x,y,z,t,e,ioff,2,p,umask);                   /* grib */
            }else if (pfi->ncflg==1) {               
               rc = gancsetup();
               if (rc) return(rc);
               ncig = (gaint)(1 + ioff%pfi->ppisiz);
            ncjg = (gaint)(1 + ioff/pfi->ppisiz);
                rc = gancrow(ncig,ncjg,z,tt,e,2,p,umask);                /* netcdf */  
            }else if (pfi->ncflg==2) {              
                ncig = (gaint)(1 + ioff%pfi->ppisiz);
            ncjg = (gaint)(1 + ioff/pfi->ppisiz);
                rc = gahrow(ncig,ncjg,z,tt,e,2,p,umask);                 /* hdf */  
            }else if (pfi->ncflg==3) {              
                ncig = (gaint)(1 + ioff%pfi->ppisiz);
            ncjg = (gaint)(1 + ioff/pfi->ppisiz);
                rc = gah5row(ncig,ncjg,z,tt,e,2,p,umask);                 /* hdf5 */  
            }else {                              
                rc = garead(pos,2,p,umask);                              /* binary */
            }
        if (rc) return(rc);

            /* Get the second two pre-projected grid values */
        if (pfi->idxflg) {
               rc = gairow(x,y,z,t,e,ioff+pfi->ppisiz,2,p+2,umask+2);   /* grib */
            }else if (pfi->ncflg==1) {
           ncjg++;
               rc = gancrow(ncig,ncjg,z,tt,e,2,p+2,umask+2);            /* netcdf */
            }else if (pfi->ncflg==2) {
           ncjg++;
               rc = gahrow(ncig,ncjg,z,tt,e,2,p+2,umask+2);             /* hdf */  
            }else if (pfi->ncflg==3) {
           ncjg++;
               rc = gah5row(ncig,ncjg,z,tt,e,2,p+2,umask+2);             /* hdf5 */  
            }else {
               rc = garead(pos+pfi->ppisiz,2,p+2,umask+2);              /* binary */  
            }
        if (rc) return(rc);

            /* Do the bilinear interpolation, as long as we have no undefs */
        if (umask[0]==0 || umask[1]==0 || umask[2]==0 || umask[3]==0) {
               *gr = pgr->undef;
               *gru = 0;
            }else {
          g1 = p[0] + (p[1]-p[0])*dx;
           g2 = p[2] + (p[3]-p[2])*dx;
           *gr = g1 + (g2-g1)*dy;
               *gru = 1;
        }
      }
      gr++; gru++;
    }
插值代码在gaio.c文件中，以上是我阅读GrADS源代码理解的，针对pdef 为lcc，验证后和GrADS读数有些差别，望大家看下，多交流


void ll2lc (gadouble *vals, gadouble grdlat, gadouble grdlon, gadouble *grdi, gadouble *grdj, gadouble *wrot) {

/*  Subroutine to convert from lat-lon to lambert conformal i,j.
    Provided by NRL Monterey; converted to C 6/15/94.

c          SUBROUTINE: ll2lc
c
c          PURPOSE: To compute i- and j-coordinates of a specified
c                   grid given the latitude and longitude points.
c                   All latitudes in this routine start
c                   with -90.0 at the south pole and increase
c                   northward to +90.0 at the north pole.  The
c                   longitudes start with 0.0 at the Greenwich
c                   meridian and increase to the east, so that
c                   90.0 refers to 90.0E, 180.0 is the inter-
c                   national dateline and 270.0 is 90.0W.
c
c          INPUT VARIABLES:
c
c  vals+0    latref: latitude at reference point (iref,jref)
c                     
c  vals+1    lonref: longitude at reference point (iref,jref)
c                     
c  vals+2    iref:   i-coordinate value of reference point
c                     
c  vals+3    jref:   j-coordinate value of reference point
c                     
c  vals+4    stdlt1: standard latitude of grid (S True lat)
c                     
c  vals+5    stdlt2: second standard latitude of grid (only required
c                    if igrid = 2, lambert conformal) (N True lat)
c                     
c  vals+6    stdlon: standard longitude of grid (longitude that
c                     points to the north)
c                     
c  vals+7    delx:   grid spacing of grid in x-direction
c                    for igrid = 1,2,3 or 4, delx must be in meters
c                    for igrid = 5, delx must be in degrees
c                     
c  vals+8    dely:   grid spacing (in meters) of grid in y-direction
c                    for igrid = 1,2,3 or 4, delx must be in meters
c                    for igrid = 5, dely must be in degrees
c
c            grdlat: latitude of point (grdi,grdj)
c
c            grdlon: longitude of point (grdi,grdj)
c
c            grdi:   i-coordinate(s) that this routine will generate
c                    information for
c
c            grdj:   j-coordinate(s) that this routine will generate
c                    information for
c
*/

  gadouble pi, pi2, pi4, d2r, r2d, radius, omega4;
  gadouble gcon,ogcon,H,deg,cn1,cn2,cn3,cn4,rih,xih,yih,rrih,check;
  gadouble alnfix,alon,x,y,windrot;
  gadouble latref,lonref,iref,jref,stdlt1,stdlt2,stdlon,delx,dely;

  pi = M_PI;
  pi2 = pi/2.0;
  pi4 = pi/4.0;
  d2r = pi/180.0;
  r2d = 180.0/pi;
  radius = 6371229.0;
  omega4 = 4.0*pi/86400.0;

  latref = *(vals+0);   
  lonref = *(vals+1);   
  iref   = *(vals+2);   
  jref   = *(vals+3);   
  stdlt1 = *(vals+4);   
  stdlt2 = *(vals+5);   
  stdlon = *(vals+6);   
  delx   = *(vals+7);   
  dely   = *(vals+8);   

/* case where standard lats are the same */
/* corrected by Dan Geiszler of NRL; fabs of the
   lats was required for shem cases */

  if (stdlt1 == stdlt2) {
    gcon = sin(d2r*(fabs(stdlt1)));
  } else {
    gcon = (log(sin((90.0-fabs(stdlt1))*d2r))
           -log(sin((90.0-fabs(stdlt2))*d2r)))
          /(log(tan((90.0-fabs(stdlt1))*0.5*d2r))
           -log(tan((90.0-fabs(stdlt2))*0.5*d2r)));
  }
  ogcon = 1.0/gcon;
  H = fabs(stdlt1)/(stdlt1);        /* 1 for NHem, -1 for SHem */
  cn1 = sin((90.0-fabs(stdlt1))*d2r);
  cn2 = radius*cn1*ogcon;
  deg = (90.0-fabs(stdlt1))*d2r*0.5;
  cn3 = tan(deg);
  deg = (90.0-fabs(latref))*d2r*0.5;
  cn4 = tan(deg);
  rih = cn2*pow((cn4/cn3),gcon);

  xih =  rih*sin((lonref-stdlon)*d2r*gcon);
  yih = -rih*cos((lonref-stdlon)*d2r*gcon)*H;
  deg = (90.0-grdlat*H)*0.5*d2r;
  cn4 = tan(deg);
  rrih = cn2*pow((cn4/cn3),gcon);
  check  = 180.0-stdlon;
  alnfix = stdlon+check;
  alon   = grdlon+check;

  while (alon<  0.0) alon = alon+360.0;
  while (alon>360.0) alon = alon-360.0;

  deg = (alon-alnfix)*gcon*d2r;
  x =  rrih*sin(deg);
  y = -rrih*cos(deg)*H;
  *grdi = iref + (x-xih)/delx;
  *grdj = jref + (y-yih)/dely;
  /* mf 20040630 -- use ftp://grads.iges.org/grads/src/grib212.f to calc rotation factor */
  windrot=gcon*(stdlon-grdlon)*d2r;
  *wrot=windrot;
}

rookie

针对计算经纬度对应的原始网格行列号，使用ArcMap添加一个站点，和使用GrADS内部的ll2lc()函数计算出的差别还是很大，计算出来到48.8,54.9，ArcMap中大概在48.5,54.5左右。

dannybear

rookie 发表于 2015-9-16 18:43
针对计算经纬度对应的原始网格行列号，使用ArcMap添加一个站点，和使用GrADS内部的ll2lc()函数计算出的差别 ...

请问楼主，pdef，是不是表示这个数据的排列是按照兰伯特投影的那种等角排列啊，就是不是规则的网格点是吗？即 x方向和y方向的间距不一样

rookie

dannybear 发表于 2016-3-25 18:24
请问楼主，pdef，是不是表示这个数据的排列是按照兰伯特投影的那种等角排列啊，就是不是规则的网格点是吗 ...

这要看pdef参数，如果后面两个参数一样，那就是正方形网格

dannybear

rookie 发表于 2016-3-28 11:06
这要看pdef参数，如果后面两个参数一样，那就是正方形网格

请问楼主有没有用grads处理过这中二进制文件

rookie

dannybear 发表于 2016-3-28 14:31
请问楼主有没有用grads处理过这中二进制文件

看过这块的代码，数据没有试过

dannybear

好的，谢谢