[ncl-talk] Extracting "dates" if a threshold is exceeded in any gridpoint of a domain
Dennis Shea
shea at ucar.edu
Wed Apr 26 11:54:53 MDT 2017
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/20170426/07ea4361/attachment.html
More information about the ncl-talk
mailing list