Re: use cd_calendar when it is necessary to read an variable from many directories

Yi-Chih,

Yes, you are overwriting the previous time coordinates during array
concatenation, xx=xs(ny,:,:,:)= etc. Preserving the right metadata can be
tricky when concatenating source arrays, each with their own metadata.
 Save the time coordinates from each file in a separate array, like this:

1. After this line:
   xx= new((/nyr,ds(0),ds(1),ds(2)/),float,"No_FillValue")

Add this:
   xtime = new ( (/nyr,ds(0)/), double)

2. In the read loop, after this line:
   xx(ny,:,:,:)= xs

Add this. The critical metadata on the time coordinate variable is also
copied:
   xtime(ny,:) = xs&$xs!0$

3. After the loop, add this. YMDHMS becomes a 2-D array, NOT 4-D:
   YMDHMS = cd_calendar(xtime, 1)
   print(YMDHMS(0:2,0:2))

Note, this method require that the "units" attribute of the time
coordinates be identical between all files. In your case I think this is a
good assumption. "hours since 1800-01-01 00:00" is a commonly used fixed
value. If this is not the case, then you can add a little more code to
perform time rebasing.

--Dave

On Wed, Nov 6, 2013 at 10:00 PM, Yi-Chih Huang <dscpln@gmail.com> wrote:

> Dave,
>
> Currently my working environment is one variable and one year per
> directory. I am new to NCL so that I can't change the data for now
> although I have a feeling it would be easier to work on one variable per
> file. My input data and the script I am developing are as below.
>
> Thanks much,
>
> Yi-Chih
>
> #####
> Copyright (C) 1995-2013 - All Rights Reserved
> University Corporation for Atmospheric Research
> NCAR Command Language Version 6.1.2
> The use of this software is governed by a License Agreement.
> See http://www.ncl.ucar.edu/ for more details.
>
> Variable: f
> Type: file
> filename: ipvhnl.gdas.2008
> path: ipvhnl.gdas.2008.grib2
> file global attributes:
> dimensions:
> initial_time0_hours = 12
> lat_0 = 361
> lon_0 = 720
> variables:
> float TMP_P8_L1_GLL0 ( initial_time0_hours, lat_0, lon_0 )
> center : US National Weather Service - NCEP (WMC)
> production_status : Operational products
> long_name : Temperature
> units : K
> _FillValue : 1e+20
> grid_type : Latitude/longitude
> parameter_discipline_and_category : Meteorological products,
> Temperature
> parameter_template_discipline_category_number : ( 8, 0, 0,
> 0 )
> level_type : Ground or water surface
> level : 0
> type_of_statistical_processing : Average of N uninitialized
> analyses, starting at reference time.
> statistical_process_duration : not applicable (intervals of 6
> hours)
> number_in_average : <ARRAY of 12 elements>
> forecast_time : 0
> forecast_time_units : hours
>
> double initial_time0_hours ( initial_time0_hours )
> long_name : initial time
> units : hours since 1800-01-01 00:00
>
> double initial_time0_encoded ( initial_time0_hours )
> long_name : initial time encoded as double
> units : yyyymmddhh.hh_frac
>
> float lat_0 ( lat_0 )
> long_name : latitude
> grid_type : Latitude/Longitude
> units : degrees_north
> Dj : 0.5
> Di : 0.5
> Lo2 : 359.5
> La2 : -90
> Lo1 : 0
> La1 : 90
>
> float lon_0 ( lon_0 )
> long_name : longitude
> grid_type : Latitude/Longitude
> units : degrees_east
> Dj : 0.5
> Di : 0.5
> Lo2 : 359.5
> La2 : -90
> Lo1 : 0
> La1 : 90
>
> string initial_time0 ( initial_time0_hours )
> long_name : Initial time of first record
> units : mm/dd/yyyy (hh:mm)
>
>
> #####
> load "$SysE/lib/ncl/helper_libs.ncl"
>
> begin
> latS = -5.0
> latN = 5.0
> lonE = 120.0
> lonW = 280.0
>
> yrStrt = 1979
> yrLast = 2008
>
> dir1= "/fs3/SysE_DB/nmm/CFSR/Monthly/Mean/"
> dir2=
> (/"1979/","1980/","1981/","1982/","1983/","1984/","1985/","1986/","1987/","1988/","1989/",\
>
> "1990/","1991/","1992/","1993/","1994/","1995/","1996/","1997/","1998/","1999/","2000/",\
> "2001/","2002/","2003/","2004/","2005/","2006/","2007/","2008/"/)
> yr= systemfunc("cd "+dir1+" ; ls ") ; input years
> print(yr)
> varName= (/"sst"/)
> nyr= dimsizes(dir2)
> print(nyr)
>
> dir = str_concat(dir1+dir2(0))
> fil= systemfunc("cd "+dir+" ; ls sst.nc") ;;;;;;;; read SST
> f= addfile (dir+fil, "r")
> x= f->$varName(0)$
> dims = getfilevardims(f,varName(0))
> xs= x(:,{latS:latN},{lonE:lonW}) ; tropical Pacific
> ds= dimsizes(xs)
> xx= new((/nyr,ds(0),ds(1),ds(2)/),float,"No_FillValue")
>
> delete (fil)
> do ny= 0,nyr-1
> ;print(ny)
> ;print(dir2(ny))
> dir = str_concat(dir1+dir2(ny))
> fil= systemfunc("cd "+dir+" ; ls sst.nc")
> f= addfile (dir+fil, "r")
>
> x= f->$varName(0)$
> dims = getfilevardims(f,varName(0))
> xs= x(:,{latS:latN},{lonE:lonW}) ; tropical Pacific
> xx(ny,:,:,:)= xs
> end do
>
> xx!0= "year"
> xx!1= "month"
> printVarSummary(xx)
> ; YMDHMS = cd_calendar(xx, 1)
> ;print(YMDHMS)
>
> ; in = addfile("/fs3/yhuang/precip.mon.mean.nc","r") ; 197901-201306
> ; YMDHMS = cd_calendar(in->time, 0) ; all times
> ; index for desired time
> ; iYear = ind(YMDHMS(:,0).ge.yrStrt .and. YMDHMS(:,0).le.yrLast)
> ; read time and region
> precip = in->precip(iYear,{latS:latN},{lonE:lonW})
> pcpClm = clmMonTLL( precip )
>
> xt= precip(0:23,{latS:latN},{lonE:lonW})
> ds= dimsizes(xt)
> nyr= yrLast -yrStrt +1
> pcpm= new((/ds(0),ds(2)/),float)
>
> yyyymmM = cd_calendar(precip&time,-1) ; verify times are correct
>
> dimp = dimsizes(precip)
> ntim = dimp(0)
> nlat = dimp(1)
> mlon = dimp(2)
>
> pcpt = dim_avg_n_Wrap( pcpClm,1 ) ; (time,lat)
> do nm= 0,11
> pcpm(nm,:) = pcpt(nm,:)
> pcpm(nm+12,:)= pcpt(nm,:)
> end do
>
>
>
> On Thu, Nov 7, 2013 at 1:07 PM, Dave Allured - NOAA Affiliate <
> dave.allured@noaa.gov> wrote:
>
>> Yi-Chih,
>>
>> It seems that you have concatenated your 30 input directories in an
>> unusual way. You created two time dimensions (year and month), where most
>> existing NCL examples use a single time dimension ("time") that encodes
>> complete dates. The two dimensional case can also be handled, but we need
>> to better understand your input data.
>>
>> If you are still interested in getting cd_calendar to work, please show
>> the complete metadata for one of the original input files. Use ncdump -h
>> on the OS command line, or else printVarSummary(f) where f is the result of
>> addfile in NCL. Showing only xx is not enough, because I think there
>> are already mistakes in how this array was created. It looks like you
>> may accidentally be overwriting the date/time information in memory, for
>> all input files except the last one.
>>
>> Also please show the complete original NCL script that "reads one
>> variable xx from 30 directories for 30 years". I think this is where some
>> corrections will be needed.
>>
>> --Dave
>>
>>
>> On Wed, Nov 6, 2013 at 7:24 PM, Yi-Chih Huang <dscpln@gmail.com> wrote:
>>
>>> Dennis,
>>>
>>> Thanks much for the reply. I am trying to make the following
>>> commands working when the variable has to be created from 30 directories.
>>>
>>> YMDHMS = cd_calendar(in->time, 0) ; all times
>>> iYear = ind(YMDHMS(:,0).ge.yrStrt .and. YMDHMS(:,0).le.yrLast)
>>>
>>> The xx is created in the way below. The information month:
>>> [39460..39795] is from the months in 2008. If it is not possible to use
>>> cd_calendar when the variable has to be created from 30 directories, maybe
>>> I should give up cd_calendar and write the script in my own way.
>>>
>>> xs= x(:,{latS:latN},{lonE:lonW}) ; tropical Pacific
>>> xx(ny,:,:,:)= xs
>>>
>>> Many thanks,
>>>
>>> Yi-Chih
>>>
>>>
>>> On Thu, Nov 7, 2013 at 10:35 AM, Dennis Shea <shea@ucar.edu> wrote:
>>>
>>>>
>>>>
>>>> Hi Yi-Chih
>>>>
>>>> Please do not send direct emails to me or Mary.
>>>> Always send to ncl-talk.
>>>>
>>>> In particular, I am not paid to answer ncl-talk questions.
>>>> I do it voluntarily to help people. However, getting personal
>>>> emails makes the communication one-on-one. Now,
>>>> your problem is my problem. I just don't have the time,
>>>> I have my own job to do. Sending to ncl-talk means
>>>> many people see the issue and, hopefully, will respomd.
>>>>
>>>> The honest truth is that I and others are not sure what
>>>> your objective is.
>>>>
>>>> ---
>>>> Naming the dimensions is not important.
>>>>
>>>> **It is the *units* of the 'time' variable that are important**
>>>>
>>>>
>>>> The documentation states that the 'time' variable
>>>>
>>>> "Must have a "units" attribute string in the format
>>>> "units since basetime", for example, "days since 1971-1-1"
>>>>
>>>> time is usually a 'one dimensional array', not 2D (year,month)
>>>>
>>>> Normal netCDF style usage is
>>>>
>>>> xx(time,lat,lon) ===> xx(360,21,321)
>>>>
>>>> where the time units would be something like
>>>>
>>>> "days since ???? ...."
>>>>
>>>> ===
>>>> Perhaps an alternative would be ...
>>>>
>>>> You created the variable. You must know the years used.
>>>> If you know the start and end year
>>>>
>>>> yrStrt =
>>>> yrLast =
>>>> year = ispan(yrStrt,yrLast,1)
>>>> year!0 = "year"
>>>>
>>>> month = ispan(1,12,1)
>>>> month!0 = "month"
>>>>
>>>> xx!0 = "year"
>>>> xx&year = year
>>>>
>>>> xx!1 = "month"
>>>> xx&month= month
>>>>
>>>> printVarSummary(xx)
>>>>
>>>>
>>>>
>>>>
>>>>
>>>> ---
>>>>
>>>> > Variable: xx
>>>> > Coordinates:
>>>> > month: [39460..39795]
>>>> lon: [120..280]
>>>>
>>>> ==============
>>>>
>>>> (39795-39460+1) = 336 ..... not 360 so not wure what this represents.
>>>>
>>>>
>>>> On 11/6/13 5:58 PM, Yi-Chih Huang wrote:
>>>>
>>>>> Dennis,
>>>>>
>>>>> I gave the names of year and month to the first two dimensions of
>>>>> xx.
>>>>> But there is no change to YMDHMS. How to give coordinates to xx to
>>>>> make cd_calendar
>>>>> working?
>>>>>
>>>>> Thanks much,
>>>>>
>>>>> Yi-Chih
>>>>>
>>>>> ######
>>>>> xx!0= "year"
>>>>> xx!1= "month"
>>>>> printVarSummary(xx)
>>>>> YMDHMS = cd_calendar(xx, 1)
>>>>> print(YMDHMS)
>>>>>
>>>>> ######
>>>>> Variable: xx
>>>>> Type: float
>>>>> Total Size: 9707040 bytes
>>>>> 2426760 values
>>>>> Number of Dimensions: 4
>>>>> Dimensions and sizes: [year | 30] x [month | 12] x [lat | 21] x [lon
>>>>> |
>>>>> 321]
>>>>> Coordinates:
>>>>> month: [39460..39795]
>>>>> lat: [-5.. 5]
>>>>> lon: [120..280]
>>>>> Number Of Attributes: 16
>>>>> DataFreq : Monthly
>>>>> center : US National Weather Service - NCEP (WMC)
>>>>> production_status : Operational products
>>>>> long_name : Temperature
>>>>> units : K
>>>>> grid_type : Latitude/longitude
>>>>> parameter_discipline_and_category : Meteorological products,
>>>>> Temperature
>>>>> parameter_template_discipline_category_number : ( 8, 0, 0, 0
>>>>> )
>>>>> level_type : Ground or water surface
>>>>> level : 0
>>>>> type_of_statistical_processing : Average of N uninitialized
>>>>> analyses, starting at reference time.
>>>>> statistical_process_duration : not applicable (intervals of 6
>>>>> hours)
>>>>> number_in_average : <ARRAY of 12 elements>
>>>>> forecast_time : 0
>>>>> forecast_time_units : hours
>>>>> _FillValue : 1e+20
>>>>>
>>>>>
>>>>> Variable: YMDHMS
>>>>> Type: double
>>>>> Total Size: 19414080 bytes
>>>>> 2426760 values
>>>>> Number of Dimensions: 4
>>>>> Dimensions and sizes: [30] x [12] x [21] x [321]
>>>>> Coordinates:
>>>>> Number Of Attributes: 1
>>>>> _FillValue : 9.969209968386869e+36
>>>>> (0,0,0,0) 9.969209968386869e+36
>>>>>
>>>>>
>>>>>
>>>>> On Wed, Nov 6, 2013 at 10:40 PM, Dennis Shea <shea@ucar.edu> wrote:
>>>>>
>>>>> Your 'xx' is a standard variable. It is not a 'time' variable as
>>>>>> required by the 'cd_calendar' function:
>>>>>> http://www.ncl.ucar.edu/Document/Functions/Built-in/
>>>>>> cd_calendar.shtml
>>>>>>
>>>>>> The documentation states that the 'time' variable
>>>>>>
>>>>>> "Must have a "units" attribute string in the format
>>>>>> "units since basetime", for example, "days since 1971-1-1"
>>>>>>
>>>>>> Hence, your usage
>>>>>> YMDHMS = cd_calendar(xx, 1)
>>>>>> is not correct.
>>>>>>
>>>>>> You should do a printVarSummary(xx)
>>>>>>
>>>>>> Is there coordinate information for the named dimensions 'year'
>>>>>> and 'month'? If there is no information, you must create.
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> On 11/5/13 11:46 PM, Yi-Chih Huang wrote:
>>>>>>
>>>>>> Hello,
>>>>>>>
>>>>>>> I read one variable xx from 30 directories for 30 years. The
>>>>>>> dimension
>>>>>>> of the variable xx is (year,month,lat,lon). I am trying to use
>>>>>>> cd_calendar
>>>>>>> to get the information of year and month. However, every element of
>>>>>>> YMDHMS
>>>>>>> is 9.969209968386869e+36. In addition, the error below occurred
>>>>>>> when the
>>>>>>> command "iYear = ind(YMDHMS(:,0).ge.yrStrt .and.
>>>>>>> YMDHMS(:,0).le.yrLast)" was executed. What is the right way to use
>>>>>>> cd_calendar when it is necessary to read an variable from many
>>>>>>> directories?
>>>>>>>
>>>>>>> YMDHMS = cd_calendar(xx, 1)
>>>>>>>
>>>>>>> (29,11,20,320) 9.969209968386869e+36
>>>>>>> fatal:Number of subscripts do not match number of dimensions of
>>>>>>> variable,(2) Subscripts used, (4) Subscripts expected
>>>>>>> ^Mfatal:["Execute.c":8128]:Execute: Error occurred at or near line
>>>>>>> 51 in
>>>>>>> file sstTime.ncl
>>>>>>>
>>>>>>> Thanks much,
>>>>>>>
>>>>>>> Yi-Chih
>>>>>>>
>>>>>>

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Nov 7 11:36:51 2013