[ncl-talk] Extracting "dates" if a threshold is exceeded in any gridpoint of a domain

[Karin Meier-Fleischer]

Hi Lyndon,

first you have to select the sub-region of your data and then get the 
time indices where your data is greater equal than the threshold.

Your script should be looking like:

begin

   f = addfile("vort_700hPa.nc","r")              ;-- open file

   vr = f->vr                                     ;-- read variable

   time = f->time                                 ;-- read time
   ntimes = dimsizes(time)                        ;-- get size of time

;-- define sub-region
   lat0 =   7.5
   lat1 =  27.5
   lon0 =  11.0
   lon1 = 135.0

   vr_region = vr(:,{lat0:lat1},{lon0:lon1})      ;-- select the data of 
the sub-region
   printVarSummary(vr_region)

   threshold = 3e-5                               ;-- define threshold

;-- loop over time to retrieve times where data is greater equal threshold
   j = 0
   do i=0,ntimes-1
      tmp := ndtooned(vr_region(i,:,:))          ;-- convert to 1d array
      if(any(tmp .ge. threshold)) then
         if(j .eq. 0) then
            tt = time(i)                         ;-- first time step 
where threshold is exceeded
         else
            tt := array_append_record (tt, time(i), 0) ;-- append next 
time step where threshold is exceeded
         end if
         j = j+1
      end if
   end do

;-- convert time to readable date string
   utc_date   = cd_calendar(tt, 0)                          ;-- convert 
date to UT-referenced date
   year       = sprinti("%0.4i",tointeger(utc_date(:,0)))   ;-- get year 
as integer value
   mon        = sprinti("%0.2i",tointeger(utc_date(:,1)))   ;-- get 
month as integer value
   day        = sprinti("%0.2i",tointeger(utc_date(:,2)))   ;-- get day 
as integer value
   str_date   = year+"/"+mon+"/"+day                        ;-- yyyy/mm/dd
   print(str_date)

end

Hope this helps,
Karin


Am 26.04.17 um 07:09 schrieb Lyndon Mark Olaguera:
> Dear fellow NCL users,
>
>
> I would like to extract the "date" when any grid point in a domain 
> exceeds a threshold value using NCL (daily time steps). I tried to 
> create a netcdf file of daily vorticity using the script below.
>
> I would like to add a line in this script that will extract the date 
> (day month year) whenever a grid point from this domain is exceeded:
>               LAT : 7.5N to 27.5N
>               LON : 110E to 135E
>
> the threshold value that I am using is 3e-5 per second.
>
>
> I am attaching the netcdf file and the image in this email.​
> test_2003_30.png 
> <https://drive.google.com/file/d/0B9faET7Bc2o8RWRXcTdKZEx0YlE/view?usp=drive_web>
> ​​​
>
> Is this possible in NCL?
> Is there any function in NCL that can extract values from every 
> gridpoint in a specified domain?
>
>
> I'll appreciate if anyone can point me to the right 
> direction/functions to implement this in NCL.
>
> Many thanks in advance!
>
>
> begin
> ;************************************************
> ; variable and file handling
> ;************************************************
>   ufile  = addfile("uwind_may_2003_700hPa.nc","r") ; open netcdf file
>   u   = ufile->uwnd(:,0,:,:)   ; pull u off file
>   lat = ufile->lat
>   lon = ufile->lon
>   vfile  = addfile("vwind_may_2003_700hPa.nc","r")
>   v   = vfile->vwnd(:,0,:,:)   ; pull v off file
> ;************************************************
> ; calculate vorticity on a Gaussian Grid
> ;************************************************
>   ;vrt   = u                                      ; retain coordinates
>   ;vrt   = uv2vrG_Wrap(u,v)
>   vr = uv2vr_cfd (u,v,lat,lon, 3)
>   copy_VarMeta(u,vr)
>   vr at long_name = "vorticity"
>   vr at units     = "per second"
> ;***********************************************
> ; save to a netcdf file
> ;***********************************************
>   ncdf = addfile("vort_700hPa.nc","c")
>   fAtt = True
>   fAtt at title         = "Vorticity at 700hPa"
>   fAtt at source_file   = "daily winds from NCEP-NCAR Reanaysis"
>   fAtt at Conventions   = "None"
>   fAtt at creation_date = systemfunc("date")
>   fileattdef(ncdf,fAtt)            ; copy file attributes
>   filedimdef(ncdf,"time",-1,True)
>   ncdf->vr = vr
> end
>
> *_Lyndz_*
>
>
> _______________________________________________
> 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/20170426/ebbc4bf3/attachment.html