Re: sorting daily precipitation data.

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Mon Apr 12 2010 - 16:26:28 MDT

Try .... a positive

       ip = dim_pqsort_n(x, 2, 0)

You can print out values at some selected grid point.

       print(x(prStrt:,jj,ii))

where jj and ii are user specified index values

---
Good Luck
On 04/12/2010 04:06 PM, Sho Kawazoe wrote:
> Dennis/NCL users,
>
>   I'm running into some sorting issues with this script. I'm working
> with daily accumulated precipitation, and looking for a method to sort
> the data to ultimately analyze the top 10% extreme precipitation events.
> I need each precipitation data to have the associated dimensions
> (time,yc,xc) along with it, to determine when these events occur.
>
> However, output of PrcExtreme data results in all zeros (probably
> indicating that it's reading the bottom 10%). Also, the time output is
> showing the last 10% of the entire time file, so the metadata is not
> being carried along with the associated precipitation variable. I'm
> aware that the dim_pqsort_n function for this script sorts both x and
> time, but I was wondering if anyone had a clever method around this.
>
> begin
>
> ;********************************************************
> ; Read original NetCDF file. Input necessary
> ; files
> ;*********************************************************
>
>    fname = ("region_output.nc <http://region_output.nc>")
>    fin = addfile (fname,"r")
>    prc = fin->pr                     ; (time,yc, xc) 3 hrly
>    lat     = fin->lats_region                 ; lat(yc, xc)
>    lon     = fin->lons_region                 ; lon(yc, xc)
>    time    = fin->time
>    Lambert_Conformal = fin->Lambert_Conformal
>
>    dimx3   = dimsizes(prc)
>    NTIM    = dimx3(0)               ; TOTAL number of time steps
>    print("NTIM="+NTIM)
>
> ;************************************************************
> ; Create daily averages.
> ;************************************************************
>
>      ntJump = 8
>
>      x = prc(2::ntJump,:,:) ; trick: create array with meta
>      printVarSummary(x) ; values will be overwritten with averages
>
>      ntStrt = 3
>      ntLast = ntStrt + ntJump - 1
>
>      do nt=0,NTIM-9,8
>         x(nt/8,:,:) = (dim_avg_n_Wrap(prc(ntStrt:ntLast,:,:), 0))
>         ntStrt = ntStrt + ntJump
>         ntLast = ntLast + ntJump
>      end do
>      printVarSummary(x)
>      printMinMax(x,True)
>
>      printMinMax(time,True)
>
>      x@info_tag = "daily average"
>
> ;**************************************************************
> ; Convert units kg/(m2-s) to mm/day
> ;**************************************************************
>
>      x = x*86400.
>      x@units = "mm/day"
>      printVarSummary(x)
>      printMinMax(x,True)
>
> ;**************************************************************
> ; Sorting subroutine. Add if your looking for extreme
> ; events!
> ;*************************************************************
>
>     ip = dim_pqsort_n(x, -2, 0)  ; 'x' will be in descending order
>     delete(ip)
>
>     dimx = dimsizes(x)
>     ntim = dimx(0)
>
>     pc   = 10                 ; 10%
>     prStrt = floattoint( ((100-pc)/100.)*ntim  )  ; index value
>
>
> ;*********************************************************
> ; Create a netCDF with the top 10% of values.
> ;*********************************************************
>
>     diro = "./"
>     filo = "Precip_upper_"+pc+".nc"
>
>     system("/bin/rm -f "+diro+filo)        ; remove any pre-existing file
>     ncdf = addfile(diro+filo ,"c")         ; open output netCDF file
>
>     PrcExtreme = x(prStrt::,:,:)
>
>     ncdf->PrcExtreme = PrcExtreme     ; top 'pc'% only
>     ncdf->lat = lat
>     ncdf->lon = lon
>     ncdf->Lambert_Conformal = Lambert_Conformal
>
>     printMinMax(PrcExtreme, True)
>
>     globeAtt              = 1
>     globeAtt@title        = "MM5: daily averages"
>     globeAtt@acronym      = "MM5: Mesoscale Model"
>     globeAtt@description  = "Data analyzed uses archived NCEP
> reanalysis. This project is in association with the NARCCAP project"
>     globeAtt@process      = "NetCDF is sorted output (top 10% extreme
> events)"
>
>     fileattdef(ncdf ,globeAtt)
>
> end
>
>
> --
> Sho Kawazoe
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-- 
======================================================
Dennis J. Shea                  tel: 303-497-1361    |
P.O. Box 3000                   fax: 303-497-1333    |
Climate Analysis Section                             |
Climate & Global Dynamics Div.                       |
National Center for Atmospheric Research             |
Boulder, CO  80307                                   |
USA                        email: shea 'at' ucar.edu |
======================================================
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Apr 12 16:26:44 2010

This archive was generated by hypermail 2.1.8 : Wed Apr 14 2010 - 09:15:22 MDT