<div dir="ltr">Dear all,<div><br></div><div><br></div><div>Many many thanks for the help!</div><div><br></div><div><br></div><div>Sincerely,<br><div class="gmail_extra"><br clear="all"><div><div class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><font face="comic sans ms, sans-serif"><b><u>Lyndz</u></b></font></div><div dir="ltr"><div><br></div></div></div></div></div></div></div></div></div></div>
<br><div class="gmail_quote">On Thu, Apr 27, 2017 at 2:54 AM, Dennis Shea <span dir="ltr"><<a href="mailto:shea@ucar.edu" target="_blank">shea@ucar.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div>Well, I guess 'great minds think alike' .... <br>Alan ... What r u doing on my turf????? :-)<br><br></div>I had basically used Alan's approach. I have some local files which had daily data for one year.<br></div>The threshold=3e-5 yielded many days exceeding this value. For the example. I used 5e-5<br><br>=======<br></div><div>SCRIPT<br>=======<span class=""><br>******************************<wbr>******************<br>; variable and file handling<br>;*****************************<wbr>*******************<br> ;ufile = addfile("uwind_may_2003_<wbr>700hPa.nc","r") ; open netcdf file<br> ;u = ufile->uwnd(:,0,:,:) <wbr> ; pull u off file<br></span> ;vfile = addfile("vwind_may_2003_<wbr>700hPa.nc","r")<span class=""><br> ;v = vfile->vwnd(:,0,:,:) <wbr> ; pull v off file<br><br></span> diri = "./"<br> ufile = addfile(diri+"<a href="http://uwnd.2008.nc" target="_blank">uwnd.2008.nc</a>","<wbr>r") ; open netcdf file<br> u = short2flt(ufile->uwnd(:,{700},<wbr>:,:)) ; pull u off file<br> vfile = addfile(diri+"<a href="http://vwnd.2008.nc" target="_blank">vwnd.2008.nc</a>","<wbr>r") ; open netcdf file<br> v = short2flt(vfile->vwnd(:,{700},<wbr>:,:)) ; pull v off file<br>;*****************************<wbr>*******************<br>; Calculate vorticity on a Rectilinear Grid<br>; Since this is a globalg grid, use spherical harmonics ["highly accurate"]<br>;*****************************<wbr>******************* <br> vr = uv2vrF_Wrap (u,v) ; globe<br> printVarSummary(vr)<br>;*****************************<wbr>*******************<br>; Extract vr over region<br>;*****************************<wbr>******************* <br> latS = 7.5<br> latN = 27.5<br> lonL = 11.0<br> lonR = 135.0<br><br> vr_region = vr(:,{latS:latN},{lonL:lonR}) <wbr> ;-- select the data of the sub-region<br> printVarSummary(vr_region)<br> printMinMax(vr_region,0)<br><br>;*****************************<wbr>*******************<br>; Find date where threshold is exceeded<br>;*****************************<wbr>******************* <br> threshold = 5e-5 <wbr> ;-- define threshold<br><br> vr_thres = dim_max_n(vr_region, (/1,2/)) ;-- max at each time step <br> nt_thres = ind(vr_thres .ge. threshold) ;-- time indices >= threshold<br><br> if (ismissing(nt_thres(0))) then<br> print("++++++++++++++")<br> print(" NO values exceed threshold="+threshold)<br> print("++++++++++++++")<br> exit<br> end if<br><br>;*****************************<wbr>******************<br>; print<br>;*****************************<wbr>******************<br> ymd = cd_calendar(u&time, -2) ;-- all dates<br> print(nt_thres+" "+ymd(nt_thres))<br><br>;*****************************<wbr>******************<br>; Save to a netcdf file<br>; Note that the indices are non-contiguous<br>;*****************************<wbr>******************<br> diro = "./"<br> filo = "<a href="http://vort_700hPa.threshold.nc" target="_blank">vort_700hPa.threshold.nc</a>"<br> ptho = diro+filo<br> system("/bin/rm -f "+ptho)<br> ncdf = addfile(ptho,"c")<br><br> fAtt = True<br> fAtt@title = "Vorticity at 700hPa >= threshold="+threshold<span class=""><br> fAtt@source_file = "daily winds from NCEP-NCAR Reanaysis"<br> fAtt@Conventions = "None"<br> fAtt@creation_date = systemfunc("date")<br> fileattdef(ncdf,fAtt) <wbr> ; copy file attributes<br> filedimdef(ncdf,"time",-1,<wbr>True)<br></span> ncdf->vr = vr(nt_thres,:,:) ; just the subset dates<br><br></div><div><br>===<br></div>OUTPUT<br>===<br><br><div>(0) 130 20080510<br>(1) 131 20080511<br>(2) 207 20080726<br>(3) 208 20080727<br>(4) 234 20080822<br>(5) 254 20080911<br>(6) 255 20080912<br>(7) 264 20080921<br>(8) 265 20080922<br>(9) 269 20080926<br>(10) 270 20080927<br>(11) 271 20080928<br>(12) 275 20081002<br>(13) 333 20081129<br>(14) 349 20081215<br><div><div>--<br></div></div></div></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Apr 26, 2017 at 11:14 AM, Alan Brammer <span dir="ltr"><<a href="mailto:abrammer@albany.edu" target="_blank">abrammer@albany.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word">Slight alternative to Karen's suggestion. Same basic idea but without the loop in the script it self.<div><br><div><br></div><div><div>load "$NCARG_ROOT/lib/ncarg/nclscri<wbr>pts/contrib/cd_string.ncl"</div><span><div>begin</div><div><br></div><div> f = addfile("vort_700hPa.nc","r") ;-- open file</div><div><br></div><div> vr = f->vr ;-- read variable</div><div> </div><div> time = f->time ;-- read time</div><div> ntimes = dimsizes(time) ;-- get size of time</div><div><br></div><div>;-- define sub-region</div><div> lat0 = 7.5</div><div> lat1 = 27.5</div><div> lon0 = 11.0</div><div> lon1 = 135.0</div><div> </div><div> vr_region = vr(:,{lat0:lat1},{lon0:lon1}) ;-- select the data of the sub-region</div><div> printVarSummary(vr_region)</div><div><br></div><div> threshold = 3e-5 ;-- define threshold</div><div><br></div></span><div> vr_max = dim_max_n(vr_region, (/1,2/)) ;; max over lat/lon per time. </div><div> time_inds = ind(vr_max .ge. threshold) ;; if the max isn’t greater than threshold then no values are. </div><div><br></div><div> if(all(ismissing(time_inds)))</div><div> print("No Times exceed threshold")</div><div> status_exit(1)</div><div> end if</div><div><br></div><div> times = time(time_inds)</div><div> print(cd_string(times, ""))</div><div>end</div><div><div class="m_-6461354879875006399h5"><div><br></div><div> </div><div><br></div><div><br></div><div><blockquote type="cite"><div>On 26 Apr 2017, at 12:20, Karin Meier-Fleischer <<a href="mailto:meier-fleischer@dkrz.de" target="_blank">meier-fleischer@dkrz.de</a>> wrote:</div><br class="m_-6461354879875006399m_-3120263377794784492Apple-interchange-newline"><div>
<div bgcolor="#FFFFFF" text="#000000">
Hi Lyndon,<br>
<br>
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.<br>
<br>
Your script should be looking like:<br>
<br>
<font color="#3333ff"><tt>begin</tt></font><br>
<br>
<font color="#3333ff"><tt> f =
addfile("vort_700hPa.nc","r") <wbr> ;-- open file</tt><tt><br>
</tt><tt><br>
</tt><tt> vr = f->vr <wbr> ;--
read variable</tt><tt><br>
</tt><tt> </tt><tt><br>
</tt><tt> time = f->time <wbr> ;--
read time</tt><tt><br>
</tt><tt> ntimes = dimsizes(time) <wbr> ;-- get
size of time</tt><tt><br>
</tt><tt><br>
</tt><tt>;-- define sub-region</tt><tt><br>
</tt><tt> lat0 = 7.5</tt><tt><br>
</tt><tt> lat1 = 27.5</tt><tt><br>
</tt><tt> lon0 = 11.0</tt><tt><br>
</tt><tt> lon1 = 135.0</tt><tt><br>
</tt><tt> </tt><tt><br>
</tt><tt> vr_region = vr(:,{lat0:lat1},{lon0:lon1}) <wbr> ;--
select the data of the sub-region</tt><tt><br>
</tt><tt> printVarSummary(vr_region)</tt><tt><br>
</tt><tt><br>
</tt><tt> threshold = 3e-5 <wbr> ;--
define threshold</tt><tt><br>
</tt><tt><br>
</tt><tt>;-- loop over time to retrieve times where data is
greater equal threshold</tt><tt><br>
</tt><tt> j = 0</tt><tt><br>
</tt><tt> do i=0,ntimes-1</tt><tt><br>
</tt><tt> tmp := ndtooned(vr_region(i,:,:)) <wbr> ;--
convert to 1d array</tt><tt><br>
</tt><tt> if(any(tmp .ge. threshold)) then</tt><tt><br>
</tt><tt> if(j .eq. 0) then</tt><tt><br>
</tt><tt> tt = time(i) <wbr> ;-- first
time step where threshold is exceeded</tt><tt><br>
</tt><tt> else</tt><tt><br>
</tt><tt> tt := array_append_record (tt, time(i), 0)
;-- append next time step </tt></font><font color="#3333ff"><tt><font color="#3333ff"><tt>where threshold is exceeded</tt></font></tt><tt><br>
</tt><tt> end if</tt><tt><br>
</tt><tt> j = j+1</tt><tt><br>
</tt><tt> end if </tt><tt><br>
</tt><tt> end do</tt><tt><br>
</tt><tt><br>
</tt><tt>;-- convert time to readable date string</tt><tt><br>
</tt><tt> utc_date = cd_calendar(tt,
0) ;-- convert date to UT-referenced
date</tt><tt><br>
</tt><tt> year =
sprinti("%0.4i",tointeger(utc_<wbr>date(:,0))) ;-- get year as
integer value</tt><tt><br>
</tt><tt> mon =
sprinti("%0.2i",tointeger(utc_<wbr>date(:,1))) ;-- get month as
integer value</tt><tt><br>
</tt><tt> day =
sprinti("%0.2i",tointeger(utc_<wbr>date(:,2))) ;-- get day as
integer value</tt><tt><br>
</tt><tt> str_date =
year+"/"+mon+"/"+day <wbr> ;-- yyyy/mm/dd</tt><tt><br>
</tt><tt> print(str_date)</tt></font><br>
<br>
<font color="#3333ff"><tt>end</tt></font><br>
<br>
Hope this helps,<br>
Karin<br>
<br>
<br>
<div class="m_-6461354879875006399m_-3120263377794784492moz-cite-prefix">Am 26.04.17 um 07:09 schrieb Lyndon
Mark Olaguera:<br>
</div>
<blockquote type="cite">
<div dir="ltr">Dear fellow NCL users,
<div><br>
</div>
<div><br>
</div>
<div>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.</div>
<div><br>
</div>
<div>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:</div>
<div> LAT : 7.5N to 27.5N</div>
<div> LON : 110E to 135E</div>
<div><br>
</div>
<div>the threshold value that I am using is 3e-5 per second.</div>
<div><br>
</div>
<div><br>
</div>
<div>I am attaching the netcdf file and the image in this
email.
<div class="gmail_chip gmail_drive_chip" style="width:396px;height:18px;max-height:18px;background-color:rgb(245,245,245);padding:5px;color:rgb(34,34,34);font-family:arial;font-style:normal;font-weight:bold;font-size:13px;border:1px solid rgb(221,221,221);line-height:1"><a href="https://drive.google.com/file/d/0B9faET7Bc2o8RWRXcTdKZEx0YlE/view?usp=drive_web" style="display:inline-block;max-width:366px;overflow:hidden;text-overflow:ellipsis;white-space:nowrap;text-decoration:none;padding:1px 0;border:none" target="_blank"><img style="vertical-align:bottom;border:none" src="https://ssl.gstatic.com/docs/doclist/images/icon_10_generic_list.png"> <span dir="ltr" style="color:rgb(17,85,204);text-decoration:none;vertical-align:bottom">test_2003_30.png</span></a><img style="opacity:0.55;float:right;display:none"></div>
</div>
<div><br>
</div>
<div>Is this possible in NCL?</div>
<div>Is there any function in NCL that can extract values from
every gridpoint in a specified domain?</div>
<div><br>
</div>
<div><br>
</div>
<div>I'll appreciate if anyone can point me to the right
direction/functions to implement this in NCL.</div>
<div><br>
</div>
<div>Many thanks in advance!</div>
<div><br>
</div>
<div><br>
</div>
<div>
<div>begin</div>
<div>;*****************************<wbr>*******************</div>
<div>; variable and file handling</div>
<div>;*****************************<wbr>*******************</div>
<div> ufile = addfile("uwind_may_2003_700hPa<wbr>.nc","r")
; open netcdf file</div>
<div> u = ufile->uwnd(:,0,:,:)
; pull u off file</div>
<div> lat = ufile->lat</div>
<div> lon = ufile->lon</div>
<div> vfile = addfile("vwind_may_2003_700hPa<wbr>.nc","r")</div>
<div> v = vfile->vwnd(:,0,:,:)
; pull v off file</div>
<div>;*****************************<wbr>*******************</div>
<div>; calculate vorticity on a Gaussian Grid</div>
<div>;*****************************<wbr>*******************
</div>
<div> ;vrt = u ;
retain coordinates</div>
<div> ;vrt = uv2vrG_Wrap(u,v)</div>
<div> vr = uv2vr_cfd (u,v,lat,lon, 3)</div>
<div> copy_VarMeta(u,vr)</div>
<div> vr@long_name = "vorticity"</div>
<div> vr@units = "per second"</div>
<div>;*****************************<wbr>******************</div>
<div>; save to a netcdf file</div>
<div>;*****************************<wbr>******************</div>
<div> ncdf = addfile("vort_700hPa.nc","c")</div>
<div> fAtt = True</div>
<div> fAtt@title = "Vorticity at 700hPa"</div>
<div> fAtt@source_file = "daily winds from NCEP-NCAR
Reanaysis"</div>
<div> fAtt@Conventions = "None"</div>
<div> fAtt@creation_date = systemfunc("date")</div>
<div> fileattdef(ncdf,fAtt) ; copy file attributes</div>
<div> filedimdef(ncdf,"time",-1,True<wbr>)</div>
<div> ncdf->vr = vr</div>
<div>end</div>
</div>
<div><br>
</div>
<div>
<div>
<div class="m_-6461354879875006399m_-3120263377794784492gmail_signature">
<div dir="ltr">
<div dir="ltr">
<div dir="ltr">
<div dir="ltr">
<div dir="ltr">
<div><font face="comic sans ms, sans-serif"><b><u>Lyndz</u></b></font></div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
<br>
<fieldset class="m_-6461354879875006399m_-3120263377794784492mimeAttachmentHeader"></fieldset>
<br>
<pre>______________________________<wbr>_________________
ncl-talk mailing list
<a class="m_-6461354879875006399m_-3120263377794784492moz-txt-link-abbreviated" href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a>
List instructions, subscriber options, unsubscribe:
<a class="m_-6461354879875006399m_-3120263377794784492moz-txt-link-freetext" href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" target="_blank">http://mailman.ucar.edu/mailma<wbr>n/listinfo/ncl-talk</a>
</pre>
</blockquote>
</div>
______________________________<wbr>_________________<br>ncl-talk mailing list<br><a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a><br>List instructions, subscriber options, unsubscribe:<br><a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" target="_blank">http://mailman.ucar.edu/mailma<wbr>n/listinfo/ncl-talk</a><br></div></blockquote></div><br></div></div></div></div></div><br>______________________________<wbr>_________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu" target="_blank">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/mailma<wbr>n/listinfo/ncl-talk</a><br>
<br></blockquote></div><br></div>
</div></div></blockquote></div><br></div></div></div>