Re: Is there an easier way to add a "data buffer" around a subset?

Hi Bridget

The attached function has not been well tested.
You may need to alter for your needs.

Happy New Year
D

Bridget Thrasher wrote:
> I have some data files that contain values only over the continental US
> and fill values in the other cells. In order to not lose data when
> interpolating to a finer grid, I decided to add some buffer data around
> the group of cells that contain values. Right now I'm doing it with
> several loops and if statements (see below if interested) that basically
> copy values first along the west coast, then the east coast, then the
> southern cells, and then the northern cells. This is the only solution I
> could come up with, and while it works, it's not particularly speedy. Is
> there a better way?
>
> -Bridget
>
>
> do i = 0,nlats-1
> do j = 1,nlons-2
> if
> (num(ismissing(aggobs(:,i,j))).gt.0.and.num(ismissing(aggobs(:,i,j-1))).gt.0.and.num(ismissing(aggobs(:,i,j+1))).eq.0)
> then
> aggobs(:,i,j) = aggobs(:,i,j+1)
> end if
> end do
> do j = nlons-2,1,1
> if
> (num(ismissing(aggobs(:,i,j))).gt.0.and.num(ismissing(aggobs(:,i,j+1))).gt.0.and.num(ismissing(aggobs(:,i,j-1))).eq.0)
> then
> aggobs(:,i,j) = aggobs(:,i,j-1)
> end if
> end do
> if
> (num(ismissing(aggobs(:,i,0))).gt.0.and.(num(ismissing(aggobs(:,i,nlons-1))).eq.0.or.num(ismissing(aggobs(:,i,1))).eq.0))
> then
> aggobs(:,i,0) = avg((/aggobs(:,i,nlons-1),aggobs(:,i,1)/))
> end if
> if
> (num(ismissing(aggobs(:,i,nlons-1))).gt.0.and.(num(ismissing(aggobs(:,i,nlons-2))).eq.0.or.num(ismissing(aggobs(:,i,0))).eq.0))
> then
> aggobs(:,i,nlons-1) = avg((/aggobs(:,i,nlons-2),aggobs(:,i,0)/))
> end if
> end do
> do j = 0,nlons-1
> do i = 1,nlats-2
> if
> (num(ismissing(aggobs(:,i,j))).gt.0.and.num(ismissing(aggobs(:,i-1,j))).gt.0.and.num(ismissing(aggobs(:,i+1,j))).eq.0)
> then
> aggobs(:,i,j) = aggobs(:,i+1,j)
> end if
> end do
> do i = nlats-2,1,1
> if
> (num(ismissing(aggobs(:,i,j))).gt.0.and.num(ismissing(aggobs(:,i+1,j))).gt.0.and.num(ismissing(aggobs(:,i-1,j))).eq.0)
> then
> aggobs(:,i,j) = aggobs(:,i-1,j)
> end if
> end do
> if
> (num(ismissing(aggobs(:,0,j))).gt.0.and.(num(ismissing(aggobs(:,nlats-1,j))).eq.0.or.num(ismissing(aggobs(:,1,j))).eq.0))
> then
> aggobs(:,0,j) = avg((/aggobs(:,nlats-1,j),aggobs(:,1,j)/))
> end if
> if
> (num(ismissing(aggobs(:,nlats-1,j))).gt.0.and.(num(ismissing(aggobs(:,nlats-2,j))).eq.0.or.num(ismissing(aggobs(:,0,j))).eq.0))
> then
> aggobs(:,nlats-1,j) = avg((/aggobs(:,nlats-2,j),aggobs(:,0,j)/))
> end if
> end do
>
>
> --
> Bridget Thrasher, PhD
> Postdoctoral Researcher
> Climate Central
> www.climatecentral.org <http://www.climatecentral.org>
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

function addBuffer3 ( x[*][*][*], opt[1]:integer )
local dimx, nt, ny, mx, ny2, mx2, xBuf
begin
  dimx = dimsizes(x)
  nt = dimx(0)
  ny = dimx(1)
  mx = dimx(2)

  ny2 = ny+2
  mx2 = mx+2

  xBuf = new ( (/nt,ny2,mx2/), typeof(x), getFillValue(x))

  xBuf(:,1:ny2-2,1:mx2-2) = (/ x /) ; main body

  if (opt.eq.0) then ; assign boundary vals
      xBuf(:,1:ny2-2,0) = (/ x(:,:,0) /) ; left
      xBuf(:,1:ny2-2,mx2-1) = (/ x(:,:,mx-1) /) ; right
      xBuf(:,0,1:mx2-2) = (/ x(:,0,:) /) ; bottom
      xBuf(:,ny2-1,1:mx2-2) = (/ x(:,ny-1,:) /) ; top
                                                     ; corners
      xBuf(:,0 ,0 ) = x(:,0,0) ; bottom left
      xBuf(:,0 ,mx2-1) = x(:,0,mx-1) ; bottom right
      xBuf(:,ny2-1,0 ) = x(:,ny-1,0) ; upper left
      xBuf(:,ny2-1,mx2-1) = x(:,ny-1,mx-1) ; upper right
  end if

  if (opt.eq.1) then ; linear extrapolation
      xBuf(:,1:ny2-2,0) = (/ 2*x(:,:,0) -x(:,:,1) /) ; left
      xBuf(:,1:ny2-2,mx2-1) = (/ 2*x(:,:,mx-1)-x(:,:,mx-2) /) ; right
      xBuf(:,0,1:mx2-2) = (/ 2*x(:,0,:) -x(:,1,:) /) ; bottom
      xBuf(:,ny2-1,1:mx2-2) = (/ 2*x(:,ny-1,:)-x(:,ny-2,:) /) ; top
                                                     ; corners
      xBuf(:,0 ,0 ) = 0.5*(xBuf(:,0,1)+xBuf(:,1,0)) ; bottom left
      xBuf(:,0 ,mx2-1) = 0.5*(xBuf(:,0 ,mx2-2)+xBuf(:,1 ,mx2-1)) ; bottom right
      xBuf(:,ny2-1,0 ) = 0.5*(xBuf(:,ny2-2,0 )+xBuf(:,ny2-1,1 )) ; upper left
      xBuf(:,ny2-1,mx2-1) = 0.5*(xBuf(:,ny2-1,mx2-2)+xBuf(:,ny2-2,mx2-1)) ; upper right
  end if

  copy_VarAtts(x, xBuf) ; contributed.ncl
  return( xBuf )
end
   dir = "/Users/shea/Data/netCDF/"
   f = addfile (dir+"uv300.nc", "r")
   u = f->U ; (ntim,ny,mx)
   printMinMax(u, True)

   U0= addBuffer3( u, 0 ) ; (ntim,NY,MX) => (ntim,ny+2,mx+2)
   printMinMax(U0, True)

   U1= addBuffer3( u, 1 ) ; (ntim,NY,MX) => (ntim,ny+2,mx+2)
   printMinMax(U1, True)

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Jan 8 15:27:54 2010