[ncl-talk] Extracting "dates" if a threshold is exceeded in any gridpoint of a domain
Lyndon Mark Olaguera
olagueralyndonmark429 at gmail.com
Wed Apr 26 20:19:39 MDT 2017
Dear all,
Many many thanks for the help!
Sincerely,
*Lyndz*
On Thu, Apr 27, 2017 at 2:54 AM, Dennis Shea <shea at ucar.edu> wrote:
> Well, I guess 'great minds think alike' ....
> Alan ... What r u doing on my turf????? :-)
>
> I had basically used Alan's approach. I have some local files which had
> daily data for one year.
> The threshold=3e-5 yielded many days exceeding this value. For the
> example. I used 5e-5
>
> =======
> SCRIPT
> =======
> ************************************************
> ; variable and file handling
> ;************************************************
> ;ufile = addfile("uwind_may_2003_700hPa.nc","r") ; open netcdf
> file
> ;u = ufile->uwnd(:,0,:,:) ; pull u off file
> ;vfile = addfile("vwind_may_2003_700hPa.nc","r")
> ;v = vfile->vwnd(:,0,:,:) ; pull v off file
>
> diri = "./"
> ufile = addfile(diri+"uwnd.2008.nc","r") ; open netcdf file
> u = short2flt(ufile->uwnd(:,{700},:,:)) ; pull u off file
> vfile = addfile(diri+"vwnd.2008.nc","r") ; open netcdf file
> v = short2flt(vfile->vwnd(:,{700},:,:)) ; pull v off file
> ;************************************************
> ; Calculate vorticity on a Rectilinear Grid
> ; Since this is a globalg grid, use spherical harmonics ["highly accurate"]
> ;************************************************
> vr = uv2vrF_Wrap (u,v) ; globe
> printVarSummary(vr)
> ;************************************************
> ; Extract vr over region
> ;************************************************
> latS = 7.5
> latN = 27.5
> lonL = 11.0
> lonR = 135.0
>
> vr_region = vr(:,{latS:latN},{lonL:lonR}) ;-- select the data of
> the sub-region
> printVarSummary(vr_region)
> printMinMax(vr_region,0)
>
> ;************************************************
> ; Find date where threshold is exceeded
> ;************************************************
> threshold = 5e-5 ;-- define threshold
>
> vr_thres = dim_max_n(vr_region, (/1,2/)) ;-- max at each time step
> nt_thres = ind(vr_thres .ge. threshold) ;-- time indices >=
> threshold
>
> if (ismissing(nt_thres(0))) then
> print("++++++++++++++")
> print(" NO values exceed threshold="+threshold)
> print("++++++++++++++")
> exit
> end if
>
> ;***********************************************
> ; print
> ;***********************************************
> ymd = cd_calendar(u&time, -2) ;-- all dates
> print(nt_thres+" "+ymd(nt_thres))
>
> ;***********************************************
> ; Save to a netcdf file
> ; Note that the indices are non-contiguous
> ;***********************************************
> diro = "./"
> filo = "vort_700hPa.threshold.nc"
> ptho = diro+filo
> system("/bin/rm -f "+ptho)
> ncdf = addfile(ptho,"c")
>
> fAtt = True
> fAtt at title = "Vorticity at 700hPa >= threshold="+threshold
> 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(nt_thres,:,:) ; just the subset dates
>
>
> ===
> OUTPUT
> ===
>
> (0) 130 20080510
> (1) 131 20080511
> (2) 207 20080726
> (3) 208 20080727
> (4) 234 20080822
> (5) 254 20080911
> (6) 255 20080912
> (7) 264 20080921
> (8) 265 20080922
> (9) 269 20080926
> (10) 270 20080927
> (11) 271 20080928
> (12) 275 20081002
> (13) 333 20081129
> (14) 349 20081215
> --
>
> On Wed, Apr 26, 2017 at 11:14 AM, Alan Brammer <abrammer at albany.edu>
> wrote:
>
>> 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 listncl-talk at ucar.edu
>> List instructions, subscriber options, unsubscribe: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
>>
>>
>>
>> _______________________________________________
>> 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/20170427/977ce16c/attachment-0001.html
More information about the ncl-talk
mailing list