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

[Mary Haley]

Steven,

Thanks for providing the script and data.

Did you actually look at the data in the guatamala_mon_rain.txt file?  The
latitude values go from -92.0772 to -88.5917.  This is not only out of
range of your original latitude data, which goes from -40 to 40, but you
shouldn't have latitudes that go outside of -90 to 90.

Second, your longitudes on this file go from 9 to 366, but the input
longitudes go from -20 to 55.  You need to put the longitudes in the same
range.

Third, and this is just an improvement, but you can replace this do loop:

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

with:

ii = ind(f(:,1).eq.year.and.f(:,2).eq.month)
lat_out = f(ii,4)
lon_out = f(ii,5)

Now you can fix the longitudes to go from -180 to 180 using the "where"
function:

lon_out = where(lon_out.gt.180,lon_out-360,lon_out)

You can now call linint2_points with:

interpolation = linint2_points(lon,lat,rain,False,lon_out,lat_out,0)

However, you still have the problem with the range of "lat_out" being
disjoint from the range of "lat". This is why you are getting all missing
values returned.

--Mary



On Tue, Jan 17, 2017 at 1:36 PM, Steven Fuhrman - NOAA Affiliate <
steven.fuhrman at noaa.gov> wrote:

> 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,Fa
> lse,subset(:,4),subset(:,5),0)   ; subset was formally x
>
> print (interpolation)
> end
>
> --
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170117/9fddb9a2/attachment.html