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

Sorry, I thought of a better solution if you do not need to have the
numeric time coordinate values [39460..39795] etc. for a later purpose.
 Just compute YYYYMM inside the loop, and the time rebasing problem is
automatically solved:

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

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

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

Add this:
   YYYYMM(ny,:) = cd_calendar(xs&$xs!0$, 1)

3. After the loop, add this. YYYYMM is a 2-D array, NOT 4-D:
   print(YYYYMM(0:2,0:2))

--Dave

On Thu, Nov 7, 2013 at 11:36 AM, Dave Allured - NOAA Affiliate <
dave.allured@noaa.gov> wrote:

> 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:48:25 2013