DaViT Tutorials

Installing and Running DaViT
How do I...
... get DaViT? We are discouraging use of these IDL routines as we are working toward moving to a python-based analysis package. If you need IDL support, please contact someone on our personnel page.
... install DaViT? After you have downloaded the tarball, follow the instructions in DaViT Installation Guide.pdf.
... run DaViT? Start the DaViT by opening a shell and typing path-to-davit/bin/davit. If you have your own set of IDL library procedures/functions that you would like to use in conjunction with DaViT, I recommend setting up a small shell script that changes into your working directory, add your IDL library path to the IDL_PATH system variable and then starts DaViT. This file could look like this: cd /home/lbnc/idl IDL_PATH=$IDL_PATH+':+/home/lbnc/idl/idl_lib' export IDL_PATH /davit/bin/davit Be sure not to change the system variable IDL_STARTUP in that startup script or in your .bashrc or equivalent. This will prevent DaViT from starting properly.
... use the RST together with DaViT? When installing the RST, it requires you to set a bunch of environment variables, including IDL_STARTUP. You might be worried that when you use DaViT that the RST will stop working. That is not so. The DaViT env_davit shell script does all the things the RST did for you. You will notice for example, that the very first thing DaViT does when starting is to compile all the RST IDL routines. Also, DaViT comes bundled with the RST so that you don't need to install the RST in order to use it.
... run DaViT at VT? All Linux machines in the SuperDARN group at Virginia Tech should be set up such that when you open a shell and enter /davit/bin/davit, DaViT should start. If not check the following: 1) Is there a local directory called /davit? Enter ls /davit. If it says "No such file or directory", contact us. Or even better, just walk into my office, that should be faster than emailing. 2) Is there an entry in /etc/fstab that mounts the /davit directory? Enter more /etc/fstab | grep davit. If nothing is printed, contact us - walk, office, you remember? 3) Is the /davit directory mounted? Enter more /etc/mtab | grep davit. If nothing is printed, contact us. 4) Is IDL installed? Enter idl and see if that starts up. If not, contact us. 5) Did you change the IDL_STARTUP system variable? 6) If you get compile errors thrown by IDL during the startup, your IDL_PATH variable might not be set properly. Enter env | grep IDL_PATH. If nothing is printed, contact us.
... get DaViT to access my data? DaViT's default settings rely on a specific structure to access radar data (though some flexibility is possible). Once you provide the path to your data directory, DaViT searches for data in 'YYYY/file_format/radar_code/YYYYMMDD.HHMM.radar_code.file_format[.bz2]', where YYYY is the 4-digit year, MM the 2-digit month, DD the 2-digit day, radar_code the 3-letter radar code (bks, wal...) and file_format the file format (fitacf, fitex, rawacf...).

Misc
How do I...
... find that function which name I forgot? DaViT includes a little function called pri. You can use it to tell you what functions there are and what their syntax is. So say you remember that there is a function to read radar fit data, but you don't remember the name. just enter pri,'fit' and DaViT will tell you all about the functions that have the word fit in them. Then you remember, right, it was called fitread but you have no idea what the syntax is. Just enter pri,'fitread' and behold the magic. Caveat: pri only works on routines written in IDL. If you provide just a part of the function name, like in the first example, only the already compiled IDL routines are found. However, pri will try and compile a function of the name that you are providing. So it will give you the syntax of fitread even if that function is not compiled yet.

Reading Radar Data
How do I...
... read fitACF/fitEX files? Use the rad_fit_read command. The syntax is rad_fit_read, yyyymmdd, 3letterradarcode [, time=[hhmm,hhmm] ] If no time span is given via the time keyword, data for the whole day is read. This command will read the fitACF/EX data found in the RAD_FIT_DATA_PATH set in path-to-davit/bin/davit. It will first look for zipped fitEX files, then fitEX files, then zipped fitACF and finally fitACF files. The data is read into a structure called rad_fit_data in the RAD_DATA_BLK common block. Type rad_fit_info to display some information about the data currently present in that variable. You can also give the command a filename as string (which can contain wildcards) of the file(s) you would like to read. See the documentation for details.
... read grd/grdEX files? Use the rad_grd_read command. The syntax is rad_grd_read, yyyymmdd [, /north ] [, /south ] If neither north nor south is set, data for the northern hemisphere is read. This command will read the grdEX/grd data found in the RAD_GRD_DATA_PATH set in path-to-davit/bin/davit. It will first look for zipped grdEX files, then grdEX files, then zipped grd and finally grd files. The data is read into the structure called rad_grd_data in the RAD_DATA_BLK common block. Type rad_grd_info to display some information about the data currently present in that variable. You can also give the command a filename as string (which can contain wildcards) of the file(s) you would like to read. See the documentation for details.
... read map/mapEX files? Use the rad_map_read command. The syntax is rad_map_read, yyyymmdd [, /north ] [, /south ] If neither north nor south is set, data for the northern hemisphere is read. This command will read the mapEX/map data found in the RAD_MAP_DATA_PATH set in path-to-davit/bin/davit. It will first look for zipped mapEX files, then mapEX files, then zipped map and finally map files. The data is read into the structure called rad_map_data in the RAD_DATA_BLK common block. Type rad_map_info to display some information about the data currently present in that variable. You can also give the command a filename as string (which can contain wildcards) of the file(s) you would like to read. See the documentation for details.

Plotting
How do I...
... plot a range-time plot of radar data? Use the rad_fit_plot_rti command. The syntax is rad_fit_plot_rti That will plot a range-time plot of the selected parameter for a given beam with title and colorbar. Of course, you need to read data before you can plot any. If you just want to plot the panel without a title and colorbar, use the rad_fit_plot_rti_panel command. The syntax is rad_fit_plot_rti_panel Suppose you want to plot Hankalsami spectral width data for March 3rd, 2003, 12:00-18:00 UT. Just enter: rad_fit_read, 20030303, 'han', time=[1200,1800] rad_fit_plot_rti, time=[1200,1800], param='width'
... redirect graphic output to a PostScript? Use the ps_open command to open and the ps_close command to close the PostScript. The syntax is ps_open, filename and ps_close command. One advantage of using these higher level procedures (as opposed to IDL's native SET_PLOT and DEVICE) is that ps_close will write a comment into the header of the PostScript file. This comment contains the name of the IDL procedure that created that plot - actually, it is a complete stack trace. That is very helpful if one day all your PostScript files are called asdf.ps and test.ps and you have no idea how you made this plots. You can then simply open the PostScript file in any text editor and look at these comments to remind yourself. Here's an example. Suppose you want to plot Hankalsami spectral width data for March 3rd, 2003, 12:00-18:00 UT. Just enter: rad_fit_read, 20030303, 'han', time=[1200,1800] ps_open, '20030303_han_width.ps' rad_fit_plot_rti, time=[1200,1800], param='width' ps_close
... plot a timeseries plot of radar data? Use the rad_fit_plot_tsr command. The syntax is rad_fit_plot_tsr That will plot a timeseries plot of the selected parameter for a given beam and gate. It will also plot a title and colorbar. Of course, you need to read data before you can plot any. If you just want to plot the panel without a title and colorbar, use the rad_fit_plot_tsr_panel command. The syntax is rad_fit_plot_tsr_panel Suppose you want to plot Hankalsami spectral width data for March 3rd, 2003, 12:00-18:00 UT, beam 4 and gate 25. Just enter: rad_fit_read, 20030303, 'han', time=[1200,1800] rad_fit_plot_tsr, time=[1200,1800], param='width', beam=4, gate=25
... plot a FFT of a timeseries of radar data? As so often in programming, there are a million ways to do this. But it all reaks down to this 1) Read the data 2) Extract the data for a given parameter, beam, and gate 3) Calculate the FFT using IDL's FFT function 4) Calculate the frequencies 5) Plot Here is some example code to accomplish this for Blackstone data on April 8th, 2009, beam 7. We're going to average over some gates (5-35) and resample the data: ; read data rad_fit_read, 20090408, 'bks', time=[200,1200] di = rad_fit_get_data_index() binds = where( (*rad_fit_data[di]).beam eq 7, bc ) juls = (*rad_fit_data[di]).juls[binds] raw_vel = (*rad_fit_data[di]).velocity[binds,5:35] ; get rid of "invalid" data ninds = where( abs(raw_vel) gt 500., gc ) if gc gt 0L then raw_vel[ninds] = !values.f_nan ; average over gates avg_vel = total(raw_vel, 2, /nan)/( total( finite(raw_vel), 2) > 1. ) ; beam 7 is the camping beam, hence the sampling time is not constant ; so we resample to 8 seconds dt = 8. ijuls = timegen( start=juls[0], final=juls[bc-1L], step=dt, unit='S') int_vel = interpol( avg_vel, juls, ijuls) ; do FFT nn = n_elements(ijuls) ; before fft-ing it's a good idea to detrend ; the data and use a window function fft_vel = hanning(nn)*detrend(int_vel) power = (abs(fft(fft_vel, -1))^2)[0:nn/2] ; and calculate frequencies (in mHz) freqs = findgen(nn/2)/(nn*dt)*1000. ; and plot plot, freqs, power, /xstyle, xrange=[0,10], /ylog, yrange=[1e-2,1e1]
... plot a fan/scan of radar data? Use the rad_fit_plot_scan command. The syntax is rad_fit_plot_scan That will overlay a scan of the selected parameter over a map in geomagnetic (AACGM) coordinates. It will also plot a title and colorbar. Of course, you need to read data before you can plot any. If you just want to plot the panel without a title and colorbar, use the rad_fit_plot_scan_panel command. The syntax is rad_fit_plot_scan_panel Suppose you want to plot Hankalsami spectral width data for March 3rd, 2003, 12:18 UT, on a geographic map. Just enter: v rad_fit_read, 20030303, 'han', time=[1200,1300] rad_fit_plot_scan, time=1218, param='width', coords='geog', /no_fill, scale=[0,50]
... plot multiple panels on one page? There's a million ways to do that! Which, of course, is not the answer you wanted. So here is a scenario: You would like to plot a panel of OMNI data on the top and then Blackstone RTI data from beam 4 and then Rankin Inlet data from beam 7 and 12. All data is from March 10th, 2008, 12:00 to 15:00 UT. Here's the code that does that: omn_read, 20080310 set_format, /portrait, /sardines clear_page omn_plot_panel, 1, 4, 0, 0, param='bz_gsm', date=20080310, time=[1200,1500] rad_fit_read, 20080310, 'bks', time=[1200,1500] rad_fit_plot_rti_panel, 1, 4, 0, 1, date=20080310, time=[1200,1500], beam=4 rad_fit_read, 20080310, 'rkn', time=[1200,1500] rad_fit_plot_rti_panel, 1, 4, 0, 2, date=20080310, time=[1200,1500], beam=7 rad_fit_plot_rti_panel, 1, 4, 0, 3, date=20080310, time=[1200,1500], beam=12, /last To quote Steve M. from his Go tutorial: As easy as falling off a bike. But alas, you realize there is no colorbar and no title and you wanted it in magnetic coordinates and not gate and you wanted to plot velocity, not power! Hence you do this: omn_read, 20080310 set_format, /portrait, /sardines set_parameter, 'velocity' set_coordinates, 'magn' clear_page omn_plot_panel, 1, 4, 0, 0, param='bz_gsm', date=20080310, time=[1200,1500], /bar rad_fit_read, 20080310, 'bks', time=[1200,1500] rad_fit_plot_rti_panel, 1, 4, 0, 1, date=20080310, time=[1200,1500], beam=4, /bar rad_fit_read, 20080310, 'rkn', time=[1200,1500] rad_fit_plot_rti_panel, 1, 4, 0, 2, date=20080310, time=[1200,1500], beam=7, /bar rad_fit_plot_rti_panel, 1, 4, 0, 3, date=20080310, time=[1200,1500], beam=12, /last, /bar plot_colorbar rad_fit_plot_title, 'Data on 20080310', ' ' Now you realize that the range of the y axis for the Rankin Inlet data is too small, so you use the following command for plotting the Rankin Inlet data: rad_fit_plot_rti_panel, 1, 4, 0, 2, date=20080310, time=[1200,1500], beam=7, /bar, yrange=[60,90] rad_fit_plot_rti_panel, 1, 4, 0, 3, date=20080310, time=[1200,1500], beam=12, /last, /bar, yrange=[60,90] I hope that answers some questions. For further pointers read the documentation for define_panel.
... plot multiple radar scans on one map? You achive that by plotting an empty map, reading the fit data for one radar, overlay the scan, read fit data for the next radar, overlay the scan... Here's some code that plots an empty map in MLT coordinates and overlays a scan from Blackstone and Fort Hays East. The scan is that closest to 20091207 at 13:32 UT: clear_page scan_jul = julday(12,7,2009, 13, 32) map_plot_panel, date=20091207, time=1332, coords='mlt' rad_fit_read, 20091207, 'bks', time=[1200,1400] rad_fit_overlay_scan, jul=scan_jul, coords='mlt' rad_fit_read, 20091207, 'fhe', time=[1200,1400] rad_fit_overlay_scan, jul=scan_jul, coords='mlt' But that isn't really nice because you can't really see anything. So lets loose the color of the continents, change the limits of the map and overlay the field-of-view of the radars: clear_page scan_jul = julday(12,7,2009, 13, 32) map_plot_panel, date=20091207, time=1332, coords='mlt', xrange=[0,41], yrange=[0,41], /no_fill rad_fit_read, 20091207, 'bks', time=[1200,1400] rad_fit_overlay_scan, jul=scan_jul, coords='mlt' overlay_fov, names='bks', jul=scan_jul, coords='mlt', /no_fill rad_fit_read, 20091207, 'fhe', time=[1200,1400] rad_fit_overlay_scan, jul=scan_jul, coords='mlt' overlay_fov, names='fhe', jul=scan_jul, coords='mlt', /no_fill
... overlay the data I downloaded from somewhere on a map? First of all, you need to have the position of the data you would like to overlay in some kind of coodinate system. Once you have that, we will do some coordinate transformations, plot a map using map_plot_panel and then overlay the data. The answers varies a little bit, depending on the coordinate system you would like to display the data in. Also, it depends in what coordinate system the data you wish to overlay is in. In this example, let's assume that your data is in geographic coordinates but you would like to plot the whole shebang in MLT. The approach is as follows: 1) Read your data and their positions. 2) Convert the positions of your data into MLT coordinates. 3) Convert the MLT coordinates into the map coordinate system used by map_plot_panel. 4) Plot the map using map_plot_panel. 5) Overlay your data using IDL's oplot, plots, polyfill or whatever. Say you have the magnetic footprint of some spacecraft in geographic latitude longitude and would like to overlay that on a MLT map, together with positions of some ground-based magnetometers and the field-of-view of the Fort Hays East radar. The date is March 12th, 1997 at 06:32 UT. Here's a script that does that, assuming you have the geographic footprint in the variables f_glat and f_glon: ; calculate the year and seconds since beginning ; of the year for the mlt routine jul = julday(3, 12, 1997, 6, 32) year = 1997 ut_sec = ( jul - julday(1, 1, 1997, 0, 0) )*86400.d nf = n_elements(f_lat) f_mlat = f_glat f_mlon = f_glon for i=0, nf-1 do begin   ; convert footprint into geomagnetic (AACGM) coordinates   ; using the RST cnvcoord IDL command   temp = cnvcoord(f_glat[i], f_glon[i], 1.)   ; convert magnetic longitude in MLT longitude   f_mlat[i] = temp[0]   f_mlon[i] = mlt(year, ut_sec, temp[1])   ; convert into map coordinates   ; use the MLT switch to indicate that the input   ; longitude is in hours MLT rather than degree,   ; which is the default   temp = calc_stereo_coords(f_mlat[i], f_mlon[i], /mlt)   f_mlat[i] = temp[0]   f_mlon[i] = temp[1] endfor map_plot_panel, jul=jul, coords='mlt' oplot, f_mlat, f_mlon, thick=4 ; overlay radar fov and gbm stations overlay_fov, name='fhe', jul=jul, coords='mlt', /no_fill overlay_gbm, ['fsim','snkq'], jul=jul, coords='mlt'
... get the footprints in the previous example from? You need the spacecraft positions in some coordinate system, say GSM. DaViT includes the THEMIS Data Analysis Software (TDAS). That software includes the Tsyganenko magnetic field models and an IDL routine called trace2iono. Its basic syntax is trace2iono, time_array, in_pos_array, out_foot_array, external_model=external_model,   in_coord=in_coord, km=km, out_coord=out_coord, par=par Input to the function is time_array and in_pos_array. time_array is in UTC seconds, i.e. seconds since January 1, 1970, 00:00:00.000. To convert the spacecraft position in GSM into footprints in GEO using the Tsyg 96 model with Dst=-10 nT, pdyn=2.0 nPa, By=0.0 nT and Bz=0.0 nT, simply enter trace2iono, time_array, in_pos_array, out_foot_array, external_model='t96',   in_coord='GSM', out_coord='GEO', par=[2.0,-10.,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] The array out_foot_array will now contain the three coordinates of the footprint, one for each time given in time_array. Note that the input parameter array for the model needs to always contain 10 entries. You can now look at the previous example to learn how to plots these footprints onto a map in MLT coordinates. Well actually, you first need to convert the footprints from the geocentric coordinates returned by trace2iono to the geocentric latitude and longitude system. But I trust you are able to figure that one out by yourself.
... plot a map, overlay THEMIS all-sky imager and radar data? DaViT includes the THEMIS Data Analysis Software (TDAS). That software allows you to download all space and ground-based data within the THEMIS mission. Say you want to plot the data in geomagnetic coordinates. So here's what we'll do: 1) Download the all-sky imager (ASI) data and its calibration 2) Read ASI data 3) Convert the positions of the field-of-view of the ASIs to geomagnetic and then map coordinates 4) Plot the data This script does just that: ; call some TDAS routines to get the right data ; for March 10th, 2008, 4:00 until 6:00 UT timespan, '2008-03-10:04:00:00.000', 2, /hours thm_load_asi, site='FSIM', /downloadonly thm_load_asi_cal, 'FSIM', /downloadonly ; this will download the data and place it somewhere ; on your hard disk. You can now either use TDAS's ; get_data to read the data into variables or you ; can locate the CDF files, use cdf_i to find the ; variable names in the CDF file and use cdf_read ; to load the data. I prefer the latter way datfile = '~/data/themis/thg/l1/asi/fsim/2008/03/thg_l1_asf_fsim_2008031004_v01.cdf' calfile = '~/data/themis/thg/l2/asi/cal/thg_l2_asc_fsim_19700101_v01.cdf' ; see what variables are in the CDF files dummy = cdf_i(datfile, zvar=zvar) print, zvar ; tells me I want datvars = ['thg_asf_fsim', 'thg_asf_fsim_time'] dat = cdf_read(datfile, datvars) ; and dummy = cdf_i(calfile, zvar=zvar) print, zvar calvars = ['thg_asf_fsim_elev','thg_asf_fsim_glon','thg_asf_fsim_glat'] cal = cdf_read(calfile, calvars) ; Now you have the brighness of the imager in dat.thg_asf_fsim, ; the corresponding times in dat.thg_asf_fsim_time and the ; locations of each pixal in cal.thg_asf_fsim_glon and ; cal.thg_asf_fsim_glat ; let convert that into magnetic coordinates and the ; map coordinates and then plot the map ; ; Set the scale for the brightness scale = [3e3,6e3] ; ; corners has 256x256 elements, one for each pixel. ; Each pixel has 4 corners and a latitude and a longitude. corners = make_array(256, 256, 2, 4, /float, value=!values.f_nan) cal_glats = reform(cal.thg_asf_fsim_glat[1, *, *, 0])') cal_glons = reform(cal.thg_asf_fsim_glon[1, *, *, 0])') cal_elevs = reform(cal.thg_asf_fsim_elev[*, *, 0])') for x=0, 255 do begin   for y=0, 255 do begin     ; ignore all data with elevations below 10 degree     if cal_elevs[x,y] lt 10. then $       continue     corner_lats = [ cal_glats[x, y], cal_glats[x+1, y], cal_glats[x+1, y+1], cal_glats[x, y+1] ]     corner_lons = [ cal_glons[x, y], cal_glons[x+1, y], cal_glons[x+1, y+1], cal_glons[x, y+1] ]     inds = where(finite(corner_lats), cc)     if cc lt 4 then $       continue     ; here we do the conversion to magnetic coordinates     ; and map coordinates     for c=0, 3 do begin       tmp = cnvcoord(corner_lats[c], corner_lons[c], 110.)       corner_lats[c] = tmp[0]       corner_lons[c] = tmp[1]       tmp = calc_stereo_coords(corner_lats[c], corner_lons[c])       corner_lats[c] = tmp[0]       corner_lons[c] = tmp[1]     endfor     corners[x,y,0,*] = corner_lats     corners[x,y,1,*] = corner_lons     endfor endfor ; Plot the map and overlay the ASI data map_plot_panel, xrange=[-30,5], yrange = [-30,5], coords='magn' ; Set color bar and levels rad_load_colortable, /whitered cin = FIX(FINDGEN(color_steps)/(color_steps-1.)*(ncolors-1))+bottom lvl = scale[0,s]+FINDGEN(color_steps)*(scale[1,s]-scale[0,s])/color_steps ; Load the first image brightn = reform(dat.thg_asf_fsim[*,*,0]) for x=0, 255 do begin   for y=0, 255 do begin     if ~finite(corners[s,x,y,0,0]) then $       continue     ; Determine color of pixel     color_ind = (MAX(where(lvl le ((brightn[x,y] > scale[0,s]) < scale[1,s]))) > 0)     col = cin[color_ind]     ; Plot the pixel     polyfill, reform(corners[s,x,y,0,*]), reform(corners[s,x,y,1,*]), color=col   endfor endfor ; Overlay the position of the ASI overlay_asi, 'FSIM', coords=coords, symsize=.6 Easy, eh?! And yes, the field-of-view of the ASI is skewed. That is due to the AACGM coordinate system. Of course, you can now overlay radar data using the rad_fit_overlay_scan command.

Radar File Issues
How do I...
... create my own fitACF/EX files? This method has become deprecated. Please refer to RST for data processing. You use one of the two binaries provided by DaViT, located in path-to-davit/bin. To create fitACF, use: path-to-davit/bin/make_fit -new file.rawacf > file.fitacf To create fitEX, use: path-to-davit/bin/make_fitex2 -new file.rawacf > file.fitex
... create my own grd files? This method has become deprecated. Please refer to RST for data processing. You use the binary provided by DaViT, located in path-to-davit/bin. To create a grd file from several fitEX files, use: path-to-davit/bin/make_grid -new -c -i integration_in_seconds file*.fitacf > file.grd
... create my own map files? This method has become deprecated. Please refer to RST for data processing. You can use the binary provided by DaViT They are located in path-to-davit/bin. To create a map file from a grd file, do the following steps: path-to-davit/bin/map_grd -new file.grd > file.map path-to-davit/bin/map_addhmb -new file.map > file.bnd.map path-to-davit/bin/map_addimf -new file.map > file.imf.map path-to-davit/bin/map_addmodel -new file.map > file.imf.map path-to-davit/bin/map_fit -new file.map > file.fit.map More comprehensive documentation can be found here.