Gridded data area average

NCL users,

 I'm currently working with University of Washington gridded precipitation
data, which is at 1/8 by 1/8 degree grid. The coordinates are in 1-D
latitude/longitude, with the precipitation dimensions being
(time,latitude,longitude).

Unfortunately, the data is too large for a separate Fortran Script to
handle, so I need to convert it to a 1/2 by 1/2 degree grid. My current
method involves simply taking the data at every 1/2 degree, but some
questions arose about how effective and accurate this really is. Instead, I
want to still have a 1/2 by 1/2 degree data reading, but instead of a simple
point analysis, I want to average values of the data from the points around
it (so technically it would be the average of 16points around one point).

 I've attached a script, and as it shows, I've tried a variety of methods,
but all with various issues (such as dimensions being incorrect, all values
are zeros, do loop not actually looping anything etc). Am I on the right
track? Do I regrid? Any tip or suggestions is greatly appreciated.

Thank You
Sho Kawazoe

;*******************************************************************
;Load NCL library files. Always load "Shea_util" AFTER contributed
;*******************************************************************

  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/contributed.ncl"
  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin

;*****************************************************
; Read original NetCDF file.
;*****************************************************

  fin = addfile("UW_region.nc", "r")

  prc = fin->prc(:,:,:)
  lat = fin->latitude(:)
  lon = fin->longitude(:)
  time = fin->time(:)

  dim = dimsizes(prc)
  ntime = dim(0)
  nlat = dim(1)
  nlon = dim(2)

  x = prc(:,0:14,0:14)

  printMinMax(x,True)

  y0 = 44.3125
  x0 = -93.4375

; x2 = linint2_points(lon,lat,prc,False,x0,y0,0)
; printVarSummary(x2) ; Another Method (only creates a 2D array)
; printMinMax(x2,True)

; do k = 0,1
; do j = 0,5
; do i = 0,5
; x(k,j,i) = prc(k,j,i) + prc(k,j+1,i) + prc(k,j,i+1) + prc(k,j+1,j+1)
doloop method 1
; prc_ave = prc(k,j,i) + prc(k,j+1,i) + prc(k,j,i+1) + prc(k,j+1,j+1)
doloop method 2
; end do
; end do
; end do

  ntStrt = 0
  ntLast = 1

  do nt = 0,5,1
    x(nt,:,:) = dim_avg_n_Wrap(prc(:,ntStrt:ntLast,ntStrt:ntLast),0)
    ntStrt = ntStrt + 1 ; Another Method (Dim on right is not the
same as the left)
    ntLast = ntLast + 1
  end do

  print("")
  printMinMax(x,True)
  printVarSummary(x)
; print(x)

; prc_ave!0 = "time"
; prc_ave!1 = "latitude"
; prc_ave!2 = "longitude"

;******************************************************************
; Create a new nc output file. create global attributes
;******************************************************************

   system("/bin/rm -f UW_domAve.nc") ; remove any pre-existing file
   fout = addfile("UW_domAve.nc","c")
   filedimdef(fout,"time",-1,True) ; Makes time an UNLIMITED
dimension (as suggested by NCL)

   globeAtt = 1
   globeAtt@description = "Data analyzed uses gridded UW observational
data. This project is in association with the NARCCAP project"
   globeAtt@title = "University of Washington gridded observational
data output prepared for NARCCAP present-day climate"
   globeAtt@project_id = "NARCCAP"
   globeAtt@contact = "Edwin Mauer: emaurer@engr.scu.edu"
   globeAtt@instituion = "University of Washington Department of
Atmospheric Sciences (Seattle, Washington)"

   fileattdef(fout, globeAtt) ; create the global [file] attributes

; fout->prc = x2
   fout->time = time

end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Nov 15 15:25:59 2010