[ncl-talk] Interpolating from a grid to a set of points

[Steven Fuhrman - NOAA Affiliate]

Hello All,

I have gridded satellite estimated rainfall and a file of precipitation
measurements from gauges. I would like to Interpolate the gridded data to
the gauge locations for comparison. I figured the function linint2_points
might work for this process, but thus far I've had no luck. I get no error
message but the return value is just an array of missing values  (-999). My
best guess is that my "rain" array isn't set up correctly for the function
when it is read in from binary.

The 2 files have been placed on the FTP

I am running ncl version 6.0.0

thanks,
Steven Fuhrman

CODE:

begin

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

;***********************************************************
; read in station data file
;***********************************************************
  filename = "guatamala_mon_rain.txt"
  rows = numAsciiRow(filename)
  colls = numAsciiCol(filename)
  f = asciiread(filename,(/rows,colls/),"float")

;***********************************************************
; read in gridded estimated precip file
;***********************************************************

  diri  = "/cpc/fews/production/ARC/MONTHLY/"   ; input directory
  fili  = "arc2_monthly_200010"     ; DSET
  fName = diri+fili                       ; path

  nlat  = 801                             ; YDEF
  nlon  = 751                            ; XDEF
  nvars = 1                               ; vars

  year  = 2000                            ; TDEF
  ntim  = 1                             ; time steps
  nmos  = 1

                                          ;create an array to contain data
  UNDEF = -999.0                          ; UNDEF
  rain     = new ( (/ntim,nlat,nlon/), float, UNDEF)   ;switched!

  rain at long_name = "precipitation"              ; VARS
  rain at units     = "mm"


  do nt=0,ntim-1                          ; read each record: store in x
     rain(nt,:,:) = fbindirread(fName, nt, (/nlat,nlon/), "float")
  end do

  rain!0 = "time"
  rain!1 = "lat"
  rain!2 = "lon"
  printVarSummary(rain)
  print ("min(rain)="+min(rain))
  print ("max(rain)="+max(rain))

;************************************************************
; Generate lat/lon that pairs with gridded data
;************************************************************
  lon = fspan(-20,55,751)
  lon!0 = "lon"
  lon at long_name = "longitude"
  lon at units     = "degrees_east"

  lat = fspan(-40,40,801)
  lat!0 = "lat"
  lat at long_name = "latitude"
  lat at units     = "degrees_north"

;************************************************************
; count number of stations
;************************************************************
count = 0
previous = f(0,0)
do i=1,rows-2
  if(f(i,0) .ne. previous)
    count = count + 1
;    print(f(i,0)+" i am counting!"+i)
  end if
  previous = f(i,0)
end do
print("count = "+count)

;************************************************************
; make subset -an array with only a single date (year/month) for each
station.
;************************************************************
subset =new((/count,colls/),float,-999)
year = 2000
month = 1
line = 0
  do i=0,rows-1
    if ((f(i,1) .eq. year) .and. (f(i,2) .eq. month))
      subset(line,0) = f(i,0)
      subset(line,1) = f(i,1)
      subset(line,2) = f(i,2)
      subset(line,3) = f(i,3)
      subset(line,4) = f(i,4)
      subset(line,5) = f(i,5)
      line = line + 1
    end if
  end do

;***********************************************************
; perform interpolation with subset
;***********************************************************
interpolation = linint2_points(lon,lat,rain,False,subset(:,4),subset(:,5),0)
; subset was formally x

print (interpolation)
end

--
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170117/917e8af1/attachment.html