[ncl-talk] climo_3.ncl + loop + conventional subscripts

zoe jacobs zoejacobs1990 at gmail.com
Fri Feb 28 14:39:22 MST 2020


Dear all,
considering the past conversation, I need to convert monthly
climatology (prcClm
  = clmMonTLL(prcMonth) ) to a nc file ( instead of plotting ). How can I
do this, please?
KInd regards,

On Mon, Nov 4, 2019 at 6:28 PM Dennis Shea <shea at ucar.edu> wrote:

> Hi Zoe
>
> This is offline. To avoid multiple emails, I downloaded two of the daily
> files [1979/1980].
>
> You will have to change
>
> ymdLast = 19801231
>
> to the appropriate value.
>
> ---
> You should look carefully at what is being done.
> Also, at the output and the documentation associated with graphic
> resources and functions..
> This information is not for me per se. It is for you to examine the data
> and utilities used.
>
> ---
> %> ncl prc.ncl_zoe
>
> Good Luck
>
>
> On Sat, Nov 2, 2019 at 11:26 PM zoe jacobs <zoejacobs1990 at gmail.com>
> wrote:
>
>> Dear Dennis,
>> Many thanks for your help. I used the attached script. I just want to
>> check with you to make sure that the script is correct. Would you please
>> check it?
>> Best wishes,
>>
>> On Mon, Oct 28, 2019 at 12:31 AM Dennis Shea <shea at ucar.edu> wrote:
>>
>>> The following is untested. However, it gives you the idea.
>>> It is your responsibility to look up documentation and look at the
>>> printed output.
>>>
>>>
>>>
>>>  fils = systemfunc ("ls precip.*.nc") ; file paths
>>>  a    = addfiles (fils, "r")
>>>  printVarSummary (a)
>>>
>>> ;************************************************
>>> ; Read the file
>>> ;************************************************
>>>
>>>   ymdStrt = 19790101                            ; climatology start
>>> year/month/day
>>>   ymdLast = 20151231                             ;      "      last
>>>  year/month/day
>>>   time    = a[:]->time
>>>   yyyymmdd= cd_calendar(time, -2)       ; note -2
>>> print(yyyymmdd)
>>>   ntStrt = ind(yyyymmdd.eq.ymdStrt)                ; start time index
>>>   ntLast = ind(yyyymmdd.eq.ymdLast)              ; last  time index
>>> print("ntStrt="+ntStrt+"  ntLast="+ntLast)
>>>    prc   = a[:]->precip(ntStrt:ntLast,:,:)               ; read precip
>>> from all files
>>>
>>>    printVarSummary (prc)
>>>    printMinMax (prc,0)
>>>   print("=====")
>>>                    ; you decide what you need
>>>    prcMonth = *calculate_monthly_values(*
>>> <http://www.ncl.ucar.edu/Document/Functions/Contributed/calculate_monthly_values.shtml>prc,
>>> "avg", 0, False)   ; monthly mean
>>>  ;;prcMonth = *calculate_monthly_values(*
>>> <http://www.ncl.ucar.edu/Document/Functions/Contributed/calculate_monthly_values.shtml>prc,
>>> "sum", 0, False)  ; monthly total
>>>
>>>    printVarSummary (prcMonth)
>>>    printMinMax (prcMonth,0)
>>>    print("=====")
>>>    opt = True                          ; examine data distribute
>>> [information only]
>>>    opt at PrintStat = True
>>>    stat_prc = stat_dispersion(prcMonth, opt )
>>>
>>>    prcClm   = *clmMonTLL*
>>> <http://www.ncl.ucar.edu/Document/Functions/Contributed/clmMonTLL.shtml>(prcMonth)
>>> ; Compute monthly climatology
>>>    printVarSummary (prcClm)
>>>    printMinMax (prcClm,0)
>>>
>>> ====
>>> Add the plot code.
>>> Look carefully at: http://www.ncl.ucar.edu/Applications/HiResPrc.shtml
>>>
>>> EG: The following may be inappropriate for your dataset but this is the
>>> idea.
>>>
>>>  res at cnLevelSelectionMode = "ExplicitLevels"
>>>  res at cnLevels             = (/ 0.01, 0.02, 0.04, 0.08, 0.16, \
>>>                                0.32, 0.64, 0.96/)
>>>  res at cnFillColors         = (/"white","cyan", "green","yellow",\
>>>                               "darkorange","red","magenta","purple",\
>>>                               "black"/)
>>>
>>>
>>> On Sun, Oct 27, 2019 at 2:48 PM Dennis Shea <shea at ucar.edu> wrote:
>>>
>>>>
>>>> [0]
>>>> As noted on ncl-talk many time the first rule of data processing is
>>>> *Look at your data!"
>>>> The temperature data were *monthly means* [K].
>>>> The precipitation data are* daily* totals [mm/day]. How do I know? I
>>>> looked at the README file at the URL. You should do that.
>>>>
>>>> [1]
>>>> a =* addfiles*
>>>> <http://www.ncl.ucar.edu/Document/Functions/Built-in/addfiles.shtml>
>>>> (fils, "r")
>>>>
>>>> You have the following which is appropriate for an *addfile*
>>>> <http://www.ncl.ucar.edu/Document/Functions/Built-in/addfile.shtml>
>>>> reference
>>>>        time   = a->time
>>>>
>>>> However, you are using *addfiles*
>>>> <http://www.ncl.ucar.edu/Document/Functions/Built-in/addfiles.shtml>
>>>> It should be:
>>>>        time   = a*[:]*->time   ; *[:]* syntax is appropriate for an
>>>> *addfiles* reference
>>>>
>>>> [2]
>>>> Did you read the documentation for *cd_calendar *
>>>> <http://www.ncl.ucar.edu/Document/Functions/Built-in/cd_calendar.shtml>as
>>>> I suggested? The whole point of using *cd_calendar* is to return
>>>> 'human-readable' start and end times from which the actual start/stop time
>>>> index values could be determined.. Your script uses the actual time values
>>>> [ 692496 , 1016808 ] mixed with human-readable yyyymm.  Further, you used
>>>>
>>>>     yyyymm = cd_calendar(time, *-1*)         ; -1 returns YYYYMM
>>>>
>>>> The precipitation files contain [as noted] daily values. You must use a
>>>> human readable time like YYYYMMDD [eg: 19790101]. The us
>>>>
>>>>     yyyymmdd = cd_calendar(time, *-2*)         ; -2 returns YYYYMMDD
>>>>
>>>> This converts the actual time values [eg:      692496] to
>>>> 'human-readable' yyyymmdd [year=>yyyy, month=>mm, day=dd>*
>>>> 692496=>197901*01] values.
>>>>
>>>>
>>>> ===
>>>> Also, your script will need to be changed to accommodate the
>>>> precipitation variable plot values
>>>> There are many examples of how to do this.
>>>>
>>>> eg: *http://www.ncl.ucar.edu/Applications/HiResPrc.shtml*
>>>> <http://www.ncl.ucar.edu/Applications/HiResPrc.shtml>
>>>> ===
>>>> In a previous email, you mentioned that you converted degrees Kelvin to
>>>> C.
>>>> It is best to self-document the  variable with its new units.
>>>>
>>>>    tmp = tmp-273.15
>>>>    tmp at units = "degC"  ; this will be automatically plotted
>>>> ===
>>>> Good Luck
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> On Sat, Oct 26, 2019 at 8:19 AM zoe jacobs <zoejacobs1990 at gmail.com>
>>>> wrote:
>>>>
>>>>> Dear all,
>>>>> I would like to have climatology precipitation like temperature
>>>>> climatology discussed through previous emails. However, there are
>>>>> precipitation files for each year from 1979- 2015 (
>>>>> ftp://ftp.cdc.noaa.gov/Datasets/cpc_global_precip/). In the first
>>>>> step, I tried to merge all those files as one, and then use Dennis script
>>>>> to plot a map. I faced some errors and need your advice please. the script
>>>>> I am running has been attached.Many thanks in advance,
>>>>> Best regardsn Sat, Oct 26, 2019 at 3:20 PM zoe jacobs <
>>>>> zoejacobs1990 at gmail.com> wrote:
>>>>>
>>>>>> I just noticed that I could use res at gsnRightString = " deg C" .
>>>>>> Once again, many thanks for your guide.
>>>>>> Best wishes,
>>>>>>
>>>>>>
>>>>>> On Fri, Oct 25, 2019 at 11:07 PM Dennis Shea <shea at ucar.edu> wrote:
>>>>>>
>>>>>>> I forgot to mention that plots can be done in 'portrait' [default]
>>>>>>> or 'landscape' mode.
>>>>>>>
>>>>>>> You should experiment to see which is appropriate for your needs.
>>>>>>>
>>>>>>> resP at gsnPaperOrientation = "landscape"        ; "portrait" is
>>>>>>> default
>>>>>>>
>>>>>>>   do nmo=0,11                                   ; loop over the
>>>>>>> months
>>>>>>>      res at gsnLeftString     = months(nmo)
>>>>>>>      plot(nmo) = gsn_csm_contour_map(wks,tmpClm(nmo,:,:), res)  ;
>>>>>>> create plot
>>>>>>>   end do
>>>>>>>
>>>>>>>   if (resP at gsnPaperOrientation .eq."landscape") then
>>>>>>>
>>>>>>> gsn_panel(wks,plot,(/3,4/),resP)                                  ; 3 rows
>>>>>>> x 4 columns
>>>>>>>   else
>>>>>>>       gsn_panel(wks,plot,(/4,3/),resP)            ; portrait
>>>>>>> ; 4 rows x 3 columns
>>>>>>>   end if
>>>>>>>
>>>>>>> On Fri, Oct 25, 2019 at 12:15 PM Dennis Shea <shea at ucar.edu> wrote:
>>>>>>>
>>>>>>>> Hello,
>>>>>>>>
>>>>>>>> We do not usually do this but this is offline from ncl-talk.
>>>>>>>>
>>>>>>>> [1]
>>>>>>>> Conventional subscripts start at 1 for Fortran and Matlab.  Hence,
>>>>>>>> for 12 months [nmos=12]
>>>>>>>>      Fortran: do nmo=1,nmos
>>>>>>>>      Matlab: for 1:nmos   or 1:1:nmos  [I am not a Matlab user so
>>>>>>>> this is a guess.]
>>>>>>>> NCL 'conventional' subscripts start at 0 [like C/C++/IDL/Python]
>>>>>>>>      NCL: do nmo=0,ntim-1
>>>>>>>>
>>>>>>>> [2] Dave illustrated 'coordinate subscripting' which uses the *{*
>>>>>>>> ...*}* syntax.
>>>>>>>>
>>>>>>>> The script I am attaching uses conventional subscripting. Just like
>>>>>>>> Matlab/Fortran except it uses a range fro 0 to 11 [12 elements].
>>>>>>>>
>>>>>>>> [3] Please carefully read the documentation for *clmMonTLL*
>>>>>>>> <http://www.ncl.ucar.edu/Document/Functions/Contributed/clmMonTLL.shtml>,
>>>>>>>> *cd_calendar*
>>>>>>>> <http://www.ncl.ucar.edu/Document/Functions/Built-in/cd_calendar.shtml>
>>>>>>>> ,* ind*
>>>>>>>> <http://www.ncl.ucar.edu/Document/Functions/Built-in/ind.shtml>
>>>>>>>>
>>>>>>>> [4] Look at the output from '*printVarSummary*
>>>>>>>> <http://www.ncl.ucar.edu/Document/Functions/Built-in/printVarSummary.shtml>'.
>>>>>>>> USE *printVarSummary *frequently.
>>>>>>>>      Note the dimensions and added attributes:
>>>>>>>>      time_op_ncl : Climatology: 29 years
>>>>>>>>      info : function clmMonTLL: contributed.ncl
>>>>>>>>
>>>>>>>> [5] NCL offer a large number of color tables [palettes]:
>>>>>>>>
>>>>>>>> *http://www.ncl.ucar.edu/Document/Graphics/color_table_gallery.shtml*
>>>>>>>> <http://www.ncl.ucar.edu/Document/Graphics/color_table_gallery.shtml>
>>>>>>>>      Using a color map that emphasizes features can be very useful.
>>>>>>>>
>>>>>>>> [6] You must invest the time to learn any new language.
>>>>>>>> Karin Meier-Fliesher and Michael Bottinger [DKRZ] wrote a wonderful
>>>>>>>> tutorial.
>>>>>>>>
>>>>>>>> *http://www.ncl.ucar.edu/Document/Manuals/NCL_User_Guide/*
>>>>>>>> <http://www.ncl.ucar.edu/Document/Manuals/NCL_User_Guide/>
>>>>>>>>
>>>>>>>> I suggest you read it.
>>>>>>>>
>>>>>>>> Good luck
>>>>>>>> ==============
>>>>>>>> *%>* ncl  zoe_jacobs.hgcn_cams.ncl
>>>>>>>>
>>>>>>>> will produce a png file.
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> On Fri, Oct 25, 2019 at 11:16 AM Dave Allured - NOAA Affiliate via
>>>>>>>> ncl-talk <ncl-talk at ucar.edu> wrote:
>>>>>>>>
>>>>>>>>> 1.  In your example, please notice that the month loop has a
>>>>>>>>> stride of 3.  Therefore it is making plots for only four months:  January,
>>>>>>>>> April, July, and October.
>>>>>>>>>
>>>>>>>>> If you want all 12 months, remove the stride.  Then make sure the
>>>>>>>>> array "plot" is dimensioned 12.  I think you can interchangeably use either
>>>>>>>>> a 1-D or 2-D plot array (i.e. 3*4), your choice.
>>>>>>>>>
>>>>>>>>> The secondary index "i" is used to write the four plots into
>>>>>>>>> positions 0, 1, 2, 3 in the "plot" graphics array.  If you are making 12
>>>>>>>>> plots, you do not really need to use a secondary index.
>>>>>>>>>
>>>>>>>>> Panel plotting uses a graphics array containing multiple plots.
>>>>>>>>> Please study basic examples and documentation for making panel plots.
>>>>>>>>>
>>>>>>>>> 2.  That data file has full coordinates.  Therefore you can use
>>>>>>>>> either conventional or coordinate subscripting, or mixed, your choice.
>>>>>>>>>
>>>>>>>>> I suggest using one of the date functions with coordinate
>>>>>>>>> subscripting, to index the time subset that you want.  Something like this:
>>>>>>>>>
>>>>>>>>>     time_units = f->air&time at units
>>>>>>>>>     time1 = cd_inv_calendar (year1, 1, 1, 0, 0, 0, time_units, 0)
>>>>>>>>>     time2 = cd_inv_calendar (year2, 12, 31, 23, 0, 0, time_units,
>>>>>>>>> 0)
>>>>>>>>>     air_subset = f->air({time1:time2},:,:)
>>>>>>>>>
>>>>>>>>> Use printVarSummary to ensure that the subset has the dimension
>>>>>>>>> sizes that you expect.  Then you can proceed to compute the climatology.
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> On Fri, Oct 25, 2019 at 10:14 AM zoe jacobs via ncl-talk <
>>>>>>>>> ncl-talk at ucar.edu> wrote:
>>>>>>>>>
>>>>>>>>>> Hi all,
>>>>>>>>>> regarding climo_3.ncl (
>>>>>>>>>> https://www.ncl.ucar.edu/Applications/Scripts/climo_3.ncl) ,
>>>>>>>>>> I have 2 questions:
>>>>>>>>>> 1. I need to plot all months on one panel (say 3*4), and cannot
>>>>>>>>>> understand the logic behind the below loop, which used in the climo_3.ncl
>>>>>>>>>> script :
>>>>>>>>>>
>>>>>>>>>> i = -1                                        ; Climatologies
>>>>>>>>>>   do nmo=0,11,3                                 ; loop over the months
>>>>>>>>>>      i = i+1
>>>>>>>>>>      res at gsnCenterString   = months(nmo)+":"+time(0)/100 +"-"+ time(ntim-1)/100
>>>>>>>>>>      plot(i) = gsn_csm_contour_map(wks,prcClm(nmo,:,:), res)  ; create plot
>>>>>>>>>>   end do
>>>>>>>>>>
>>>>>>>>>> How  does it  work??/
>>>>>>>>>>
>>>>>>>>>> 2. I would like to show climotology temperature from 1987- 2015. Data which I am using is air.mon.mean.nc  ( https://www.esrl.noaa.gov/psd/data/gridded/data.ghcncams.html)  and it is using conventional subscripts. I am not familiar with conventional subscripts. So how can I convert coordinate subscripting to conventional subscripts?!!
>>>>>>>>>>
>>>>>>>>>> Please kindly advice me .
>>>>>>>>>>
>>>>>>>>>> Many thanks in advance,
>>>>>>>>>>
>>>>>>>>>> Best regards,
>>>>>>>>>>
>>>>>>>>>> _______________________________________________
>>>>>>>>> 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/20200229/88c195fb/attachment.html>


More information about the ncl-talk mailing list