Difference between revisions of "GAMAP tips and tricks"

Revision as of 20:45, 3 April 2008

Date and time

Julian Day vs. Astronomical Julian Day

The Julian Day (JD) is used to denote the number of days that have elapsed since the start of a year. For example, the Julian Day number of 2008/01/10 is 10 (since it is the 10th day of 2008).

However, another term that is sometimes confused with the Julian Day is the Astronomical Julian Day (AJD). This is the number of days that have elapsed since 12:00 GMT on January 1, 4713 BC. The Astronomical Julian Day at 00:00 GMT on 2008/01/10 is 2454475.5.

The Astronomical Julian Day is mostly used in astrophysics and space sciences to compute long intervals of time that may span many years (e.g. the period of a variable star, etc.). You can find algorithms for computing AJD in several textbooks, including Practical Astronomy with Your Calculator by Peter Duffett-Smith, Cambridge Univ. Press, 1992.

The chief advantage of using the Astronomical Julian Day is that its computation accounts for the end-of-month, end-of-year, end-of-century, and end-of-millenium transitions. This makes it very useful for computing future or past dates in scripts and computer programs.

For example, if you want to compute what the date will be 1000 days from 2008/01/01, then all you have to do is:

1. compute the Astronomical Juliay Day for 2008/01/01
2. add 1000
3. convert back to a calendar date

We'll return to this example a little further below.

IDL has two Astronomical Julian Day functions

• JULDAY -- converts a year, month, day (, hour, min, sec) to Astronomical Julian Date
• CALDAT -- converts an Astronomical Julian Day back to year, month, day (, hour, min, sec)

and they are called as follows:

; Compute the Astronomical Julian Day for 2008/01/10
; Note: the day is required, the hours, mins, seconds are optional
IDL> print, julday( 1, 10, 2008, 0, 0, 0 )
       2454475.5

; Convert the Astronomical Julian Day back to a calendar date
IDL> caldat, 2454475.5, y, m, d, h, mi, s
IDL> print, y, m, d, h, mi, s
           1          10        2008           0           0       0.0000000

To compute the date 1000 days after 2008/01/01 (from the example above), you would do the following:

IDL> jd = julday( 1, 1, 2008 ) 
IDL> jd = jd + 1000
IDL> caldat, jd, y, m, d
IDL> print, y, m, d
          9          27        2010

We find that gives us the date 2010/09/27. NOTE: This has also accounted for the leap-year-day on Feb 29, 2008.

For your convenience, GAMAP has a function called ADD_DATE which will do this for you in one fell swoop:

IDL> print, add_date( 20080101, 1000 )
    20100927

TIP: In order to avoid confusing Julian Day with Astronomical Julian Day, we recommend that you refer to the number of days in a year as "Day of Year". This will remove all ambiguity.

--Bob Yantosca 15:40, 3 April 2008 (EDT)

Computing the day of year from a calendar date (and vice-versa)

A quick way to compute the day of the year is w/ IDL's JULDAY function:

Day_of_Year = JULDAY( Month, Day, Year ) - JULDAY( 1, 0, Year )

For example:

IDL> print, julday( 1, 10, 2008 ) - julday( 1, 0, 2008 )
      10

JULDAY( 1, 0, 2008 ) is essentially 12/31/2007. This is used so that you don't have to add one to the above difference. (This is equivalent to JULDAY( 1, 10, 2008 ) - JULDAY( 1, 1, 2008 ) + 1.)

This is essentially done for you in the GAMAP function ~/IDL/gamap2/date_time/day_of_year.pro. You could replace the above code with:

IDL> print, day_of_year( 1, 10, 2008 )
      10

To compute the calendar date from the day of the year (the inverse operation), you can use IDL's CALDAT function. However, you must also add the Astronomical Julian day for the first of the year, for example:

IDL> caldat, julday( 1, 0, 2008 )+10, y, m, d
IDL> print, y, m, d
        1          10        2008

You can also do this with the GAMAP functions DAY_OF_YEAR and ADD_DATE, as follows

; Convert month/day/year to day of year
IDL> doy = day_of_year( 1, 10, 2008 )
IDL> print, doy
      10

; convert day of year back to month/day/year
IDL> print, add_date( 20080101, 10-1 )
    20080110

The only thing to remember is you have to subtract 1 from the day of year due to the way that ADD_DATE is written.

NOTE: Starting in GAMAP v2-12 (not released yet), you will be able to specify a single YYYYMMDD argument instead of the month, day, year arguments:

; Convert month/day/year to day of year
IDL> doy = day_of_year( 20080110 )

--Bob Yantosca 15:40, 3 April 2008 (EDT)

Combining output from timeseries files

Ray Nassar (ray@io.as.harvard.edu) wrote:

I just have a quick question, is there any GAMAP routine that can average all hourly data blocks in a timeseries file to make a single daily average bpch file?

It appears like I could use the average option of ctm_sum.pro but I do not quite understand where the averaged data goes and afterwards would still have to write to bpch.

Philippe Le Sager (plesager@seas.harvard.edu) replied:

Check the

   /gamap2/timeseries/gc_combine_nd49.pro
   /gamap2/timeseries/gc_combine_nd48.pro 

routines, which combines daily bpch files (nd48 & nd49 output, but also met field input files) into 4D data blocks. It has many options. You can extract a subset of data according to time and/or location, or process the data (moving average, daily max, shift to local time). You can either save the data into a new bpch file or just get an array of data in output.

I wrote a tutorial that gives some pointers here.

-Philippe

--Bob Yantosca 09:43, 1 April 2008 (EDT)

Memory management and CTM_MAKE_DATAINFO

When using CTM_MAKE_DATAINFO, you create few pointers and allocate memory to pointed data. Three scenarios to free that memory are possible.

(1) By default, GAMAP keeps track of the pointers in a global structure, and you can clean up the memory by calling CTM_CLEANUP (with or without the keyword /NO_GC):

    ctm_make_datainfo(data, ...)
    CTM_WriteBpch, DataInfo, FileInfo, FileName=file
    ctm_cleanup

(2) If you use the keyword /NO_GLOBAL when calling CTM_MAKE_DATAINFO, the story is a little bit more subtle. Now, GAMAP has no idea of the created pointers. If you use CTM_CLEANUP without the keyword /NO_GC, everything is clean up because there is a call to heap_gc. But you are also loosing **all other** pointers and objects.

    ctm_make_datainfo(data, ..., /No_global)
    CTM_WriteBpch, DataInfo, FileInfo, FileName=file
    ctm_cleanup

Note if you do not call ctm_cleanup or call it with /No_GC, then the memory allocated by CTM_MAKE_DATAINFO is still allocated and the pointers that refers to it are alive... until you exit the routine (unless they are passed back). Once you are out of the routine, this memory remains allocated but is useless since unaccessible, in other words you have memory leak.

(3) So, if you want to keep some objects and/or pointers alive in your code (i.e., you do not want to call heap_gc or ctm_cleanup,/no_gc), you need to free only the created pointers as follows:

    ctm_make_datainfo(data, datainfo, fileinfo, ..., /No_global)
    CTM_WriteBpch, DataInfo, FileInfo, FileName=file
    ptr_free, DataInfo.data
    ptr_free, fileinfo.gridinfo

To free only the heap memory created by multiple calls to CTM_MAKE_DATAINFO, the procedure is:

  for D=0L, NTracers-1L do begin
     
     Success = CTM_Make_DataInfo( Data[*,*,D], DataInfo, FileInfo, ..., /No_Global )
     
     ArrDataInfo = D eq 0l ? [ DataInfo ] : [ ArrDataInfo, DataInfo ]
     
     if D ne Ntracers-1l then ptr_free, Fileinfo.gridinfo
     
  endfor
     
  CTM_WriteBpch, ArrDataInfo, FileInfo, FileName=OutFileName
     
  for d=0, n_elements(ArrDataInfo)-1l do ptr_free, ArrDataInfo[d].data
  ptr_free, fileinfo.gridinfo

--Philippe Le Sager 10:28, 13 February 2008