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

[Alan Brammer]

Slight alternative to Karen's suggestion.  Same basic idea but without the loop in the script it self.


load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
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

  vr_max = dim_max_n(vr_region, (/1,2/))     ;; max over lat/lon per time. 
  time_inds = ind(vr_max .ge. threshold)       ;;  if the max isn’t greater than threshold then no values are. 

  if(all(ismissing(time_inds)))
    print("No Times exceed threshold")
    status_exit(1)
  end if

  times = time(time_inds)
  print(cd_string(times, ""))
end

  


> On 26 Apr 2017, at 12:20, Karin Meier-Fleischer <meier-fleischer at dkrz.de> wrote:
> 
> 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 <mailto:ncl-talk at ucar.edu>
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk <http://mailman.ucar.edu/mailman/listinfo/ncl-talk>
> _______________________________________________
> 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/aa715c73/attachment.html