Dies ist eine Sammlung meiner idl Entdeckungen This is a collection of my IDL discoveries

In the web you may find many IDL tutorials or helb pages …

IDL current reference @ NV5 Geospatial software :
IDL 8.9. Reference @ NV5
IDL 8.9. funclisting @ NV5
IDL 8.9. whats new @ NV5

IDL 8.9. references @ l3harris geospatial

IDL 8.8 referecnes (dead, Certifacte expired ??? ) It was once Exelis, but they deviate to harris:
refefrences … links to newest version at harris

A site which provides the IDL-help pages which you also cen get at the IDL prompt with IDL>? functionname
https://www.irya.unam.mx/computo/sites/manuales/IDL/idl_Left.html

At FZ-Jülich you can find a manual for IDL programming:
Praxisbezogene IDL Programmierung

You can get Help from IDL by just typing ? at the IDL prompt or for a certain idl function by typing ? <functioname>. This will open your web browser with a local help page.

We once decided to use a revision control system (RCS) for the idl routines under /home/hatpro/idl

IDL uses some Environment variables containing the pathes to its libraries etc. see here . Internally some of them can be accessed as system variables: The environment varibles usually have the form IDL_<NAME>. Inside IDL they are accesible via !<NAME> Basic path to IDL is in environamet variable IDLDIR inside IDL it is accesible via !DIR which is currently (idl version 9.0) /opt/nv5/idl90 .

There is an alternative for IDL based on the same syntax, function names etc.:
Gnu Data Language (GDL)

http://gnudatalanguage.sourceforge.net/
http://linux.die.net/man/1/gdl
weiterführende hilfe für die installation
list of functions missing in GDL
… auf der qeb ist das standard ubuntu paket gnudatalanguage version 0.9 installiert.
Call with
$>gdl
But it is since years in development and many functions are missing or not working properly.
See this list: /home/jschween/ideen/gdl/readme.txt

Starten eines Scirptes von der Shell geht in GDL 0.9 nicht wie in IDL mit \$> idl @idl_script sondern über eine bash pipe:
\$> cat @idl_script | gdl

Einige Besonderheiten: siehe diese einfachste Einleitung in GDL

• shell Befehle müssen nicht per spawn gestartet werden sondern können in einer Zeile die mit \$ beginnt stehen, z.B.:
  GDL> $ cp idl_script.pro gdl_script.pro
• …

start from the terminal prompt:

~>idl

Help with

IDL> ?

or

IDL> ? command-name

idls own routines are located in

/usr/local/itt/idl64/lib

print, 'hello world !'

  ; number of data points
  N = 100
  ; an array of x-values in interval 0..1
  x = findgen(N)/N
  ; an array of y-values
  y = sin(x*2*!Pi)
  ; plot the data
  plot, x, y

add axis titles and title to plot:

  plot, x,y , title='this is the title' , xtitle='x-axis', ytitle='ordinate'

manipulate the range

  plot, x,y , xrange = [ -0.1 , 1.5 ] , xstyle=1 , yrange = [ -1.2 , +1.5 ] , ystyle=1

use linestyle and colors

  plot, x,y , linestyle=1, color=3

plot a second function over the first

  oplot, x, 0.6*sin(2*x) , linestyle=2, color=255

  openw, file, 'data.txt',  /Get_LUN
  printf, file , 'x=', x 
  printf, file , 'y=', y
  close, file
  free_lun, file

  ; add two to all values
  x2 = x+2
  ; multiply all by three
  x3 = x*3
  ; square and squareroot
  x_square = x^2
  x_root = sqrt(x)
  ; mean etc.
  ym = mean(y)
  yvar = variance(y)
  ysigma = stddev(y)
 

<command> [ , <arguments> ] [ , <name>=<variable> ] [ , /<flagname> ]

• be aware of the comma !
• lowercase/uppercase does not matter
• <arguments> is a comma separated list like x,y,z which must come in the specified order
• the last arguments can be omitted, e.g. x,y instead of x,y,z (if <command> can handle this )
• <name>=<variable> are named variables inside command, these parameters can usually be omitted
• <name> arguments must not come in a special order
• /<flagname> is an abbreviation for flagname=1

cd , "<path>" [ ,current=<var> ]

path in single (') or double (“) quotation marks !

SPAWN [, "<Command>" [, Result] [, ErrResult] ]

command in single (') or double (”) quotation marks !

e.g. pwd = print working directory:

IDL>SPAWN , 'pwd'

In some cases you may get error message from your program. For example curl:

IDL> spawn, 'curl -V'
curl: /opt/nv5/idl90/bin/bin.linux.x86_64/libcurl.so.4: no version information available (required by curl)
curl 8.5.0 (x86_64-pc-linux-gnu) libcurl/8.4.0-DEV OpenSSL/3.0.11 zlib/1.2.13 zstd/1.5.5 nghttp2/1.58.0
Release-Date: 2023-12-06, security patched: 8.5.0-2ubuntu10.6
Protocols: dict file ftp ftps gopher gophers http https imap imaps mqtt pop3 pop3s rtsp smb smbs smtp smtps telnet tftp
Features: alt-svc AsynchDNS HSTS HTTP2 HTTPS-proxy IPv6 Largefile libz NTLM SSL threadsafe TLS-SRP UnixSockets zstd
WARNING: curl and libcurl versions do not match. Functionality may be affected.

The first line and the last lines tell you that libcurl.so.4 from the /opt/nv5/idl90/... path makes problems.

Another example is ghostscript (gs) - which does not even work:

IDL> SPAWN, 'gs --version'
gs: /opt/nv5/idl90/bin/bin.linux.x86_64/libstdc++.so.6: version `GLIBCXX_3.4.32' not found (required by /lib/x86_64-linux-gnu/libgs.so.10)

To get around you have to understand where the system seeks libraries. This is explained e.g. here: https://askubuntu.com/questions/1500913/find-default-library-paths-location-in-ubuntu, or in the man pages of ld.so. To ensure that IDL only uses libraries of a certain version they set the variable LD_LIBRARY_PATH to the idl library folders. It is not set in standard bash environment, but in the SPAWN environment as it is a subrpocess of IDL:

IDL> SPAWN, 'echo $LD_LIBRARY_PATH'
/opt/nv5/idl90/bin/bin.linux.x86_64:/opt/nv5/idl90/bin/bin.linux.x86_64/dm/lib:/opt/nv5/idl90/bin/bin.linux.x86_64/jre/lib/server:/opt/nv5/idl90/bin/bin.linux.x86_64/jre/lib/server/..:/opt/nv5/idl90/bin/bin.linux.x86_64/jre/lib/server/../native_threads

We can ask Linux where we can find the libraries needed:

$whereis libcurl.so.4
libcurl.so.4: /usr/lib/x86_64-linux-gnu/libcurl.so.4

and

$whereis libstdc++.so.6
libstdc++.so.6: /usr/lib/x86_64-linux-gnu/libstdc++.so.6

We have several possibilities:

• adapt the LD_LIBRARY_PATH from within IDL before invoking SPAWN. This was proposed to me : https://www.nv5geospatialsoftware.com/Support/Forums/aft/8650 in the IDL forum for the PAtH variable but maybe IDL will have problems with that:

SETENV, 'LD_LIBRARY_PATH=/usr/lib/x86_64-linux-gnu/:'+GETENV('LD_LIBRARY_PATH')

• Alternatively adapt LD_LIBRARY_PATH within the SPAWN environment when calling your program

SPAWN, 'LD_LIBRARY_PATH=/usr/lib:$LD_LIBRARY_PATH; your_command_here'

• or even simpler unset this path variable in SPAWN: it is not needed as we are not doing IDL specific things:

SPAWN, 'unset LD_LIBRARY_PATH; your_command_here'

Check with curl:

  SPAWN, 'unset LD_LIBRARY_PATH ; curl -V'
  curl 8.5.0 (x86_64-pc-linux-gnu) libcurl/8.5.0 OpenSSL/3.0.13 zlib/1.3 brotli/1.1.0 zstd/1.5.5 libidn2/2.3.7 libpsl/0.21.2 (+libidn2/2.3.7) libssh/0.10.6/openssl/zlib nghttp2/1.59.0 librtmp/2.3 OpenLDAP/2.6.7
  ...

Check with ghostscript

  SPAWN, 'unset LD_LIBRARY_PATH ; gs --version'
  10.02.1 

⇒ WORKS!

to just compile do

.COMPILE <filename_1> [,<filename_2> [...]]

to compile and run:

.RUN <filename_1> 

with <filename_1> without the extension “.pro”

The semicolon ';' is the comment sign: everything after it is ignored by the interpreter

a PROGRAM is a text file with IDL commands ending with an END statement (lower and uppercase is allowed, dont mix cases !)

print, 'Hello world !'
print, '3*4=', 3*4
x = 3
y = 2.5
z = x/y
print, 'x=', x, ' y=',y,' x/y=z=', z
print, 'its time to say goodbye'
END

decisions are made with IF THEN ELSE

IF x eq 4 THEN print, 'x equals 4' ELSE print, 'x is not 4'

this must be in one line but can be extended of several lines with the '$' character:

IF x eq 4 $
  THEN print, 'x equals 4' $
  ELSE print, 'x is not 4'

or with BEGIN END blocks

IF x eq 4 THEN BEGIN
  print, 'x equals 4' $
END ELSE BEGIN
  print, 'x is not 4'
END

for safety you may extend the END statements with blocked statement i.e. ENDIF or ENDELSE

you may nest i.e. build structures like IF … THEN … ELSE IF … THEN … or use the CASE statement.

Loops can be implented with FOR, WHILE or REPEAT as one-liners:

FOR i = 0,5 DO PRINT, i

i = 0  
WHILE i lt 6 DO i = i+1

i = 1
REPEAT i = i + 1 UNTIL i gt 5

for more complex loops you may use BEGIN … END

FOR i = 0,5 DO BEGIN
  PRINT, i
  END

i = 0
WHILE i lt 6 DO BEGIN
  i = i+1
  PRINT, i
  END

i = 1
REPEAT BEGIN
  i = i + 1 
  PRINT, i
END UNTIL i gt 5

to make coding a bit safer you may append the name of the loop command to the END statement, i.e.

ENDFOR 
ENDWHILE
ENDREP 

a PROCEDURE is a subroutine that may take parameters and change them. Its definition is started with the word PRO and ends with a END statement

PRO my_proc, text_str, x_value=x_value
  print, 'text_str=', text_str
  print, 'x=', x_value
  x_value = x_value + 2
  END ; pro

it is called with its parameter list separated by commas:

x = 2
my_proc, 'some text', x_value=x

In this case x will have the value x+2 after the call

if a procedure is stored in a file with the same name as the procedure IDL will compile this file …

A FUNCTION is a subroutine that may return a result

FUNCTION my_func, text_str, x_value=x_value
  print, 'text_str=', text_str
  print, 'x=', x_value
  x_new = x_value + 2
  RETURN, x_new
  print, 'this is never printed !'
  END ; function

return means that the result is returned but also that the execution returns from the sub subroutine, for that reason the last print command is ignored.

A function is called with its parameterlist in braces

x = 2
y = my_func( text_str, x_value=x )

When finding a new name which might be a function procedure or another idl program idl searches in the directories given by the idl sytem variable !PATH. To make own programs visible to idl the directory where these programs could be found must be included in the path:

!PATH = !PATH + ':/home/hatpro/idl/lib'

To ensure that the directory is not added again an again use the following code:

s = '/home/hatpro/idl/lib'
if strpos( !PATH , s ) lt 0 then !PATH = !PATH + ':' + s

Constants are defined by value, the type can be specified

<name> = <Value>[<Type Letter[s]>]

e.g.:

• i = 0b ; creates a 1byte=8bit integer constant of value zero (0..255)
• i = not 0b ; creates a 8bit integer with all bits set to 1 i.e. value 255decimal='ff'x='11111111'b
• i = 0S ; creates a signed 16bit integer constant of value zero (-32767..32767)
• i = 0US ; creates an unsigned 16bit integer constant of value zero (0..65535)
• z = 0L ; creates a signed 32bit integer constant (Long) of value zero (+/-2.2mrd)
• z = 0UL ; creates an unsigned 32bit integer constant (Long) of value zero (0…4.mrd)
• z = 0LL ; creates a signed 64bit integer constant (Long) of value zero (+/-9.2e18)
• z = 0ULL ; creates an unsigned 64bit integer constant (Long) of value zero (0…1.8e19)
• x = '0F'XU ; an unsigned integer (U) written as Hex (X) with (decimal) value 15
• x = '08'OU ; an unsigned integer (U) written as Octal (O) with (decimal) value 8
• y = 7.3E1 ; creates single precision 73.0
• z = 7.3D1 ; creates double precision 73.0
• novalue = !VALUES.F_NAN ; creates a single precision floating point NAN ('not a number')
• novalue = !VALUES.D_NAN ; creates a double precision floating point NAN ('not a number')
• infty = !VALUES.F_INFINITY ; creates a single precission infinity (positive, a minus in front makes it negative)
• infty = !VALUES.D_INFINITY ; creates a double precission infinity

…

see 'defining constants' and 'data types in IDL help.

For NAN and infinity see 'Special Floating-Point Values' or Constant System Variables

single precission is a 4 Bit floting point number (~range: -1e38..1E38 accuracy ~1e-7),
double precision is a 8Bit floating point number (~range: -1E308..1E+308 )
see wikipedia floating point.

can be generated with

x1d = fltarr( N ) ; type float
x2d = dblarr( Ni, Nj ) ; type double precision
...

The arrays will be filled with zeros, see e.g. help for fltarr or the info page about array creation

To fill the array with a certain value:

x2d[*] = -1.0

To fill it with the value of its running index:

x1d = findgen( N ) ; float type
x2d = sindgen( Ni, Nj ) ; string type
...

For details see here

You may also use makearray:

s1d = MAKE_ARRAY( Ni, /string, /index )
x2d = MAKE_ARRAY( Ni, Nj , /float, value=-1.0 )
...

Multidimensional arrays can be accessed via indices for all the different dimensions but also via less indices:

x3d = fltarr( Ni, Nj, Nk )
print, x3d[i1d] ; returns the i1d scalar element of the array
print, x3d[i,j] ; returns the i+j*Ni scalar element of the array
print, x3d[i,j,k] ; returns the i+j*Ni+k*Ni*Nj((correct ?)) scalar element of the array
...

For two dimensional arrays the relation between one and two dimensional indices is:

i1d = i + j*Ni
i = i1d mod Ni
j = i1d / Ni

This can also be done with the idl function ARRAY_INDICES

A structure is a dataype which joins variables of different type in one variable. It can be generated with (see http://www.harrisgeospatial.com/docs/Creating_and_Defining_St.html)

var = { name : 'bill' , int1:1 , a2:2.4 }

or with the CREATE_STRUCT function (see http://www.harrisgeospatial.com/docs/create_struct.html)

var = CREATE_STRUCT( 'name', 'bill' , 'int1' , 1 , 'a2' , 2.4 }

usage is

var = CREATE_STRUCT(  [Name=<string>, ] <name1>, value1, [<name2>, value2, ... ] ...    )

where Name registers this strucure to IDL and forces you to not deviate from the list of names and types … or

var = CREATE_STRUCT(  [Name=<string>, ] <namelist>, value1, [ value2, ... ] ...    )

with <namelist> a string array e.g. <namelist> = [ 'name', 'int1', 'a2' ]

CREATE_STRUCT allows to append structures

var = CREATE_STRUCT(  [Name=<string>, ] [struct_1, ...] <name1>, value1, [<name2>, value2, ... ]  [struct_2,...]    )

Where struct_1, and struct_2 are previosly defined structures.

reference is made with

print, var

or with the . operator:

print, var.name
print, var.int1
...

if you do not know the elements of a strucuture you may use help with the flag /strucutre

help, var, /structure

or the TAG_NAMES function

tn = TAG_NAMES(var)
print, tn

You can construct arrays of structures with the REPLICATE function and manipulate its content and find certain elements (see http://www.harrisgeospatial.com/docs/Arrays_of_Structures.html):

; create array of 20 structures var as e.g. defined above:
var_list = REPLICATE( var , 20 )
; change name of entry name of structure number three
var_list[3].name='ff'
; find this entry
print, where(var_list.name eq 'ff' )
   

With version 8.0 IDL introduced LIST and HASH.

• LIST is an array like datatype which may contain a mixture of different variable types and structures. Elements can be accessed via their index (see LIST help page at irya.unam.mx)
• HASH is a structure like dataype which elements can be accessed via strings containing the names and which is more flexible than the STRUCUTRE data type (see HASH help page at irya.unam.mx)

In the following value can be of any IDL datatype.

Create a LIST:

x_list = LIST( Value1, Value2, ... )

Add a single element

x_list.ADD, value
...

Create a HASH:

y_hash = HASH( [Key1, Value1, Key2, Value2, ... )

Add another HASH

y_hash = y_hash + HASH( [Key21, Value21, ... )
...

Both datatypes are objects with several methods to get information or manipulate its contents - se help page referenced above.

The FOREACH statement can be applied to both data types. Assuming that you can print all their elements you can do this:

FOREACH x in x_list do print, x
FOREACH y in y_hash do print, y
...

IDL comes with two matrix multiplication operators # and ##.

• lets assume we have A=fltarr(NAcol,NArow) and B=fltarr(NBcol,NBrow)
• operator # multiplicates rows of the left by columns of the right matrix
  ⇒ it must be NArow=NBcol and we get C = A # B = fltarr(NAcol,NBrow)
  it is C[i,j] = sum(k=0..NArow-1) A[i,k] * B[k,j]
  scheme:

                          bbb
                aaaaa     bbb     ccc
C  =  A # B  =  aaaaa  #  bbb  =  ccc
                aaaaa     bbb     ccc
                          bbb

• operator ## multiplicates columns of the left by rows of the right matrix
  ⇒ it must be NAcol=NBrow and we get C = A ## B = fltarr(NBcol,NArow)
  it is C[i,j] = sum(k=0..NAcol-1) A[k,i]*B[j,k]
  scheme:

                 aaa               cccc
                 aaa      bbbb     cccc
C  =  A ## B  =  aaa  ##  bbbb  =  cccc
                 aaa      bbbb     cccc
                 aaa               cccc 

• you may use this to calculate the scalar product of two vectors:
  s = transpose(x) # y or s = x ## transpose(y) with x = fltarr(N) and y = fltarr(N)

• IDL is somewhat sloppy with the relation between vectors and matrices, i.e. one and two dimensional arrays:
  • lets assume X = fltarr(N) and Y = fltarr(N)
  • the following gives the same N*N matrix:
    • C = X # Y
    • C = X # transpose(Y)
  • and these give the transpose of C (D = transpose(C))
    • D = X ## Y
    • D = X ## transpose(Y)

• you may generate Matrices from vectors, i.e. from one-dimensional arrays:
• with X=fltarr(Nx) and Y = fltarr(Ny) you get (again using the flopiness of IDL)
  • C = X # Y or C = X # transpose(Y) ⇒ C = fltarr(Nx,Ny)
  • C = X ## Y or C = transpose(X) ## Y ⇒ C = fltarr(Ny,Nx)
• in linear algebra a product between two vectors giving a matrix is called the 'outer product' or in german 'Dyade'

• you may use this to generate matrices with equal rows or columns:
  • to generate a Matrix with Ny equal rows x = fltarr(Nx) use
    XX = x # replicate(1,Ny)
    i.e. replicate(1,Ny) is interpreted as a 1*Ny matrix
  • to generate a Matrix with Nx equal columns y = fltarr(Ny) use
    YY = replicate(1,Nx) # y
    i.e. replicate(1,Nx) is interpreted as a Nx*1 matrix - odd, he ! - try to be not confused !

                         1     xxx
                         1     xxx
XX  =  X # E  =  xxx  #  1  =  xxx
                         1     xxx
                         1     xxx

                         y     yyy
                         y     yyy
YY  =  E # Y  =  111  #  y  =  yyy
                         y     yyy
                         y     yyy 

• One can use this to generate matrices with a variable increasing only in one dimension and arrays of function values:

; generate Nx*Ny Matrix xx with Nx values increasing in x-direction from xa to xb 
x = xa + (xb-xa) * findgen(Nx)/(Nx-1) 
xx = x # replicate(1.0,Ny)

; generate Nx*Ny Matrix yy with values increasing in y-direction from ya to yb
y = ya + (yb-ya) * findgen(Ny)/(Ny-1) 
yy = replicate(1.0,Nx) # y

; Nx*Ny matrix with distances from a point with coordinates (cx,cy)
r = sqrt( (xx-cx)^2 + (yy-cy)^2 )

; Nx*Ny matrix with function values on every point 
f_xy = f( xx, yy )

• As matrix operations are faster than for-loops this may speed up your code a lot …

if an unknown procedure name <fname> is found idl looks in

1. the current directory
2. in the directories specified in IDLs '!PATH'-variable

for a file <fname>.pro or <fname>.sav

letters of the filename should be all lowercase.

the !PATH variable ist constructed as

<dir1>:<dir2>:...

if idl encounters an expression like

y = f_phi2( a )

It assumes firstly an array variable f_phi2 of which the index a is desired.
If f_phi2 is an external function it should be anounced with the command

FORWARD_FUNCTION f_phi2

see also http://www.physics.nyu.edu/grierlab/idl_html_help/F33.html#wp872909

Blocks are generated with:

BEGIN
...
END

where END can expanded to

ENDIF   	

if it terminates an IF THEN … block

ENDWHILE

if it terminates a WHILE DO … block …

This naming of the END key word is helpful to find errors in the BEGIN…END structure

soll laut idl mit widgets gemacht werden - ist aber aufwändig !!

es gibt zwei einfache methoden mit der Tastatur zu Arbeiten:

• der Read-Befehl erwartet eine Eingabe die mit RETURN abgeschlossen werden muss
• der GET_KBRD erlaubt eine direkte Tastaturabfrage:

der Read-Befehl:

c = ' '
read,'continue with RETURN',c

mit GET_KBRD kann man Tastaturabfragen machen:

  Result = GET_KBRD([Wait] [, /ESCAPE | , /KEY_NAME])

Zum Beispiel:

key = GET_KBRD( 0 )

liefert in key die letzte gedrückte Taste oder '' (nullstring) wenn keine Taste gedrücktwurde …

key = GET_KBRD( 1 )

wartet auf tastendruck und liefert in key die letzte gedrückte Taste

key = GET_KBRD( 1 , /escape )

wartet auf tastendruck und liefert in key die letzte gedrückte Taste oder die escape sequenz die der Taste entspricht - also zb. string(byte([27,91,66])) für die 'pfeil nach unten' taste

Liste verschiedener escape sequenzen:

damit lässt sich ein einfache block aufbauen mit dem man ein Programm steuern kann:

• das Programm soll eine variable i hoch zählen
• anhalten wenn ' ' (space) gedrückt wird
• eins weiter (i+1) gehen wenn '>' oder 'cursor nach rechts' gedrückt wird
• eins zurück (i-1) gehen wenn '<' oder 'cursor nach links' gedrückt wird
• an den anfang (i=0) gehen wenn '0' oder 'pos 1' gedrückt wird
• an das Ende (i=N-1) gehen wenn 'E' oder 'Ende' gedrückt wird
• beenden wenn ESC gedrückt wird

i = 0
finished = 0
key = ' '
while not finished do begin
 
  ; ... do something ...
  print, i,' ... '
 
  ; check keyboard - return escape sequences 
  ; if the last key was a space - wait
  if key ne ' ' then key = GET_KBRD( 0 , /ESCAPE ) $
                else key = GET_KBRD( 1 , /ESCAPE )
 
  ; if we got a key
  if key ne '' then begin
    if key eq ' ' then begin
      ; ' '(space) => increment counter
      if i lt N-1 then i += 1
    endif else if (key eq '<') or (key eq string(byte([27, 91, 68]))) then begin
      ; '<' or cursor left => go back
      if i gt 0 then i -= 1
      key = ' '
    endif else if (key eq '>') or (key eq string(byte([27, 91, 67]))) then begin
      ; '>' or cursor right => go further
      if i lt N-1 then i += 1
      key = ' ' 
    endif else if (key eq '0') or (key eq string(byte([27, 91, 72]))) then begin
      ; '0' or pos 1 => go to start
      i = 0
      key = ' ' 
    endif else if (key eq 'e') or (key eq string(byte([27, 91, 70]))) then begin
      ; 'e' or Ende => go to end
      i = N-1
      key = ' ' 
    endif else if key eq string(byte([27,27])) then begin
      ; ESC-key => set stop flag
      finished = 1
    endif
  endif else begin ; no key pressed
    ; increment counter
    if i lt N-1 then i += 1
    ; set 'stop flag' if at end of data
    if i eq N-1 then key = ' '
    print,'End of data ...'
  endelse
 
endwhile 

Dieser code findet sich als funktion key_control.pro in hatpro/idl/lib.

mit der funktion cursor es ist möglich die Mausposition abzufragen:

  ; wait for mousebutton down event
  cursor, mx,my, /down

oder

  ; get mouse position imediately - point coordinates in data-coordinates
  cursor, mx,my, 0, /data

Information über die gedrückte Maustaste erhält man über die systemvariable !Mouse.button (linke maustaste == 1, rechte ==2, mitte == 4).

Eine Anwendung findet sich in mouse_clip.pro in hatpro/idl/lib/.

can be split by

s2 = STRMID( s , index , [length , ... ] )

position of substrings with

i = STRPOS( s , sub [ ... ] )

split string at character sep into substrings, return array of substrings (sx) and give number of substrings (N)

sx = STRSPLIT( s , sep , count=N, /extract )

can be used to separate columns - see text files below…

you can also use regular expressions - e.g. lets assume you have strings with variable names and unit in braces:

s = 'wind speed (m/s)'

then you may want to split at the opening and closing braces:

ss = strsplit( s, '(\(|\))', /extract, /regex )

will deliver

ss[0] = 'wind speed'
ss[1] = 'm/s'

join array of strings into one string with separating character sep

s = STRJOIN( s , sep )

filenames can be separated into path and name by

file_dirname( <path to file> )
file_basename( <path to file> )

Formatted output is generated with a syntax similar to FORTRAN:

PRINT, 'ABC..', 1, 2, 1.2, FORMAT='(A8,I," / ",I03," * ",F5.1)'

generates:

 ABC..      1  / 002 *   1.2

see widhts and positions:

 123456781234  / 123 * 12345

where the FORMAT='( ... )' specifier said:

• A8 : print as alpha-numerc-string of 8 characters width
• I : print as integer with standard width (8 digits, leading blanks)
• “ / ”: insert string ' / '
• I03: print as integer with 3 digits width and leading zeroes
• “ * ”: insert string ' * '
• F5.1: print as flaoting point with total width 5 and one digit after the decimal point

To print an array x=fltarr(N) with N ≥ 10 into a table with colums separated by tabs (ascii code 9) use:

TAB=STRING(9b)
PRINT, 'x=',x,FORMAT='(A,10("'+TAB+'",F8.2))'

Time as Julian date can be formatted as well. To get a date as DD.MM.YYYY hh:mm:ss use:

PRINT, t, FORMAT='(C(CDI02,".",CMoI02,".",CYI04," ",CHI02,":",CMI02,":",CSI02))'

Where format included:

• C(…) : print in calendar format
• CDI02 : calendar format for day as 2 digit integer with leading zero
• CMoI02 : calendar format for Month as 2 digit integer with leading zero
• CYI04 : calendar format for Year as 4 digit integer with leading zeroes
• CHI02, CMI02, CSI02 : calendar format for hour, minute, second as 2 digit integers with leading zeroes

It is possible to print seconds as float:

PRINT, t, FORMAT='(C(CHI02,":",CMI02,":",CSf05.2))'

General syntax of a Format specifier is

[<repeat>]<Format code>[+|-]<width>[.<precision>[E<exponent>]]

Where <repeat> <width> and <precision> are integer numbers.

• <repeat> says that the format code is repeatedly used.
• <width> is the width of the whole expression, <width>=0 means that the width is adapted to the number, with <width>=03 leading zeros are set to fill a three digit field.
• <precision> for floating point numbers gives the digits after the point, with integers the resulting number is filled to this width with leading Zeros.
• E<exponent> gives the digits of the exponent for scientific notation
• Format specifiers can be grouped with brackets which may have also a <repeat> field in front.

<Format code> is a usually a single letter:

There are further formatting codes → type '? format codes' at the idl prompt. Or look at using explicitly format, format codes fortran or format codes cprintf.

function mean calculates the average of a complete array.

xm = mean( x , /nan )

Flag /nan makes idl to ignore nans. If you also need variance and higher moments use instead moment which does the calculation faster:

m = moment( x, /nan )
x_mean = m[0]
x_var  = m[1]
s_skew = m[2]
s_kurt = m[3]

With version 8 idl can restrict moment and mean to (only) one dimension:

; generate some data
N = 1000
x = randomu(rs,N) ; this gives a fltarr(N)
y = (randomu(rs,N))^2; this gives a fltarr(N)
xy = [ [x], [y] ] ; this gives a fltarr(N,2)
m_x_y = mean(   xy, dimension = 1 ) ; this gives a fltarr(2) equal to [ mean(x), mean(y) ]
m_xy  = mean(   xy, dimension = 2 ) ; this gives a fltarr(N) equal to m_xy[i] = mean([ x[i], y[i] ])
mom = moment(   xy, dimension = 1 ) ; this gives a fltarr(2,4) equal to  [ [moment(x)], [moment(y)] ]

with variable dimension indicating the number of the dimension over which to average. The first dimension has number 1.

The same can be done with the max and min functions.

max_x_y = max( xy, dimension=1 ) ; this gives a fltarr(2) equal to [ max(x), max(y) ]
max_xy  = max( xy, dimension=2 ) ; this gives a fltarr(N) equal to max_xy[i] = max([x[i], y[i]]) 
...

IDL becomes slow if you use loops like for, while or repeat and conditional statments like if or case and even slower if you combine them.

This webpage provides some tips to speed your code: http://www.danidl.co.uk/idl.speed.tips.shtml

You can speed up your code if you use idl routines, matrix operations etc. Some examples below:

We want a running average over 2*M points:

this is slow

M = 11
for i = 0, N-1 do begin
  ia = max([i-M,0])
  ib = min([i+M,N-1])
  xm[i] = mean( x[ia:ib] , /nan )
endfor

use instead

xm = smooth( x, 2*M , /edge_truncate , /nan )

The /nan key word let idl ignore NaN¹⁾

or in two dimensions if you want to average in x direction only.

Mx = 2*M
My = 1
xm = smooth( x, [Mx,My] , /edge_truncate , /nan )

alternatively:

Mx = 2*M
My =  1
kern = make_array( Mx, My, /float, value = 1.0/(Mx*My)  ) ; equal to kern=fltarr(Mx,My) & kern[*]=1.0/(Mx*My)
xm = convol( x , kern , /edge_truncate , /nan )

Adapt Mx and My if you want more … This opens also possibilities to more sophisiticated smoothing techniques. Eg. with a gaussian smooother equivalent to the arithmetic mean over M elements:

; sigma derived from arithmetic average width M
sigma_gauss = M/3.
; width should be at least 3sigma ... we take here 4 sigma
N_gauss = floor(sigma_gauss*4)*2+1
; kernel is gaussian bell 
kern = exp( - 0.5 * ( (findgen(N_gauss)-0.5*(N_gauss-1)) / sigma_gauss )^2 ) 
; nomrlaize it
kern = kern/total(kern)
; convolute
xm = convol( x , kern ,  /edge_truncate , /nan )

a running standard deviation is also easily done

x_sdev = sqrt( smooth( (x-xm)^2 , 2*M , /edge_truncate , /nan ) )

Instead of for-loops use the variance or moment routines of idl.

var_x = variance(x , /nan )
; or 
moms_x = moment( x , /nan  )
mean_x = moms_x[0]
var_x = moms_x[1]

For the cross covariance use

xm = mean(x,/nan)
ym = mean(x,/nan)
covar_xy = mean( (x-xm)*(y-ym) , /nan )

If your data set is small you should multiply a factor N_data/(N_data-2) where N_data = total(finite(x-xm) AND finite(y-ym))

for the time laged cross covariance use:

N_data = n_elements(x)
xp = x - mean(x,/nan)
yp = y - mean(y,/nan)
N_lag = 10 ; must be smaller than N_data !
covar_xy_pos = fltarr(N_lag) ; for positive lags = indgen(N_lag)
covar_xy_neg = fltarr(N_lag) ; for negative lags =-indgen(N_lag)
for i_lag = 0, N_lag-1 do begin
  covar_xy_pos[i_lag] = mean( xp[i_lag:N_data-1      ] * yp[  0  :N_data-1-i_lag] , /nan )
  covar_xy_neg[i_lag] = mean( xp[  0  :N_data-1-i_lag] * yp[i_lag:N_data-1      ] , /nan )
  endfor ; i_lag

And similar for the auto covariance by replacing y by x.

If you want to evaluate it for N_lag in the order of N_data it is much more time efficient to use the identity²⁾:

cov(x,y) = FFT( FFT(x) * CONJ(FFT(y)), /INVERSE )

but then there must be *not any NaN* in your dataset !
Be aware that:

• covariances for negative timelags are at the end of the iFFT array.
  
• this identity assumes periodicity of the time series i.e. x[i] = x[i+N], … etc.

If you want to use data only if it is above a threshold this will be slow:

for i = 0, N-1 do begin
  if x[i] lt threshold then begin
    x[i] = !values.f_nan
  endif ; x[i]
endfor ; i

use instead:

i_ts = where( x lt threshold , N_ts )
if N_ts gt 0 then x[i_ts] = !values.f_nan

the if N_ts gt 0 is necessary because i_ts will be equal to -1 if no data is found which will result in an error prior to version IDL 8. But IDL 8 treats negative indices as relative indices - so you might not even notice that something wrong happened.

And similar if you do not want to juggle with NaNs create a new dataset:

i_good = where( x gt threshold, N_good )
if N_good gt 0 then x_good = x[i_good] $
               else x_good = [] ; this gives !NULL (a pointer pointing nowhere ..) which may result in strange results

have the basic structure:

PRO <name> [,<p1> [,<p2> [,<p3> [...]]]] [,<v1>=<v1> [,<v2>=<v2> [,<v3>=<v3> [...]]]]
...
command block 1
...
END
[...
command block 2
...
END  ]

• <name> name of the program, it must be identical with the filename, except for the extension.
  The code above must go into the file <name>.pro
• <p1>,<p2>,… names of parameters
• <v1>=<v1>,<v2>=… names of named parameters
• command block 1 will be executed if the program is called
• command block 2 is optional and will be executed if the program is compiled…

Example:

PRO test ,x ,y ,dings=dings

IF dings THEN PRINT,x,y
END

must be saved as test.pro
and be compiled with

.RUN test

it can then be called with

test,'aha',3.14,/dings
test,'aha',3.14,dings=0
test,'aha',3.14,dings=1
test,'aha',/dings
...

some of these examples will produce errors…

To make parameters optional you should test whether they are present when the programm is called. This can be done with:

IF NOT ARG_PRESENT(dings) THEN dings = 0

ARG_PRESENT checks whether the variable dings was passed to the program
Nearly the same can be achieved with KEYWORD_SET which checks whether a variable exists and whether its value is not zero

IF NOT KEYWORD_SET(dings) THEN dings = 0

We have only a limited number of licenses for IDL (60 but every user needs 6 at once thus 10effective, of these are 6 for idl8.2 and 4 for 6.4). To reduce the time that IDL is blocked you may either build a standalone executable or use batch jobs:

In principle it is possible to generate standalone executables. They run in IDL virtual machine which executes the code ³⁾. Additional the virtual machine always asks for a mouse click - also if your application has no graphical user interface. The advantage is that it does not block a license.

Details how to generate and execute standalones are describe under “The idl virtual machine”:
http://www.physics.nyu.edu/grierlab/idl_html_help/distributing6.html#wp1004661

And a more compact description at idl coyote:
https://www.idlcoyote.com/tips/idl_icon.html

Here is how to do it:

• assume you have a main program main.pro which you want to run. In IDL execute the following command:

    .COMPILE main  ; compile main program
    RESOLVE_ALL    ; resolve all dependencies ...
    ; RESOLVE_ALL, CLASS=<object_names> ; ... if you use objects provide their names as a strarr in <object_names> 
    SAVE, /Routines, Filename='main.sav' ; generate the ''.sav'' file which may run in the virtual machine

• in windows you may just click on main.sav to run it
• in linux call from the command prompt:

    $ idl -vm=main.sav

If you want to run the idl from a bash script you may use an IDL batch job.

An idl program which name starts with a '@' character is interpreted as a batch job. It can be started from the comandline and occupies an idl license only as long it is running. You may edit the code between calls without blocking IDL.

Here is how to do it:

• assume your idl stuff is done in idl_program.pro
• write a batch file '@call_idl_program' (you can choose any name, but it must start with a '@' character)

; @call_idl_program
; compile and call idl_program.pro
.COMPILE idl_program ; compile your stuff
idl_program          ; run it
EXIT                 ; leave idl

• call it from command prompt with

$ idl @call_idl_program

You cannot pass paramters from the command line to your program. You either have to

• provide paramters in your @call… file, e.g. generate it from your calling script
• or call an idl program from within the @call… batch job
• or provide parameters via an input file you read in idl_program.pro

a simple idl program to provide paramters could look like

a = 1
time = julday( 1,1,1970, 0,0,0 )
x = 3.14*a
end

save as set_my_params.pro and in @call_idl_program just include a line

; @call_idl_program
; compile and call idl_program.pro
.COMPILE idl_program ; compile your stuff
set_my_params        ; set parameters
idl_program          ; call your program
EXIT                 ; leave idl

have in principle the same structure as programs, except that they return a value:

FUNCTION <name> [,<p1> [,<p2> [,<p3> [...]]]] [,<v1>=<v1> [,<v2>=<v2> [,<v3>=<v3> [...]]]]
...
command block
...
RETURN, <variable>

END

Example:

FUNCTION sinx,x
RETURN sin(x)/x
END

open (text)files for input into idl:

s = '' ; this line predefines s as a string and the file is read textline by textline

OPENR, file_unit, filename, /Get_LUN
...
  ; read one textline
  READF, file_unit, s
...
CLOSE, file_unit
FREE_LUN, file_unit

if the file is compressed with gzip it can be directly read by setting the flag /compress in OPENR

open (text)files for output from idl:

OPENW, file_unit, filename, /Get_LUN
...
  PRINTF, file_unit, s
...
CLOSE, file_unit
FREE_LUN, file_unit

Read a text file containing a table with known format

OPENR, file_unit, filename, /Get_LUN
...
  READF, file_unit, var_int, var_float, var_sci_float, format='(I3, F7.2, E12.2)'
...
CLOSE, file_unit
FREE_LUN, file_unit

Read a text file containing a table with columns separated by TAB's

  TAB = string(9b)
  OPENR, file_unit, filename, /Get_LUN
  ...
    READF, file_unit, s
    col = STRSPLIT( s , TAB , count=N_col , /extract )
  ...
  CLOSE, file_unit
  FREE_LUN, file_unit
 

Count lines in a text file can be done with function file_lines, filename

The above are simple alternatives to the very complex idl routine READ_ASCII :

data = read_ascii( filename )

will read the data into a strucutre containing one array:

help, data
Structure <1485998>, 1 tags, length=<N_col>*<N_row>*<size_of_float>, data length=<N_col>*<N_row>*<size_of_float>, refs=1:
  FIELD1          FLOAT     Array[ <N_col>, <N_row> ]

it is possible to exclude comments, a header, or read only parts of the file by defining a TEMPLATE.

This can be done with the function ASCII_TEMPLATE which opens an interactive dialog.

template = ASCII_TEMPLATE(  data_file_name )

In principle this function can detect several features of the file (like column separator, data type per column, …) But you have to click several buttons.

The function returns a structure with the following fields:

Structure <179f0e8>, 10 tags, length=120, data length=117, refs=1:
 VERSION         FLOAT           1.00000
 DATASTART       LONG                 0
 DELIMITER       BYTE        <delimiter ascii code>
 MISSINGVALUE    FLOAT               NaN
 COMMENTSYMBOL   STRING    '<comment symbol>'
 FIELDCOUNT      LONG         <N_col>
 FIELDTYPES      LONG      Array[<N_col>]
 FIELDNAMES      STRING    Array[<N_col>]
 FIELDLOCATIONS  LONG      Array[<N_col>]
 FIELDGROUPS     LONG      Array[<N_col>]

Accordingly it should be possible to construct it in the code to avoid calling the dialog:

; delimiter between columns
; delimiter = ',' ; comma for csv files
delimiter = string(9b) ; tab for tsv files

; read header
; the line where to find the column header , we assume here the very first line
i_header = 0
header = '' ; predefine as a string => file can be read textline by textline
OPENR, file_unit, filename, /Get_LUN
; read first textlines until we have the column header
for i = 0, i_header do READF, file_unit, header
;   one may alternatively search for a certain key in the text line with
;   key = '# time'+delimeter
;   N_header = 0
;   WHILE (strpos(header,key) lt 0) and (~ eof(file_unit)) DO BEGIN
;     READF, file_unit, header
;     N_header += 1
;     ENDWHILE ; key in header
;   i_header = N_header-1
; check whether there remains anything in the file
if eof(file_unit) then message, 'no data in file !'
; close file
CLOSE, file_unit
FREE_LUN, file_unit

; split header, and count columns separated by delimiter.
header = strsplit( header, delimiter, /extract, /preserve_null, count = N_col )

; define type of variables in different columns
; type: 0=skip field, 1=byte, 2=int16, 3=long32, 4=single, 5=double, 7=string, ... see idl help for size
; i assume here that the first 3 columns are strings, i.e. cannot be converted into numeric values and the rest can be converted to single precision floating point
N_col_str = 3
fieldtypes = [ replicate(7,N_col_str), replicate(4,N_col-N_col_str)  ]

template = { $
  VERSION         : 1.0, $ ; must be 1.0 ! IDL 6.0 does not know anthing else ...
  DATASTART       : i_header+1 , $  ; index to line where to start reading. First line=0. We assume here that data starts directly after the header
  DELIMITER       : delimiter, $ ; one character - must be defined above e.g. as TAB=string(9b) or comma=',' ...
  MISSINGVALUE    : !values.f_nan, $ ; idl inserts this if it cannot read
  COMMENTSYMBOL   : '#', $ ; one character ...
  FIELDCOUNT      : N_col, $ ; number of columns, number of elements of following fields must fit
  FIELDTYPES      : fieldtypes, $ ; type: 1=byte, 2=int, ... 4=single, 5=double, 7=string, see idl help for size
  FIELDNAMES      : header, $ ; these names will be used for the variables, i.e. entries in header must follow the rules for variable names
  FIELDLOCATIONS  : replicate(  0, N_col ), $ ; indices to column begin, if no delimiter is given (?), or if columns are of constant width
  FIELDGROUPS     : indgen(N_col) $ ; every field an independent group, i.e. one element in the struct data we are going to read
  ; FIELDGROUPS     : replicate(0,N_col) $ ; every field the same group - only possible if the same type ! => data = { header[0] : array(N_col:N_lines) }
  ; FIELDGROUPS     : [ replicate(0,N_col_str) , replicate(1,N_col-N_col_str) ] ,  $ ; one group for the strings, one for the floats => data = { header[0] : strarr(N_col_str,N_lines), header[N_col_str] : fltarr(N_col-N_col_str,N_lines) }
  }

; read the data
data = read_ascii( filename, template=template )

; show info about the structure of your data
help, data, /struct

...

• FIELDCOUNT and number of elements of FIELDTYPES, -NAMES, -LOCATIONS and -GROUPS must be the same !
• if you group your data
  • you will get for every group an M x N array with M the number of group members and N the number of records (or lines)
  • you nevertheless need to provide FIELDTYPES, -NAMES, -LOCATIONS and -GROUPS with FIELDCOUNT elements !
  • the name of every group is then the name of the first FIELDTYPENAME going into it

with READU binary data can be read:

...
; define x as long variable
x = 0l
; read it from file
READU, file_unit, x
...

with WRITEU binary data can be written:

...
WRITEU, file_unit, x
...

check for existence

file_info_rec = FILE_INFO( filename )
if file_info_rec.exists then print,'exists' $
                        else print, 'does not exist'
                        

or in one line

if (FILE_INFO( filename )).exists then print,'exists' else print, 'does not exist'

get file size

file_info_rec = FILE_INFO( filename )
print, 'file size=', file_info_rec.size

set file pointer to a certain position

OPENR, file, filename, /Get_LUN
...
POINT_LUN, file , file_position

get position of file pointer (⇒ negative file unit number!) :

POINT_LUN, -file , file_position

count lines in a text file

N_lines = FILE_LINES( file_name ) 

We prefer to save data as NETCDF files (=Network Common Data Format)
NetCDF is a platform independent binary format which includes meta information about its content.
It can be easily read in many programming platforms like matlab, python, C, …
As a binary format reading is very fast and access to a subset of the content is easy.

• NETCDF was developed by Unidata@UCAR see also en.wikipedia.org …
• a description can be found at https://docs.unidata.ucar.edu/netcdf-c/current/index.html
• beside the values they should contain meta-information called attributes
  this can be units of the values, instrument location, model run number etc.
  unidata provides conventions about attributes
• the data is organised in arrays i.e. every data record in a netcdf file has at least one dimension
• IDL netcdf overview: https://www.l3harrisgeospatial.com/docs/NCDF_Overview.html
• NetCDF error codes https://docs.unidata.ucar.edu/netcdf/documentation/4.8.0/nc-error-codes.html

I wrote functions nc_read.pro and read_nc_group.pro to read netcdf files into a structure comparable to the Python function xarray.open_dataset.
You can find them in /home/jschween/exp/ideen/idl/netcdf.

For details of IDLs netcdf functions see the IDL help pages for a NCDF Overview.
Below example codes rely on the basic NETCDF functions IDL provided since long. Since version 8.4.1 IDL provides a simplified interface which also provides compact access to netcdf data.

; create netcdf output file
; use flag /clobber => overwrite exisiting file (the verb to clobber means to hit or punish someone, what strange... )
nc_file = NCDF_CREATE( nc_filename , /clobber )

; set global attributes
NCDF_ATTPUT, nc_file, /Global, 'Creator' , 'Jan Schween (jschween@uni-koeln.de)'
NCDF_ATTPUT, nc_file, /Global, 'creation_date' , STRING( SYSTIME(/utc,/jul), f='(C(CDI02,".",CMoI02,".",CYI04,"/",CHI02,":",CMI02))' )

; define time as dimension
; the corresponding variable shall have the same name as it is intended to be a coordinate variable
N_time_id = NCDF_DIMDEF( nc_file, 'time', N_time )

; we use NaN for missing values
novalue = !values.f_nan

; define time variable with dimension time, 
;   as time is a coordinate variable its name should be the same as the dimension
;   additionally as coordinate variable it should be in ascending order and have no gaps (missing values, filled values)
;   we use here the IDL - time format julian day which allows the determination of calendrial dates or can be inversely determined from calendar dates.
time_id = NCDF_VARDEF( nc_file, 'time', N_time_id, /DOUBLE )
; set attribute for the units of time
NCDF_ATTPUT, nc_file, time_id, 'units' , 'days since -4712-1-1T12:0:0' ; format according to cfconventions, see https://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/ch04s04.html , this should work with   "ncdump -iv time <nc_file>"   ... not allways but sometimes(?) ...  
NCDF_ATTPUT, nc_file, time_id, 'calendar', 'gregorian' ; this is the standard calendar - in principle you can omit this.
NCDF_ATTPUT, nc_file, time_id, 'long_name', 'time as julian day'
NCDF_ATTPUT, nc_file, time_id, 'missing_value', double(novalue)
NCDF_ATTPUT, nc_file, time_id, '_fillvalue', double(novalue) ; allthough unidata gives "_FillValue" with capital letters, you may get an error when using them... (see [[http://www.rhinocerus.net/forum/lang-idl-pvwave/575253-ncdf_attput-_fillvalue-problem-string-arrays.html|here]])
NCDF_ATTPUT, nc_file, time_id, 'comment', 'Julian day are the fractional days since 1.1.4713BC 12GMT. It can be converted to unix epoch (t_unix:seconds since 1.1.1970) with t_unix=(time-2440587.5)*86400.'

; alternativley time as unix epoch - wide spread time format usually of type integer, but then with overflow problem...
tunx_id = NCDF_VARDEF( nc_file, 'tunix', N_time_id, /LONG ) 
NCDF_ATTPUT, nc_file, tunx_id , 'units' , 'seconds since 1970/1/1 0:0:0'  ; format according to cfconventions, , evtly. '-' as separator in date? see https://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/ch04s04.html, no time zone info means automatically UTC

; create id for a variable dependent on time, 
; possible types are: /BYTE , /CHAR, /DOUBLE , /FLOAT , /LONG , /SHORT, /STRING, /UBYTE, /UINT64, /ULONG, /USHORT
; where /SHORT is a two byte integer, /LONG is a 4byte integer, and /UINT64 is an unsigned 8byte integer,
; and you may already guess: /U... means that the variable is an unsigned integer
var_id = NCDF_VARDEF( nc_file, 'some variable name', N_time_id, /FLOAT )
; multidimensional arrays like fltarr(N_time,N_space) are defined via
var_2d_id = NCDF_VARDEF( nc_file, '2D variable', [ N_time_id, N_space_id ], /FLOAT )
; single vars are defined by omitting the dimension
var_one_value_id = NCDF_VARDEF( nc_file, 'one_value', /FLOAT )

; set attributes
; these here are standard attributes (including the underscore) as defined in the unidata conventions ...
NCDF_ATTPUT, nc_file, var_id, 'long_name', 'long name'
NCDF_ATTPUT, nc_file, var_id, 'units', 'arb.units'
NCDF_ATTPUT, nc_file, var_id, 'missing_value', novalue
NCDF_ATTPUT, nc_file, var_id, '_fillvalue', novalue ; allthough unidata gives "_FillValue" with capital values, you may get an error when using them... (see [[http://www.rhinocerus.net/forum/lang-idl-pvwave/575253-ncdf_attput-_fillvalue-problem-string-arrays.html|here]])



; further variables and attributes ...
; ...

; STRINGS are stored as array of (unsigned) bytes called CHAR
; we therefore need the length of the string
; and the netcdf vaiable has to dimensions:
; lets assume we have str_var = strarr(N_str_var) 
max_str_len = max( strlen(str_var) )
; define dimension
str_len_id = NCDF_DIMDEF( nc_file, 'str_len', max_str_len )
str_var_id = NCDF_VARDEF( nc_file, 'str_var_name', [ str_len_id, N_str_var_id ], /CHAR )

; close define mode
NCDF_CONTROL, nc_file, /endef

; put variable into file ...
NCDF_VARPUT, nc_file, time_id, time
NCDF_VARPUT, nc_file, tunx_id, (time-2440587.5)*86400.

; STRING variable is as simple as the others ...
NCDF_VARPUT, nc_file, str_var_id, str_var
  
; close netcdf file
NCDF_CLOSE, nc_file

For time.units I follow here cfconventions.
But they are not very clear - they just give examples which differ from each ohter.
Accordingly time.units could have the form:

• '-' as separator in date
  • <base_unit> since <[-]YYYY-MM-DD hh:mm:ss>
• clock time with no leading zeros
  • <base_unit> since <[-]YYYY-MM-DD [h]h:[m]m:[s]s>
• '/' as separator in date
  • <base_unit> since <[-]YYYY/MM/DD hh:mm:ss>

Commandline tool ncdump cannot deal in every case with the clock-time in time.units.
what seem to work with ncdump is:

• you need dashes '-' as separators in date
• blank ' ' or 'T' as separator between date and time
• Leading zeros can be omitted
• that means
  • [-][Y[Y[Y]]]Y-[M]M-[D]D{ |T}[h]h:[m]m:ss

In the above example 'time' is the name of the dimension and of a variable. That is intended to make clear the close relation between N_time and time = dblarr(N_time) and part of the CF netcdf dataformat⁴⁾. But it is also possible to use different names, in this case ncview allows you to plot an array against other variables…

In netcdf 64 it is possible to have one 'UNLIMITED' dimension. I.E. this dimension does not need to be specified at creation time and data can be appended to the file. In idl this is done by using the keyword '/unlimited' (see also the example 'appending data …' below):

id_dim_time = NCDF_DIMDEF( nc_file, 'time', /unlimited )

There is a list of netcdf attribute conventions at the unidata pages

IDL help is not 100% clear about the meaning of the datatypes in netcdf :

• /BYTE 8bit signed/unsigned(?) integer … ncview assumes signed i.e. numbers between -128 and 127
• /CHAR this is also byte, but assumed to be characters i.e. a string can be stored as an array of char but you have to convert it in idl to string. Might be also interpreted as unsigned 8bit integers i.e. numbers between 0 and 255
• /SHORT 16bit integer (2 byte) probably signed, i.e. numbers between -32768 and +32767
• /LONG longword integers i.e. probably 32 bit signed integers i.e. numbers from − 2,147,483,648 to 2,147,483,647
• /FLOAT floating point numbers - probably IEEE single precision (32bit)
• /DOUBLE IEEE double precision (64bit)

; open netcdf file
nc_file = NCDF_OPEN( filename , /NOWRITE )

; get global attribute
NCDF_ATTGET, nc_file, /GLOBAL, att_name, glob_att

; get id of dimension
dim_id = ncdf_dimid( nc_file , 'dim_name' )

; get id of variable via its name
var_id = NCDF_VARID( nc_file, 'some variable name' )

; get attribute for variable
NCDF_ATTGET, nc_file, var_id, 'units', var_unit_str

; get dimension
NCDF_diminq, nc_file , dim_id, name_dim , N_dim 

; elements of a metadata group
; see on Harris
;   http://www.harrisgeospatial.com/docs/ncdf_ncidinq.html
;   http://www.harrisgeospatial.com/docs/ncdf_groupname.html#netCDF_2618656010_1006782
metadata_grp_id = NCDF_NCIDINQ( nc_file , Metadata_group_name )
mata_data_id = NCDF_VARID( metadata_grp_id, meta_data_name )

; get data
NCDF_VARGET, nc_file, var_id, variable

; get string data
str_var_id = NCDF_VARID( nc_file, 'str_var_name' )
NCDF_VARGET, nc_file, str_var_id, str_var
str_var = string( str_var )

; close netcdf file
NCDF_CLOSE, nc_file

Appending data is possible if the dimension is defined as unlimited. Defining the variables is as usual with one exception:

; define unlimited dimension time 
id_dim_time = NCDF_DIMDEF( nc_file, 'time', /unlimited ) ; only one inlimited dimension allowed (?)
id_dim_2D = NCDF_DIMDEF( nc_file, 'N_2D', N_2d )
; connect variables to it
id_time = NCDF_VARDEF( nc_file, 'time', id_dim_time , /DOUBLE)
id_data = NCDF_VARDEF( nc_file, 'data', id_dim_time , /FLOAT )
id_data_2D = NCDF_VARDEF( nc_file, 'data', [ id_dim_2D, id_dim_time ], /FLOAT ) ; for 2D data the unlimited dimension must be at the end

If the file remains open just append the data at the desired position

; append data assuming that variables //time// and //data// contain only new data
NCDF_VARPUT, nc_file, id_time, time, offset=[ N_time ]
NCDF_VARPUT, nc_file, id_data, data, offset=[ N_time ]
NCDF_VARPUT, nc_file, id_data_2D, data, offset=[ 0, N_time ]

If the file is reopened the number of elements already in the file must be determined:

nc_file = NCDF_OPEN( nc_file_name , /write )
; get id for variables time and data
id_time = NCDF_VARID( nc_file , 'time' )
id_data = NCDF_VARID( nc_file , 'data' )

; get number of elements in time in two steps
info_time = NCDF_VARINQ( nc_file , id_time )
NCDF_DIMINQ, nc_file, info_time.dim[0], name, N_time

; alternatively
id_time = NCDF_DIMID( nc_file , 'time' )
NCDF_DIMINQ, nc_file, id_time, name, N_time

; append data assuming that variables //time// and //data// contain only new data
NCDF_VARPUT, nc_file, id_time, time, offset=[ N_time ]
NCDF_VARPUT, nc_file, id_data, data, offset=[ N_time ]

With NETCDF 4 it is possible to compress data (see unidata.ucar.edu). IDL can read NETCDF4 file as of version 8.0 (see IDL ref pages for NCDF_CREATE). You can check the NETCDF version implemented in IDL with HELP, 'ncdf', /DLM.

With our IDL9 implementation this does not work anymore: you get an netcdf error -31 after NCDF_CREATE !?!?
Although it say it has version 4 implemented:

IDL> help, 'ncdf' , /DLM
** NCDF - IDL Network Common Data Format (NetCDF) support DLM (loaded)
  Version: netcdf-c 4.9.2, Build Date: NOV  8 2023, Source: NV5 Geospatial Solutions, Inc.
  Path: /opt/nv5/idl90/bin/bin.linux.x86_64/idl_netcdf.so

To use NETCDF compression you need to set the NETCDF4_FORMAT keyword when generating the file,
and the GZIP and SHUFFLE keywords when defining a variable, eventually also CHUNK_DIMENSIONS :

; create netcdf output file
nc_file = NCDF_CREATE( filename , /clobber, /NETCDF4_FORMAT )
...
g_zip_level = 5 ; a value between 0..9 for low...high compression but also fast...slow access
; chunk dimensions: can help the compression but can also be omitted
cd_time  = N_time / 4 ; i guess that is the way to use it
cd_range = N_range / 4 ; i guess that is the way to use it

; get id for a variable dependent on time, 
; possible types are: /BYTE , /CHAR, /DOUBLE , /FLOAT , /LONG , /SHORT
data_id = NCDF_VARDEF( nc_file, 'data', N_time_id, /FLOAT , GZIP=g_zip_level, CHUNK_DIMENSIONS=cd_1, /SHUFFLE )
; multidimensional arrays are defined via
data_2d_id = NCDF_VARDEF( nc_file, 'data_2d', [ N_time_id, N_range_id ], /FLOAT , GZIP=g_zip_level , CHUNK_DIMENSIONS=[cd_time, cd_range ], /SHUFFLE )
...

• GZIP determines the compression rate. Allowed values are 1..9 for lowest to highest compression or fastest to slowest compression (see help for NC_VARDEF#GZIP).
• CHUNK_DIMENSIONS determines how large blocks are compressed, it is by default the whole size of the variable except for variable of undefined length (see help for NC_VARDEF#CHUNK_DIMENSIONS). This means you can omit it. I dont know whether this paramter can improve compression.
• SHUFFLE tries to sort bytes before compression (see help for NC_VARDEF#SHUFFLE).
• with just GZIP=9 i achieve compression ratios around only Size_comp/Size_org = 85% for some smooth data, when appliying /SHUFFLE this improves to ~40%

Data in Netcdf files can be organized in groups.

• to get a list of the ID's all groups in a netcdf file use NCDF_GROUPSINQ
• to get the names of the groups use NCDF_groupname
• to get the ID of a group from its name use NCDF_NCIDINQ

Sample code:

; open file
file_Id = NCDF_OPEN(fileName)
; get list of goup-id's
group_id_list = NCDF_GROUPSINQ( file_id )
for i = 0, n_elements( group_id_list )-1 do begin
  group_id_i = group_id_list[i]
  group_name_i = NCDF_GROUPNAME( group_id_i )
  print,  i, '->', group_id_i, '->', group_name_i
  ; in principle a group can again contain groups ...
  ; note that we now use the groups_id to reference the parent group where we used before the file id.
  sub_group_id_list = NCDF_GROUPSINQ( group_id_i )
  ... 
endfor ; i

One idea of netcdf is to read a file directly and automatically into variables with names equal to the names in the file. That means that netcdf-variable names have to follow certain rules and not every character is allowed in a variable name in an netcdf file.

In general: a valid variable name in a programming language like fortran is an allowed variable name for a netcdf file. There is an [incomplete] description of the rules at the Unidata pages:

A valid identifier can be:

ID = ([a-zA-Z0-9_]|{UTF8})([^\x00-\x1F\x7F/]|{UTF8})*

which says:

• the first char must be alphanumeric, an underscore, or a UTF8 multibyte encoded char.
• subsequent characters may be anything except control characters [\x00-\x1F] , forward slash [/], and del [\x7F]. So the excluded ascii decimal values are 0-31, 47, and 127.
• A UTF-8 multibyte encoding always has the high order bit set …

to be complete: the following are not allowed as well:

• '(' and ')'
• '/'
• '%'
• white space
• ':' (colon)
• '“' (double quotes)
• '*' '+' … (operators)
• etc. …

The reason for this strict variable naming is that the variable names from inside the file should also be used as variable names in a program. Reading would be thus very simple: just read the file and you have the variables with the same name as in the file. Unfortunately this is not impelemented in idl…

PLOT, <x-array>, <y-array>

generates a rectangle (box) with automatically scaled abscissa (bottom) and ordinate (left side) and plots <y-array> versus <x-array>. There are many optional parameters and many addtional keywords to modify the plot.

The scaling is by default something like 10% wider than the range of the data but this can be modified by the style keyword. You set essentially single bits of it and may combine (at least some of) them:

[x|y|z]style
1 -> Force exact axis range.
2 -> Extend axis range.
4 -> Suppress entire axis
8 -> Suppress box style axis (i.e., draw axis on only one side of plot)
16 -> Inhibit setting the Y axis minimum value to 0 (Y axis only) 

suppress entire x-Axis (top and bottom of the box) with style option

xstyle=4

suppress top axis and right y-axis

xstyle=8, ystyle=8

Suppress annotation of axis-ticks

xtickname=replicate(' ',30)

Die Achseenbeschriftung besteht aus dreisssig Leerzeichen, IDL beschriftet die Achse nicht …
… oder multiplot.pro von mir verwenden.

The rectangle around a plot is in idl called AXIS - there is a command AXIS to just define this rectangle without plotting it. But there is no way to return to a formerly defined one.

Contour plots (isolines and/or color shading ) can be done with the command countour:

CONTOUR, z, x, y

I takes in principle the same options and keywords as the plot command plus some other.

Polar contours can be done with:

POLAR_CONTOUR, field, az_list, R_list

you can define the length of intervals between ticks and the number of sub intervals with {x|y}tickinterval and {x|y}minor where minor gives the number of subintervals - not ticks ! If you want to have ticks every 2.5 units and minor ticks every 0.1 then:

plot, ... , xtickinterval = 2.5, xminor=5

or

!X.xtickinterval = 2.5
!X.minor = 5

if you do not want minor ticks set

xminor = 1

you may want to annotate an axis with your own words. You then have to set {x|y}ticks with the number of intervals (not ticks !) {x|y}tickname with an array of strings for the annotations and {x|y}tickv for the locations where to put the tics. Note that {x|y}ticks is the number of intervals i.e. n_elements({x|y}tickname)-1. If you omit xticks IDL assumes 5 intervals or 6 ticks !.

Here is an example with a logarithmic time axis:

hours_per_day = 24d
days_per_year = 365d

; basic unit of our x-axis is years - thus it is:
hour = 1./(hours_per_day*days_per_year)
day =  24*hour
month = 30*day  ; ... well roughly - average length of a month is 30.416.. and 30.5days for a normal and a leap year, respectively, or ~30.44days on average

x_tikvals = [ hour,   day,  month,   1.,   10. ,  100. ,  1000.,  10e3,   100e3,   1e6   ]
x_tiknams = [ '1h',  '1d',  '1Mo',  '1y', '10y', '100y',  '1ky', '10ky', '100ky', '1My'  ]
N_ticks = n_elements(x_tiknams)

plot, [1], [1], $
  xstyle = 1, $ 
  xrange = [ 0.5 * hour, 2*1e6 ], $
  /xlog, $
  xticks = N_ticks-1, $  ; holy IDL - why do you need things like this ?
  xtickname = x_tiknams, $
  xtickv    = x_tikvals, $
  ystyle  = 1, $
  yrange = [ 0.5, 4.5 ], $
  yminor = 1, $
  /nodata

Similar for an x-axis based on julian day:

hours_per_day = 24d

; basic unit of our x-axis is (julian)days - thus it is:
hour = 1./hours_per_day
day =  24*hour
month = 30*day  ; ... roughly: average length of a month is 30.416.. and 30.5days for a normal and a leap year, respectively, or ~30.44days on average
year = 365.*day ; ... roughly ...
minute = hour/60.
second = minute/60.
milisec = second/1000.

x_tikvals = [ milisec, second, minute,  hour ,  day ,  month,    year   ]
x_tiknams = [ 'msec',   'sec',  'min', 'hour', 'day', 'month',  'year'  ]
N_ticks = n_elements(x_tiknams)

plot, [1], [1], $
  xstyle = 1, $ 
  xrange = [ 0.5*x_tikvals[0], 2*x_tikvals[N_ticks-1] ], $
  /xlog, $
  xticks = N_ticks-1, $  ; holy IDL - why do you need things like this ?
  xtickname = x_tiknams, $
  xtickv    = x_tikvals, $
  ystyle  = 1, $
  yrange = [ 0.5, 4.5 ], $
  yminor = 1, $
  /nodata

IDL provides in procedure PLOT the parameter XTICKUNITS which allows automatic annotation of a time axis. If t is time in julianday then

plot, t, x, XTICKUNITS = 'Months'

generates a plot where every first of a month a tick is set. Annotations are the three first letters of the name of month (format = 'CMOA' )

Other possible values for XTICKUNITS are:

• 'Days'
• 'Hours'
• 'Minutes'
• 'Seconds'

A tick is set in fixed intervals optimal for the length of the axis (N days, or N hours etc.). Annotation is the day of the month or hour of day etc. but can be formatted with XTICKFORMAT. For example:

plot, t, x, XTICKUNITS = 'Days', XTICKFORMAT = '(C(CDI02,",",CMOI02))'

Will give annotations of the form DD.MO instead of the month name.

Another example:

plot, t, x, XTICKUNITS = 'Hours', XTICKFORMAT = '(C(CHI02,":",CMI02))'

Will give annotations of the form HH:MM instead of just hour of the day as one or two digit number.

Another way is to explicitly define format and intervals, e.g. :

plot , t, y , XTICKFORMAT = '(C(CHI02,".",CMI02,"."))', XTICKINTERVAL = 3/24., XMINOR = 3

annotates HH.MM every 3 hours (=3/24. of a day) and generates two minor ticks for 3 subintervals i.e. minor ticks at full hours.

You may want to have annotations in two lines: use code !C for a newline

plot, t, y, XTICKFORMAT = '(C(CDI02,".",CMoI02,"!C",CYI04))'

{x|y|z}style = 
 1 => force exact range
 2 => extended range
 4 => supress entire axis (dont draw line and tic marks ...)
 8 => supress box style - axis only one one side of plot
16 => supress y=0 value for y axes - you can also use the /YNOZERO flag in the plot command...

linestyle = 
 0 => Solid
 1 => Dotted 
 2 => Dashed 
 3 => Dash Dot 
 4 => Dash Dot Dot 
 5 => Long Dashes 

PSYM =
 1 => Plus sign (+) 
 2 => Asterisk (*) 
 3 => Period (.) 
 4 => Diamond 
 5 => Triangle 'up'
 6 => Square 
 7 => X 
 8 => User-defined. See USERSYM procedure. 
 9 => Undefined 
10 => Histogram mode. Horizontal and vertical lines connect the plotted points, 
      as opposed to the normal method of connecting points with straight lines. 
      Vertical lines are located in the middle between subsequent x.
      See Histogram Mode for an example. 
      

Size of the symbol is controlled via the symsize keyword. A value of one is about one charactersize.

IDL does not know filled symbols, but one can use alternatively the userdefined symbol PSYM=8 and the USERSYM procedure. To generate a filled triangle pointing to the top use:

  N_sym = 3
  w0_sym = 0.
  w_sym = findgen(N_sym+1)/N_sym*2*!Pi
  usersym, sin(w_sym+w0_sym) , cos(w_sym+w0_sym), /fill
  plot, x , y, psym=8

To make different symbols set:

More complex symbols like stars might be generated by using something like

usersym, sin(w_sym*N_arms)*sin(w_sym+w0_sym) , sin(w_sym*N_arms)*cos(w_sym+w0_sym), /fill

Or you use programs from the internet like plotsym.pro from http://idlastro.gsfc.nasa.gov/

To plot multiple data rows use the OPLOT command:

  ; plot y versus x
  PLOT, x, y
  ; plot y2 versus x on same axis
  OPLOT, x, y2

If you want to use several different axes use the AXIS command. To plot some data y and y2 on the same x-axis but on the different y-axes on the (standard) left and the right side use the following code:

  ; define test data
  N = 1000
  x = findgen(N)/N
  y = sin(2*!pi*x)
  y2= 100.*exp(-((x-0.5)/0.1)^2)
 
 
  ; plot y data on (standard) left y-axis ...
  ; ... leave additional space of 10 charwidths for x-axis-labels on both sides
  PLOT, x, y1, ytitle='y1', XMARGIN=[10,10]
  ; define second y-axis on the right side , and save it for use
  AXIS, yaxis=1, yrange=[min(y2),max(y2)], ytitle='y2', /save
  ; plot data y2 onto this axis
  OPLOT, x, y2

The keyword yaxis=0 refers to a y-axis on the left side of the plot with tics poiting to the right side (i.e. inward). The keyword yaxis=1 refers to an y-axis on the right side of the plot with tics poiting to the left side (i.e. inward). Nevertheless the position of the axis can be set by giving the x-position of the axis. Here is an example for a third axis:

  ; additional data...
  y3 = atan(10.*(x-0.3))/(!pi/4.)
  ; position for third y-axis - left to the first axis
  x_icpt = -0.13
  ; define third y-axis placed at x=x_icpt, and save it for use
  AXIS, x_icpt, yaxis=0, yrange=[min(y3),max(y3)], ytitle='y3', /save
  ; plot data y3 onto this axis
  OPLOT, x, y3

Of course you can define additonal x-Axes you then have to use the parameter XAXIS=1 (Xaxis on top) or XAXIS=0 (x-axis at the bottom).

Sometimes the title/label for the second axis is not readable (especially if you are creating .ps-plots). An easy solution is to type: !p.position=[0.1,0.1,0.9,0.9] in front of the plot-commands. (Stefan) … More elaborated is to use the XMARGIN keyword in the plot command (added in the example above (Jan)) - or by setting it globally with:

  !X.MARGIN = [10,3]

X margin gives the width of the margin in characters on the left (XMARGIN[0]) and right (XMARGIN[1]) side of the plot axis. Similarily the bottom and top margins can be set with

  !Y.MARGIN = [4,2]

or with the YMARGIN keyword in the plot command.

To get multiple sub plots arranged in a matrix use !P.multi :

!P.multi = [ $
  N_remaining , $ ; number of plots remaining on page , a zero starts a new page
  N_column , $ ; number of columns
  N_row , $ ; number of rows
  N_z_stack , $ ; number of plots stacked in z dimension
  top_2_bottom $ ; Flag: first fill columns before proceeding to next, default is to first fill rows ('column major')
   ] 

Usage:

!P.multi = [ 0 , 2, 3 ]

plots come in a new (erased) 2×3 array, next plot will go into top left position, following one in top row second column, etc.

!P.multi = [ 0 , 5, 5, 1, 1 ]

plots come in a 5×5 array, next plot will go into top left position, next one in second row, first column, etc.

!P.multi = [ 0, 3, 3 ]
...
!P.multi = [ 5 ]

Next plot goes into row 2 colomn 2

XYOUTS, X, Y, String, alignment=alignment, charsize=charsize, width=width, clip=clip, noclip=noclip, ...

print string at pos x,y

• x,y in data coordinates (axes) - but that can be changed with flags /normalized or /device
  • omitting x,y position will be last (string?) position. in a new plot this is at 0,0 normalized coordinates (i.e. lower left corner)
• alignment is a floating point number:
  • 0.0 : left bounded, i.e. lower left corner of text is at x,y or text starts at x,y
  • 0.5 : centered, i.e. center of text at baseline is at x,y
  • 1.0 : right bounded, i.e. lower right corner is at x,y
  • yon can also give numbers below 0 or above one - if you know what to do with that
• charsize is size of character in multiples of character base size (independent of !P.charsize)
• charthick is linethicknes of drawn fonts (Hershey…)
• width returns the width of the plotted text in normalized coordinate units, if no other parameter is given the text is not printed.
• clip=[x0,y0,x1,y1] defines a rectangle where to clip the output text. You need to set noclip=0 to take effect
• noclip a flag that specifies whether the current clipping region (→ parameter clip or !P.clip) is used, default is noclip=1 i.e. you have to give noclip=0 to take clipping effect.
• more details xyouts@harris

• it might be useful to give positions in units of character height or width they can be determined from the system in device coordinates:

; char height and (average) width in device coordinates 
ch = float(!D.y_ch_size) 
cw = float(!D.x_ch_size) 
; ... in normalized coordinates  
ch = float(!D.y_ch_size) / !D.y_size
cw = float(!D.x_ch_size) / !D.x_size
; ... in data coordinates (only in linear coordiante system !)
ch = float(!D.y_ch_size) / (!D.y_size * !Y.S[1])
cw = float(!D.y_ch_size) / (!D.y_size * !X.S[1])

• it seems as if for pixel based devices it is !D.y_ch_size = 12 pixel. To achieve a certain charsize q_ch relative to bitmap size do:

; charsize in normalized units
; ch = !P.charsize * float(!D.y_ch_size) / !D.y_size
; I want ch = q_ch% <=> q_ch = charsize * float(!D.y_ch_size) / !D.y_size * 100
; => charsize = q/100. * !D.y_size / float(!D.y_ch_size)
q_ch = 3. ; percent of bitmap height
!P.charsize = q_ch/100. * !D.y_size / float(!D.y_ch_size)

Um ein einzelnes Symbol an bestimmter Stelle in einem Plot zu positionieren:

PLOTS,  xpos, ypos , psym=<symbol id>, ... 

there are three font systems in IDL:

• !P.Font = -1 Hershey-Vector Fonts (scalable and rotatable Vectorfonts, default of idl)
• !P.Font = 0 'device Fonts' - in the Postscript device these are Postscript Fonts - they are better readable than true-type fonts
• !P.Font = 1 True Type Fonts (TTF if installed). With TTF some of the below described formatting stuff is not possible - I (jan) find them rather 'pixelig'

IDL provids the possibility to change the font within the xyouts routine or the [x,y,z]title of plots.

• put an escape sequence of the form !<command> in the string to be printed
• different fonts can be selected with !<digit> here comes a selection
  • !3 is the Hershey Helvetica Font which is the default. IDL calls it “simplex Roman” - obviously the people at IDL mixed up the names
  • !6 is the Hershey Font “complex Roman” where every lines is made 3times to make the letters more bolder.
  • !4 is the Hershey Font “simplex Greek” greek letters single line: greek letters in principle in standard encoding, i.e. a ist alpha etc. but q instead of r is rho … there is somewhere a shift in the encoding table! (see "Using vector fonts" ).
  • !7 is Hershey Font “complex Greek” greek letters double lines
  • !9 oder !M is “special Math” - Integral and other special Math-signs this font is used for only one letter
  • !X returns to the 'Entry font'
  • for examples of the available fonts see IDL help under Idl help "Using vector fonts". There you will also see the encoding. You can use the indices from this tables directly. To e.g. print the degree sign which you find in the default font 3 under octal '260' use:
    xyouts, x,y, '15'+string("260b)+'C'
    (the notation is correct: quotation marks do not close and b stands for ocal !)⁵⁾

Examples with simplex greek

• a small greek phi
  !4f!X (take font greek simplex (!4) set an f return to default font (!X))
• a large greek phi
  !4F!X

You may also use the complex greek (2lines per stroke to make letter more bold)

• a greek psi
  !7w!X

One may define constants for the different letters:

c_greek_font = '!4' ; simplex greek
;c_greek_font = '!7' ; complex greek
c_alpha   = c_greek_font+'a!X'
c_beta    = c_greek_font+'b!X'
c_gamma   = c_greek_font+'c!X'
c_delta   = c_greek_font+'d!X'

c_epsilon = c_greek_font+'e!X'
c_zeta    = c_greek_font+'f!X'
c_eta     = c_greek_font+'g!X'
c_theta   = c_greek_font+'h!X'

c_iota    = c_greek_font+'i!X'
c_kappa   = c_greek_font+'j!X'
c_lambda  = c_greek_font+'k!X'
c_mu      = c_greek_font+'l!X'

c_nu      = c_greek_font+'m!X'
c_xi      = c_greek_font+'n!X'
c_omikron = c_greek_font+'o!X'
c_pi      = c_greek_font+'p!X'

c_rho     = c_greek_font+'q!X'
c_sigma   = c_greek_font+'r!X'
c_tau     = c_greek_font+'s!X'
c_upsilon = c_greek_font+'t!X'

c_phi     = c_greek_font+'u!X'
c_xi      = c_greek_font+'v!X'
c_psi     = c_greek_font+'w!X'
c_omega   = c_greek_font+'x!X'

c_infinity= c_greek_font+'y!X'

; usage
xyouts, 0,0, c_phi+'('+ c_theta+') = '+ c_delta

Hersehy vector fonts are one of the oldest computer fonts. They were developed in the late sixties by A.V. Hershey, see wikipedia, Paul Bourke or Idl help "Using vector fonts". They were designed for a plotter i.e. a device which draw lines with a pen. The letters consist solely of polygons or straight lines. To achieve bold(er) letters two or three lines are drawn close to each other. The respective fonts are called then duplex or triplex respectively. For examples see https://www.nv5geospatialsoftware.com/docs/Using_Hershey_Vector_Fon.html. To find out the byte code of the symbols shown in the tables add the (decimal) values at the left and top of the respcetive rows and columns. To get e.g. an $\alpha$ , i.e. the symbol in the row labelled with 96 and the column labelled with 1 you can use string(byte(96+1)).

setting the post script device up with

SET_PLOT,'ps' ; use postscript device
!P.font = 0 ; use postscript device fonts 
DEVICE, /ENCAPSULATED, /portrait, /COLOR, /HELVETICA , filename='test_font.eps'

; set some font indices (overwrite the original ones )
; index must be in [3..19]
; odd indices are normal , even are italic/oblique
; we need one without serifs (helvetica) one with (times) a typerwriter (courier) and symbol 
DEVICE, /HELVETICA ,           FONT_INDEX=3
DEVICE, /HELVETICA , /OBLIQUE, FONT_INDEX=4
DEVICE, /TIMES ,               FONT_INDEX=5
DEVICE, /TIMES ,     /ITALIC,  FONT_INDEX=6
DEVICE, /COURIER ,             FONT_INDEX=7
DEVICE, /COURIER ,   /OBLIQUE, FONT_INDEX=8
DEVICE, /SYMBOL ,              FONT_INDEX=9 

; char height in normalized coordinates
ch = float(!D.y_ch_size) / !D.y_size

; some letters...  
abc_rst = 'abcdefghijklmnopqrst...'
; character size factor
cs = 2.5
; sample output: switch to font #4 with !4 and back to 'standard' with !X
; gor through fonts 3-9 defined above
for i = 3, 9 do begin
  si = string(i,f='(I0)')
  s = '!'+si+'font #'+si+':'+abc_rst+'!X'
  print, s
  xyouts, 0.01, 1-(i-3+1)*cs*ch, s, charsize=cs, /normal
endfor

device, /close

see the end of the Device fonts section.

You may generate super- and sub-scripts and further primitve type setting command based on Grandle and Nystrom (1980)⁶⁾. This is also used in other software packages see e.g. the documentation of NAG e.g. here or here — sorry all the links before are dead. Here are some still alive NAG docs pages: NagText, NAGText@gatech.edu, NAGGraph@gatech.edu.

The following commands are available (see also IDL help Embedded Formatting Commands)⁷⁾:

• supersuperscript !U ('Up' faktor 0.62)
• subscript !D ('Down' faktor 0.62)
• subsubscript !L ('Low' faktor 0.62)
• index !I (faktor 0.44)
• exponent !E (faktor 0.44)
• back to Normal !N
• numerator position !A ('above' division line)
• denominator position !B ('below' division line)
• new line !C ('carriage return')
• special characters !Z(b1[,b2[,…]]) where b1,b2,… are the 16 bit hexadecimal unicode codes of the desired character[s]
• Save position !S
• Restore Position !R
• exclamation mark !!

Most of these commands only affect the vertical position except for !R which restrores also the horizontal position.

Examples:

• super and subscripts
  • $A^b$ (A to the power of b)
    A!Eb!N \\(A !Exponent b !Normal )
  • To get both super and sub script directly after your letter you need to save the horizontqal position:
    $A^b_i$ (A superscript b index i)
    A!S!Eb!N!R!Ii!N
    (A !Save !Exponent b !Normal !Restore !Index !Normal)
    If super and subscript have different lengths you first need to set the shorter one.

• X-axis Annotation in two lines: use code !C for a newline

  plot, time, y, $
    xtitle = '!CDate (!S!Add.mm!R!B YYYY !N)', $
    xtickformat='(C(CDI02,".",CMoI02,"!C",CYI04))', $
    ...

• the xtitle starts with a newline (!C) because IDL does not recognize that the annotations have two line now.
  !S!A saves the position after the ”(“ and writes “dd.mm” as nominator,
  !R!A returns to the position after ”(“ and writes the denoniator
  !N ends the nominator and returns to normal before the closing ”)“ is set.
  This generates something like:

  ------+-----------+----------> 
      24.12       24.01
      2019        2020

• something similar can be achieved by using xtickformat='(C("!S!A",CDI02,".",CMoI02,"!R!B",CYI04,"!N"))' ... with no centering of the YYYY content

• special characters:
  • the degree sign:
    • a small letter 'o' in the exponent: $15^oC$ : 15!Eo!NC
    • degrees as the special character with octal code 260: string("260b)
    • via its decimal code 176: string(byte(176))
    • or via its hexcode 0xB0: !Z(B0)

• special characters: greek letters in simplex greek and complex greek
  • $\beta$: !4b!X (switch to simplex greek (!4), draw b, go back to previous font (!X))
  • $\beta$: !7b!X (switch to complex greek (!7), draw b, go back to previous font (!X))
    note that encoding of the greek letters follows only for the four first letters the latin alfabet i.e. a hives $\alpha$ b give $\beta$, but r gives not $\rho$ you need q instead, and s not $\sigma$, you need t instead.

• cdot (center dot for multiplication)
  • $\cdot$: '!9'+string(byte('056'O))+'!X'
    select special and math font (!9), select the cdot symbol (string(byte(…))), return to previous font (!X)

• integral sign with integration boundaries
  • $\int_a^b dx$ : '!9'+string(byte(105))+'!X!S!Da!N!R!Ub!N dx'
    select special math font (!9), integral sign (string(byte(105))), restore previous font (!X), save position (!S), lower boundary (!Da!N), restore position, upper boundary (!Ub!N), dx.
    the inteegral sign used here is rahter small alternatively one can use string(byte(73)) also in the MAth and special font (!9) but it is rather large …

• a vector and its absolutum
  • $\vec{v}$: '!Sv!R!9!U'+string(byte(54))+'!N!X'
    (save position(!S), draw letter v, go back to start position (!R), select special font (!9), upper position (!U), arrow (string(byte(54))), normal position (!N), previos font (!X) )
  • $||\vec{v}||$: '!9'+string(byte(35))+'!X'!Sv!R!9!U'+string(byte(54))+'!N'+string(byte(35))+'!X'
    where the aboslutum symbol is generated with '!9'+string(byte(35))

• more complex constructions show the limits:
  • fraction $\frac{a}{1-x}$
    !S!A a !R!B1-x!N
    i.e.: save position (!S), draw nominator (!A), go back to saved position and start denominator (!R!B), return to normal when done (!N at end)
    as there is no automatic centering the a in the nominator is surrounded here with two spaces which is not exact. And there is no fraction bar: you have to do it on your own which means you need to know the length of the object. This could look like this (with no guarantee):

  xyouts, x,y , '!S!A a !R!B1-x!N' , width=wdt 
  plots , [x, x+wdt ] , [ y+dy , y+dy ]%%

• The position of the fraction line must be corrected a bit upward therefore the dy⁸⁾. You may try to calculate dy but at a certain point it becomes incredible complex …

For more examples have a look at Formatting Command Examples.

With a call to plot (and similar functions) a grafic window is opened (except you have changed the output device, see below). This x-terminal (linux/linux) does not perform a refresh when covered by another window and it uses only 256 colors.

For an automatic refresh set the retain flag to two:

SET_PLOT,'X'
DEVICE, RETAIN=2

This works only if the x-window is opened for the first time. If a x-terminal is already open close it before you call the device.

If you want “true colors” (i.e. 24bit colors) use:

SET_PLOT,'X'
DEVICE, /DECOMPOSED

Colors are then defined by hexadecimal numbers. Red, green, blue and combinations can be defined by:

red   = '0000ff'x
green = '00ff00'x
blue  = 'ff0000'x
yellow= '00ffff'x
magenta='ffff00'x
cyan  = 'ff00ff'x
black = '000000'x
white = 'ffffff'x
grey  = '888888'x

Farbe wird in idl üblicherweise als Index in eine Farbpallete (=colortable) interpretiert. Diese Tabellen können maximal 256 Einträge enthalten. Einige plotdevices können aber auch true color (also 3*8Byte=24bit farbwerte).

Alle Ausgaben in das Plotdevice können farbig erfolgen - dies geschieht mit der color = <wert> option, wobei <wert> entweder ein index in eine Farbpalette ist oder ein true colorwert.

From version 8 on IDL takes truecolor rgb values which you give as 3byte integer values. You can write them as Hex-numbers like color='0000ff'x, color='00ff00'x and color='ff0000'x for red, green and blue, respectively.

Before version 8 IDL used colortables, the color appearing on the screen was defined in a table or palette, i.e. a list with 256 entires for the r,g,b truecolor values. There are several predefined colortables (currently 73). You can load them with

LOADCT, ct_index

where ct_index is the index of the color table. A command PLOT,x,y,color=4 plots data with the color defined in the palette for index 4. The default palette ist a grayscale grayscale with 256 entries.

To print a list of the color table names use

  loadct, get_names=s
  for i = 0, n_elements(s)-1 do print, i, ':', s[i]

Current list is

     0:B-W LINEAR
     1:BLUE/WHITE
     2:GRN-RED-BLU-WHT
     3:RED TEMPERATURE
     4:BLUE/GREEN/RED/YELLOW
     5:STD GAMMA-II
     6:PRISM
     7:RED-PURPLE
     8:GREEN/WHITE LINEAR
     9:GRN/WHT EXPONENTIAL
    10:GREEN-PINK
    11:BLUE-RED
    12:16 LEVEL
    13:RAINBOW
    14:STEPS
    15:STERN SPECIAL
    16:Haze
    17:Blue - Pastel - Red
    18:Pastels
    19:Hue Sat Lightness 1
    20:Hue Sat Lightness 2
    21:Hue Sat Value 1
    22:Hue Sat Value 2
    23:Purple-Red + Stripes
    24:Beach
    25:Mac Style
    26:Eos A
    27:Eos B
    28:Hardcandy
    29:Nature
    30:Ocean
    31:Peppermint
    32:Plasma
    33:Blue-Red
    34:Rainbow
    35:Blue Waves
    36:Volcano
    37:Waves
    38:Rainbow18
    39:Rainbow + white
    40:Rainbow + black
    41:CB-Accent
    42:CB-Dark2
    43:CB-Paired
    44:CB-Pastel1
    45:CB-Pastel2
    46:CB-Set1
    47:CB-Set2
    48:CB-Set3
    49:CB-Blues
    50:CB-BuGn
    51:CB-BuPu
    52:CB-GnBu
    53:CB-Greens
    54:CB-Greys
    55:CB-Oranges
    56:CB-OrRd
    57:CB-PuBu
    58:CB-PuBuGn
    59:CB-PuRd
    60:CB-Purples
    61:CB-RdPu
    62:CB-Reds
    63:CB-YlGn
    64:CB-YlGnBu
    65:CB-YlOrBr
    66:CB-BrBG
    67:CB-PiYG
    68:CB-PRGn
    69:CB-PuOr
    70:CB-RdBu
    71:CB-RdGy
    72:CB-RdYlBu
    73:CB-RdYlGn
    74:CB-Spectral

Die aktuell geladene Farbpalette kann mit

TVLCT, rgb , /Get

in das array rgb geladen werden, oder mit

TVLCT, red, grn, blu , /Get

in die arrays red, grn, und blu.

Der folgende Beispielcode liest die aktuelle Farbpalette und definiert ein paar Farben um und macht diese neue Palette zur aktuellen:

; get actual colortable
TVLCT, rgb , /Get
; size of this colortable
n_col = n_elements( rgb[*,0] )
; define some indices into the table
col_black = 0
col_red = 1
col_green = 2
col_blue = 3
col_yellow = 4
col_gray = 5
col_white = n_col-1
; define the rgb components of the colors related to this indices
rgb[ col_red   , * ] = [   0 ,   0 , 255 ]
rgb[ col_green , * ] = [   0 , 136 ,   0 ] 
rgb[ col_blue  , * ] = [ 255 ,   0 ,   0 ]
rgb[ col_yellow, * ] = [   0 , 255 , 255 ]
rgb[ col_gray  , * ] = [ 136 , 136 , 136 ]
; IDL uses the first and the last entry of the colortable 
; as background and primary drawing color i.e. as black an white 
; redefine it here ...
rgb[ col_black , * ] = [   0 ,   0 ,   0 ] 
rgb[ col_white , * ] = [ 255 , 255 , 255 ]
; make this new colortable active
TVLCT, rgb
...
; plot red axis with red line
plot, x,y , color=col_red
...

Ähnlich wirkt der IDL befehl TEK_COLOR der allerdings keine Namen definiert. First nine colors are: 0=black, 1=white, 2=red, 3=green, 4=blue, 5=cyan, 6=magenta, 8=orange.

Klassiche Farbpalleten die man auf diese Weise definieren könnte sind z.B. die CGA-Palette, die EGA-palette oder die VGA-pallette.

Einen Farbverlauf kann man ähnlich definieren. Der folgende code definiert ein Spektrum in dem, orientiert an vereinfachten Empfindlichkeitskurven der Farbrezeptoren des Auges die Kurven für rot, grün und blau gesetzt werden:

; calc. a spectrum starting with red going via yellow and green to blue and violet
; get current color table
TVLCT, RED, GRN, BLU , /Get
; size of this colortable
n_colors = n_elements( RED )
; maximum value of one channel
max_rgb = 255
; the very first and last points (col=0 and col=N-1) are excluded because idl uses them as background and foreground colors
for col = 1 , N_colors-2 do begin
  RED(col) = fix((cos(2.*!pi*0.8*(float(col-1)/(N_colors-3)    ))+1.0)/2.0*max_rgb)
  GRN(col) = fix((cos(2.*!pi*    (float(col-1)/(N_colors-3)-0.3))+1.0)/2.0*max_rgb)
  BLU(col) = fix((cos(2.*!pi*0.5*(float(col-1)/(N_colors-3)-1.0))+1.0)/2.0*max_rgb)
  endfor
; load this color table
TVLCT, RED, GRN, BLU

In the 'true color' or 'direct color' system are the 3 lowest bytes of an integer number interpreted as the intensity values for red, green and blue. If you write this as a hex number you may define rather easily arbitrary colors. Here are some basic colors as true colors (compare also with 'html colors' eg in http://www.w3schools.com/html/html_colors.asp):

col_black  = '000000'x
col_red    = '0000ff'x
col_green  = '008800'x
col_blue   = 'ff0000'x
col_yellow = '00ffff'x
col_gray   = '888888'x
col_white  = 'ffffff'x

To use true colors the plotdevice must be configured accordingly. This is done with the /decomposed keyword, e.g. for UNIX X-windows device:

SET_PLOT, 'X'
DEVICE, /decomposed

One may also define colors directly as 24bit numbers from red, green and blue values. Assume that red, grn and blu are byte values (8-bit numbers) in the range 0..255, the respective '24bit true color' values is then:

true_color_rgb = long(red) + ishft(long(grn),8) + ishft(long(blu),16)

or if you have your red, grn blu values as floats in the range 0..1

true_color_rgb = long(red*255) + ishft(long(grn*255),8) + ishft(long(blu*255),16)

!P is a structure containing system parameters describing the plot

!P.MULTI → multiplot

• !P.MULTI[0] number of plots remaining on the page
• !P.MULTI[1] plot columns per page.
• !P.MULTI[2] rows of plots per page.
• !P.MULTI[3] plots stacked in the Z dimension.
• !P.MULTI[4] is 0 to make plots from left to right (column major), and top to bottom,
  1 to make plots from top to bottom, left to right (row major).

!P.charsize → relative size of all characters, default: 1=normal size

• 2 would double the size of all text in the plot
• 0.75 would reduce the size of all characters by about 25%

!D is a structure containing parameters about the current plot-Device
!D.x_size and !D.y_size give the width of the device in device coordinates
!D.x_ch_size and !D.y_ch_size give the base size of characters in device coordinates
⇒ !D.y_ch_size / !D.y_size gives the base size of characters in normalized coordinates

!X, !Y and !Z are structures containing parameters determining the appearance of the X,Y, and Z axes
!X.Margin = [left, right] give the width of the left and right margin between the axes and the plot region in characterwidths
!Y.Margin = [bottom, top] give the width of the left and right margin between the axes and the plot region in characterheights

Some variables contain information about the last plot

!X.S and !Y.S give the scaling factor for conversion from data to normalized coordinates. It is

x_norm = !X.S[0] + !X.S[1] * x_data
x_norm = !X.S[0] + !X.S[1] * log_10(x_data)

for linear and logarithmic scaling respectivley

!X.crange = fltarr(2) contains the range used by the x-axis.

!X.window = fltarr(2) contains the normalized x-coordinates of start and endpoint of the axis.

!X.region = fltarr(2) contains the normalized x-coordinates of the plotting region

Every single plot goes into a region which should contain the plot itself, the axis labels, the axis titles and the plot title. !X.Margin determines the distance of the x-axis to the bottom margin of this plot region.

If you want to know where ticks are made then use keywords xtick_get and ytick_get when calling plot:

plot, x, y, xtick_get=xticks, ytick_get=yticks 

See also the Graphics System variables in the IDL help. A (horrible) description of the different coordinate systems of IDL can be found under Direct Graphics Coordinate Systems .

For the pixel based Windows, X and Z Device it is possible to define a bit-operation for every pixel drawn. I.e the color bits for every pixel of a line are first connected with the background of its position. This can be used to make lines contrasting with the underlying background. The operation can be defined with the SET_GRAPHICS_FUNCTION parameter of the DEVICE function :

DEVICE, GET_GRAPHICS_FUNCTION = old_gfct, SET_GRAPHICS_FUNCTION = 6

with graphics_function = 6 = GXxor an XOR operation between pixel color and the content of the pixel. The GET_GRAPHICS_FUNCTION variable stores the current state. The default is GXcopy = 3 - i.e A list of the possible grafics functions can be found under idl help → Keywords Accepted by the IDL Devices -> SET_GRAPHICS_FUNCTION

The result of this pixel operations can give weird results …

IDL does not know transparent overlays but on pixel based devices you can implement it with TVRD and TVRD. Lets assume you want to put a bitmap BMP with size Nx, Ny at a position x0 y0 (pixels) with a tranparency factor of q_tra:

; get current content of canvas as a bytarr( N_channel, Nx, Nx )
bgnd = TVRD( x0,y0, Nx, Ny, true=1 )
; blurr it a bit
N_blurr = 8
bgnd = SMOOTH( bgnd, [ 1, N_blurr, N_blurr ] , /edge_truncate )
; create a new bitmap by combining bmp and bgnd 
bmp = byte( q_tra * float(bmp)  +  (1-q_tra) * float(bgnd) )
; put the new bmp into the canvas
TV, bmp, x0, y0, true=1, /device

This is implemented in mk_legend.pro

in classical 2D plots IDL knows three coordinate systems. In plot commands as plot, oplot, xyouts, plots you can refer to them via key words

• /data the coordinate system defined by the current axes
• /normal normalized coordinates within the axes running from zero to one
• /device device coordinates in pixel or postscript points

to convert between them idl uses several parameters (see Direct Graphics Coordinate Systems ):

• /data to /normal is determined via !X.S :
  x_norm = !X.S[0] + !X.S[1]*x_data for !X.type=0 (linear) and
  x_norm = !X.S[0] + !X.S[1]*log(x_data) for !X.type=1 (logarithmic)
  and the respective forms for y (→ !Y.S) and z (→ !Z.S)
• /normal to device is determined via !D.X_VSIZE and !D.Y_VSIZE:
  x_device = x_norm*!D.X_VSIZE

If you want scalable high quality plots with great fonts use encapsulated postscript (eps). Postscript is a vectorgrafics format and encapsulated means that the grafic can be easily included in other documents (if your software can deal with eps ).

But be aware that if your plot is made up by many plot commands the result will be an incredible large file. Nevertheless you can convert eps to pdf (which is still a vector grafics format but has an internal compression, see below) or to png (see below). If an eps files is very large the pdf will be large (about 1/4) and a png will be 'normal' (what you may used to). Never use jpeg for line plots or everyting which is not a photograph.

to produce an encapsulated postscript plot:

SET_PLOT,'ps' ; select postscript device
!P.font = 0 ; select postscript fonts
DEVICE, /ENCAPSULATED, /portrait, /HELVETICA, /COLOR ; basic setup for colored EPS 
DEVICE, /COLOR, /decomposed ; in idl 8 is true color = /decompose possible 
; DEVICE, FILENAME="idl_plot.eps" ; put a filename here if you do not ike the default from idl ('idl.ps') but note that this name is reset with every DEVICE,/CLOSE call
; DEVICE, XSIZE=17.8, YSIZE=12.7 ; size in centimeter ! these are the default values corresponding to a 640*480 pixel image printed at 96pixel per inch...
DEVICE, XSIZE=29.4, YSIZE=21.0 ; size in centimeter ! DIN A4
DEVICE, /SYMBOL ,  FONT_INDEX=4 ; make greek letters available as font #4

; device keyword output is intended to add own postscript code,
; BUT it generates also ps-code for a new page.
; See also comment in http://www.physics.wisc.edu/~craigm/idl/down/pxperfect.pro
; DEVICE, OUTPUT = '%%Title: plot by Jan Schween'

; in post script we need thicker axis lines
!X.thick = 3
!Y.thick = 3

; plot in black on white background
!P.color = '000000'x
!P.background = 'ffffff'x

; for a good readable plot charsize should be larger
!P.charsize = 2.

...
plot commands
...

DEVICE,/CLOSE

; if we did not provide a plotname with DEVICE, plotname+'.eps' we can rename it here...
SPAWN, 'mv idl.ps '+plotname+'.eps' ; rename outputfile

• mit /portrait wird das koordinatensystem nicht gedreht und die bbox stimmt…
  
• mit /landscape wird das koordinatensystem gedreht und die bbox stimmt nicht mehr …
  

in idl 6.4 ist unter postscript farbe nur aus einer farpalette mit 256 enträgen möglich. Mit:

SET_PLOT,'ps'
DEVICE, ... /COLOR, BITS_PER_PIXEL=8 ...

ist es möglich 2 dimensionale arrays intarr(Nx,Ny) als farbige bitmaps in PS darzustellen. Die Farbe ergibt sich dabei aus der aktuellen colortable - die hat maximal 256 einträge also einen 8Bit langen indexraum.

Im Postscript device kann man zwar true color jpegs darstellen:

READ_JPEG, filename, img
TV, img, TRUE=1

Allerdings wird das jpg bild anscheinend auf eine 8Bit = 256Farb-Palette reduziert…

In idl 8.2 kann man mit true color arbeiten:

SET_PLOT,'ps'
DEVICE, /COLOR, /decomposed, ...

ermöglicht die verwendung von true color - z.B.:

plot, x,y, color='00ffff'x 

benutzt Gelb (=max red + max green) als Farbe.

Device coordinates in the postscript device are in 1/1000cm (1000pixel per cm). To make a plus sign at 1cm,1.5cm use

plots, 1*1000,1.5*1000, psym=1, /device

Postscript - non encapsulated :

SET_PLOT,'ps'
DEVICE [  ,/portrait | ,/landscape ]
...
plot commands
...
DEVICE,/CLOSE

(→ …/idl/eps_test/ps_landscape.pro)

• IDL produces per default a 'seascape' plot erzeugt (in gv muss man seascape oder landscape+swap landscape anwählen)
• /landscape thus has no meaning !
• /portrait turns the coo-system by 90° to the left …
• in contrast to /ENCAPSULATED the bbox is ok then !
• for details of the postscript device see here

Basic commands, select device driver

SET_PLOT, Device

e.q. for Postscript

SET_PLOT,'ps'

or for colored EPS

SET_PLOT,'ps',/ENCAPSULATED,/COLOR

with only the SET_PLOT command output will be send to the file 'idl.ps'
this can be changed by

DEVICE, FILENAME='<filename>', ...

the plot-file is then finalized with

DEVICE, /CLOSE

or

DEVICE, /CLOSE_FILE

• postscript is a vector grafics format the predecessor of pdf.

That means grafs can be enlarged or shrinked without loosing quality.

• there is no big difference between PS and EPS:
  • EPS means that the graph can be imported into other applications.
  • That means only one single page is allowed - no multiple pages
  • the graf must have an information about its size - the bounding box (bbox)

Convert postscript to png with imagemagick:

convert -density 120 <ps-filename> <filename>.png

density 120 means 120 dpi, increasing this values increases quality and size of the png. You must have ghostscript installed. For more control you may want to directly use the ghostscript call for this:

gs -dSAFER -dBATCH -dNOPAUSE -sDEVICE=png16m -r120 -dTextAlphaBits=4 -dGraphicsAlphaBits=1 -dEPSCrop -sOutputFile=<filename>.png <ps-filename>

Parameters -dSAFER -dBATCH -dNOPAUSE makes gs to run this as batch job with no questions to the user,
-sDEVICE=png16m defines true color rbg png as output format,
you may use png256 or png16 for 256 or 16 element palette images which give substantilly smaller images but you may experience ditherung artefacts (for details on the gs devices see gostscrip webpages).
-r120 means 120dpi and gives a 1388×992 pixel images for a DINA4 landscape postscript,
-dTextAlphaBits=4 ensures smooth characters after rastering,
-dGraphicsAlphaBits=1 avoids gridlines after rastering (http://www.idlcoyote.com/ps_tips/psstripes.html),
-dEPSCrop crops the image at the bounding box.

One can convert eps files into pdf using the program ps2pdf of the ghostscript package. The following command creates a pdf file quelle.pdf from the file quelle.eps where the 'papersize' is adapted to the BoundingBox of the eps file.

ps2pdf -dEPSCrop quelle.eps

Since ps2pdf is based on ghostscript all gs options can be used.

To fit the plot to a standard paper size (i.e. Letter or DINA4 depending on your system settings) use

ps2pdf -dEPSFitPage quelle.eps

One can also invoke ghost script (gs) to control more precisely what you want.

To create a PDF with a rgb colors use :

gs -dBATCH -dNOPAUSE -sColorConversionStrategy=RGB -sDEVICE=pdfwrite -dEPSCrop -sOutputFile=quelle.eps.pdf  quelle.eps

To create a PDF based on PDF-A standard (the definitions file PDFA_def.ps should be in the gs/lib dirictory … no guarantee that this works …)

gs -dPDFA=1 -dBATCH -dNOPAUSE -sColorConversionStrategy=RGB -sDEVICE=pdfwrite -sOutputFile=quelle_a.pdf  PDFA_def.ps  quelle.eps

Scalable Vector graphics can be displayed by web browsers and might get the standard vector graphic format for web applications

There are converters from eps to SVG:

is a simple converter from pdf to svg. Simple because it converts fonts to vector grafics… The conversion is done in two steps:
1: convert eps to pdf

ps2pdf -dEPSCrop wind_rose_all.eps

2: convert pdf to svg

pdf2svg wind_rose_all.pdf wind_rose_all.svg

but fonts look really ugly !

is a package to convert postscript to something. To convert to svg use:

pstoedit -f plot-svg -psarg "-dEPSCrop" <file.eps>  <file.svg>

unfortunately the -dEPScrop option does not work properly. The graph sits in the lower corner and is partially cropped.

Raster-Oriented (or Bitmap) graphics can be created with:

  ; set the z-device => plot commands go into the Z-Buffer - this is just a bitmap
  SET_PLOT,'Z'
  ; configure this bitmap: true color (24bit=3byte) and size 640*480 (default)
  DEVICE, SET_PIXEL_DEPTH=24, /DECOMPOSED, set_resolution=[ 640 , 480 ]
  ...
  plot commands 
  ...
  ; read actual graf into array img as BYTARR( N_color , width , height )
  img = TVRD( true=1 )
  ; write into png image
  WRITE_PNG, 'idl.png', img
  ; or write into jpg image
  WRITE_JPEG, 'idl.jpg', img

This creates true color graphic files.

If a 256 color palette shall be used - which reduces the size of the grafics file and which is the default for idl and for gif but also possible with png the code must be changed:

  ; configure the bitmap: 256 color palette (8bit=1byte) and size 640*480 (default)
  DEVICE, SET_PIXEL_DEPTH=8, set_resolution=[ 640 , 480 ]
  ...
  ; read actual graf into array img as BYTARR( N_color , width , height )
  img = TVRD( )
  ; read actual palette
  TVLCT, red, grn, blu , /Get
  ; write into png image
  WRITE_PNG, 'idl.png', img, red, grn, blue

It is also possible to copy an image from the X-Device (i.e. the screen) but this might be disturbed by redraw commands or overlap from other windows (see /home/jschween/exp/ideen/idl/mk_png/…) ⇒ Use the Z-Buffer instead.

It is possible to create animated gifs by repeatedly writing to the same gif file with the /MULTIPLE flag set. GIF knows only palette based color with in maximum 256 colors. We thus either need to work in a 8bit color device, or we quantizie the image colors before providng it to the WRITE_GIF :

From a 8 bit color device :

  ; plot into virtual device = bitmap
  SET_PLOT, 'Z'
  ; configure this bitmap: 8 bit colors and size is 640*480 (default)
  DEVICE, SET_PIXEL_DEPTH=8, set_resolution=[ 640 , 480 ]
 
  ; get current color palette
  TVLCT, rgb , /get
 
  ; go through images
  FOR i_img = 0, N_img-1 DO BEGIN
 
    ; create plot 
    ...
 
    ; read actual graf into array img
    bmp = TVRD( )
    ; append img to gif-file
    WRITE_GIF, 'idl.gif', bmp,  rgb[*,0],rgb[*,1],rgb[*,2],  /MULTIPLE, delay_time=6, repeat_count=0
 
  ENDFOR ; i_img
 
  ; close gif
  WRITE_GIF, 'idl.gif', img, /CLOSE

From a true color device

  ; plot into virtual device
  SET_PLOT, 'Z'
  ; configure this bitmap: 24 bit colors and size is 640*480 (default)
  DEVICE, SET_PIXEL_DEPTH=24, set_resolution=[ 640 , 480 ]
 
  ; go through images
  FOR i_img = 0, N_img-1 DO BEGIN
 
    ; create plot 
    ...
 
    ; get current bitmap
    bmp = TVRD( true=1 )
    ; quantize colors - derive gif color translation table from first image
    if i_img eq 0 then gif = COLOR_QUAN( bmp, 1, red, grn, blu, get_translation = gif_col_trans ) $
                  else gif = COLOR_QUAN( bmp, 1, red, grn, blu,     translation = gif_col_trans )
    ; append img to gif-file
    WRITE_GIF, 'idl.gif', gif,  red, grn, blu,  /MULTIPLE,  delay_time=6, repeat_count=0
 
  ENDFOR ; i_img
 
  ; close gif
  WRITE_GIF, 'idl.gif', /CLOSE

• The /MULTIPLE keyword in the call of WRITE_GIF creates at the first call a GIF-file but does not close it. During subsequent calls the images are appended. Accordingly the GIF-file must be closed by a final call with the /CLOSE keyword.
• The delay_time is given in 1/100sec; delay_time=6/100sec here corresponds to a frame rate of ~17images/sec which is about the frame rate of the old film standard super8. Television has usually 25frames/sec corresponding to a delay of 4/100sec.
• The parameter repeat_count=0 makes the animation to be repeated forever.
• Graphics commands must go into a bitmap device in this example the Z-device (it might be possible to do this also in the X-device (UNIX/LINUX), or the win-device (windows). Be aware that the color setup is different than for printer devices (like e.g. postscript): background is black, foreground color (axes of a plot etc.) are white.

image files can be read into arrays, displayed, analyzed and manipulated

  ; read image
  img = read_png( 'img.png' )
  
  ; print out dimensions
  help, img
  
  ; determine dimensions
  img_size = size(img)
  ; number of color channels, columns and rows
  Nc = img_size[1]
  Nx = img_size[2]
  Ny = img_size[3]
  
  ; alternative: query
  qr = query_png( 'img.png' , info )
  Nx = info.dimensions[0]
  Ny = info.dimensions[1]
  
  ; black out lower left corner
  img[ 0 , 0:10, 0:10 ] = 0
  img[ 1 , 0:10, 0:10 ] = 0
  img[ 2 , 0:10, 0:10 ] = 0

  ; display
  x = 10
  y = 10 
  tv, img, x, y , true=1
  
  ; mark center
  plots, Nx/2, Ny/2 , psym=1, /device

The color scheme for bitmap devices as screen (SET_PLOT,'X' or SET_PLOT,'WIN') or the 'Z'-device (SET_PLOT,'Z') is different than for printer devices like e.g. postscript:

• In bitmaps background is black an drawing is done in white
• for printers background is white an drawing in black

This can be changed with:

  SET_PLOT, 'Z'
  ; set background color to white (last color in pallete)
  !P.BACKGROUND = 255
  ; set foreground color to black (first color in pallete)
  !P.COLOR = 0

If you are working with 'true colors'

  SET_PLOT, 'Z'
  ; set background color to white
  !P.BACKGROUND = 'ffffff'x
  ; set foreground color to black
  !P.COLOR = '000000'x

• JPG is good for images with continous color gradients like photographs:
  • (you remember: jp[e]g stands for joint photographic experts group).
  • Sharp edges (e.g. lines) in the image create echo like artefacts in jpg.
• PNG is good for line graphics and images with sharp edges:
  • images with large areas with the same color have good compression rates
  • scaling may lead to moiree efects along edges
  • it is possible to give every pixel in the image a certain transparency
• EPS / PS ([Encapsulated] Post Script) is good for pure line graphics
  • PS is a vector graphis format. I.E. the file contains plot commands like moveto, lineto, textout, …
  • if a plot contains a lot of drawing commands the file can become very large
  • bitmap images can be incorporated but are not compressed ⇒ files become large…
  • Graphics can be arbitrarily scaled. Scaling is mainly limited by the resolution of the output device.
• PDF (portable document format) is an extension of PS
  • it has additional compression included
  • EPS and PS images can be converted to PDF …
  • JPG, PNG etc. can be included in PDF
• GIF is only good for animations of line graphics - otherwise use PNG or JPG instead !
  • PNG or JPG have more/better colors (GIF has only 256 colors instead of the up to 16million 'true colors' of PNG or JPG)
  • compression is not as good as with PNG or JPG
  • a gif image may have one transparent color i.e. every pixel with this color are totally transparent
  • animation is the only advantage

The plot procedure can work with polar coordinates - you need to set the /polar keyword

; define some data points
N = 100
theta = findgen(N)/N*360.
r_max = 7.4 ;-)
radius = r_max*findgen(N)/N

plot, radius, theta*!DTOR, $
  /polar, $
  xstyle = 1, $
  xrange = [ -r_max, +r_max], $
  ystyle = 1, $
  yrange = [ -r_max, +r_max], $
  /isotropic

If you do not fix range and isotropy of the plot it may look weird due to automatic scaling.

You may want to have axes through the center:

plot, radius, theta*!DTOR, $
  /polar, $
  xstyle = 1+4, $ ; exact scaling but invisible (box)axis
  xrange = [ -r_max, +r_max], $
  ystyle = 1+4, $ ; exact scaling but invisible (box)axis
  yrange = [ -r_max, +r_max], $
  /isotropic

; Draw x- and y-axis through 0,0 with ticks to the top or right, respectively.
AXIS, 0, 0, XAXIS=0
AXIS, 0, 0, YAXIS=0

You may want to have a radius and angle grid. IDL cant do this for you. You have to do it step by step:

plot, radius, theta*!DTOR, $
  /polar, $
  xstyle = 1+4, $ ; exact scaling and invisible (box)axis
  xrange = [ -r_max, +r_max], $
  ystyle = 1+4, $
  yrange = [ -r_max, +r_max], $
  /isotropic

; draw and anotate spokes
N_spk = 8
for i = 0, N_spk-1 do begin
  wi = (2*!pi/N_spk)*i
  oplot, /polar, [0,90], replicate( !pi/2-wi, 2 ), linestyle=1
  r = 90.
  xyouts, r*sin(wi), r*cos(wi), string(wi/!dtor,f='(I0)')
endfor ; i

; draw and annotate angle circles
N_w = 181
w = findgen(N_w)/(N_w-1)*2*!pi
N_r = 3
for i = 0, N_r-1 do begin
  ei = (90.*i)/N_r
  oplot, /polar, replicate( 90-ei, N_w), w, linestyle= i gt 0 ? 1 : -1
  xyouts, 0, ei, string(ei,f='(I0)')
endfor ; i

Some routines like XYOUTS and ARROW do not accept the polar keyword (or ignore it). You may convert your coordinates to cartesian as above or you can use the CV_COORD function to do this:

polar_coo = [[angle],[radius]]
xy_coo = CV_COORD( FROM_POLAR=polar_coo, /DEGREES, /TO_RECT )

A contour plot with color shading is made with:

; prepare number of levels 
min_data = min( data, /nan )
max_data = max( data, /nan )
N_iso = 10
; we need one more levels value because level defines interval boundaries, i.e. color[i] goes into interval level[i] ... level[i+1]
levels = findgen(N_iso+1)/N_iso*(max_data - min_data) + min_data
; prepare true colors : e.g. a grayscale from white to black
gray = long(255*(1-findgen( N_iso )/(N_iso-1)))
iso_colors = gray or ishft( gray , 8 ) or ishft( gray, 16 ) ; colors as true color 3byte/24bit code
; do the contour plot
contour, data, x_coo, y_coo, $ ; data nad coordinates
  levels = levels, $ ; levels for the colors
  c_color = iso_colors, $ ; 24bit colors
  /cell_fill, $  ; avoid problems with mssing data (nan), etc. (you may use /fill instead but ...)
  xtitle = 'x_coo (units)', $ 
  xstyle = 1, $ ; use exact range of x_coo for axis range
  ytitle = 'y_coo (units)', $
  ystyle = 1 ; use exact range of y_coo for axis range

Routine conotur takes all the keywords like plot and some more - see help for contour and the descritption of graphics keywords.

Because levels give the lower border of the interval for the colors we need one more level than colors (levels=findgen(N_iso+1)...). For values below level[0] nothing is drawn (no color !), for values above level[N_iso], i.e. the last level the first color is recycled. To get colors for the values below the lowest level and above the largest level you can add additional colors and levels to what is defined above:

; expand coloring levels above and below range(levels)
iso_colors = [ col_below_min, iso_colors, col_above_max  ]
levels = [     extrm_lev_min, levels,     extrm_lev_max  ]

you may want to use the most extreme values your system can provide

extrm_lev_min = -(machar()).xmax
extrm_lev_max = +(machar()).xmax

See function function machar.

contour, data, x_coo, y_coo, levels = levels

If you want iso-lines with annotations give the keyword /follow :

contour, data, x_coo, y_coo, $
  levels = levels, $
  c_annotation = string( levels , format = '(f5.1)' ), $
  /follow, $
  /overplot

Keyword c_annotations is used here toprovide a format - this is usually not necessary. You can define also different thicknesses, colors etc. with the other c_... key words of contour see help for contour.

You can use contour also to make a color bar:

A horizontal color bar at the bottom

; charsize in normalized coordinates for position 
ch = float(!D.y_ch_size) / !D.y_size
; the color bar
contour, [ [levels] , [levels] ]  , $ ; a two row array with just the levels from the contour plot above
  levels, [0,1], $ ; levels on the x-axis and 2 lines in y direction
  levels = levels, $ ; levels for the shading
  c_color = iso_colors, $ ; the colors itself
  /cell_fill, $
  xtitle = 'data', $
  xstyle = 1, $ ; scale x-axis precisely
  xticklen = 5*ch, $ ; ticks must be elongated
  ytitle = '', $ ; no title on y-axis 
  ystyle = 1+4, $ ; do not draw y-axis 
  charsize = 1.0, $ ; smaller charsize ...
  position = [ 0.7 , 2*ch, 0.95, 3*ch ], $ ; postion in lower right corner below axis ...
  /normal, $ ; positon in normalized coordinates
  /noerase ; do not erase beckground - plot over the contour plot

A vertical color bar at the right side:

; charsize in normalized coordinates for position 
ch = float(!D.y_ch_size) / !D.y_size
contour, transpose( [ [levels] , [levels] ] ) , $ ; a two column array with the levels from the contour plot
  [0,1], levels, $ ; two columns in x-direction and levels in y-direction
  levels = levels, $ ; levles for the shading
  c_color = iso_colors, $
  /cell_fill, $
  ytitle = '', $ ; title will be in second y-axis on right side
  ystyle = 1, $
  yticklen = -5*ch, $ ; longer ticks on this axis
  ytickname = replicate(' ', 10), $ ; no annotations on this side 
  xstyle = 1+4, $ ; do not draw x-axis
  charsize = 1.0, $ ; smaller charsize ...
  position = [ 1-4*ch , 7*ch, 1-3*ch, 1-4*ch ], $ ; position at right side ...
  /normal, $ ; position in normalized coordinates
  /noerase ; do not erase background - plot over the contour plot

; y-axis on the right side 
axis, yaxis=1, $
  ytitle='data', $
  ystyle = 1, $
  yticklen = -5*ch, $
  yrange=[ min(levels) , max(levels) ], $
  charsize = 1.0, $ ; smaller charsize ...
  /save
  

The position of the vertical color bar in this example is not perfect. You may have to define the position of the contour plot itself to provide space for the color bar.

In case of data defined on a irregular grid you may need a triangulation, i.e. a kind of ordering of the data in triangles. This triangulation allows you to do CONTOUR plot and further interpolation of this data to a regular grid with TRIGRID.

lets Assume you have data Z at N points defined by coordinates X_coo and Y_coo. To get the triangulation just use:

TRIANGULATE, X_coo, Y_coo, triangles_idx, boundary_idx

triangles_idx is then of type lonarr(3,N_tri) containing the indices to the coordinates of the points defining the N_tri triangles.
boundary_idx is the list of indices of the coordinates of points of the convex boundary.
If the boundary of the region covered by your points has concave stretches TRIANGULATE produces triangles which cover them.
I found no parameter or function to avoid this or remove them.
If your coordinates are given in different, non isotropic units you should make them isotropic (e.g. X_coo as geographic coordinates like degrees east, and y_coo in meter above see level → use TRIANGULATE, X_coo*deg_2_meter, Y_coo, triangles_idx )

The output can be fed into CONTOUR with

CONTOUR, Z, X_coo, Y_coo, triangulation=triangles_idx, ...

You can interpolate triangulated data to a regular grid with TRIGRID:

If you just want a grid with Nx x Ny points

Z_grid = TRIGRID(  x_coo, y_coo, Z, triangles_idx, Nx=Nx, Ny=Ny, xgrid = x_grid, ygrid = y_grid, missing=!Values.F_NAN )

Z_grid is then fltarr(Nx,Ny) and y_grid and y_grid are of type fltarr(Nx) and fltarr(Ny).

If you want a section defined by its limits and a step size

grid_spacing = [ dx, dy ] 
grid_limits = [ x_min, y_min, x_max, y_max ]
Z_grid = TRIGRID(  x_coo, y_coo, Z, triangles_idx, grid_spacing, grid_limits, xgrid = x_grid, ygrid = y_grid, missing=!Values.F_NAN )

You can also specify coordinates x_nodes and y_nodes to which to interpolate:

Z_grid = TRIGRID(  x_coo, y_coo, Z, XOUT = X_nodes, YOUT = Y_nodes, missing=!Values.F_NAN )

Sometimes it might be preferable to have just colored tiles around the nodes of your data. Here is an example how to do it: The coordinate system must be defined - you may use the contour command as above with the /nodata keyword

; field sizes
Nx = n_elements(x_coo)
Ny = n_elements(y_coo)
; predefine neighbour nodes
x2 = x_coo[0]
x1 = x2
; intermediate border
x12 = (x1+x2)/2.
for ix = 0, Nx-1 do begin
  ; shift through nodes and borders ...
  x0 = x1
  x1 = x2
  x01 = x12
  ; new node and border
  x2 = x_coo[min([ix+1,Nx-1])]
  x12 = (x1+x2)/2.
  ; predefine ...
  y2 = y_coo[0]
  y1 = y2
  y12 = (y1+y2)/2.
  for iy = 0, Ny-1 do begin
    ; shift through ...
    y0 = y1
    y1 = y2
    y01 = y12
    ; new node ...
    y2 = y_coo[min([iy+1,Ny-1])]
    y12 = (y1+y2)/2.
    ; value
    zij = data[ix,iy]  
    ; color index
    if zij le min_data then i_col = 0 else $
    if zij le max_data then i_col = fix(float(zij-min_data)/(max_data-min_data)*(N_iso-1)) $
                       else i_col = N_iso-1
    ; color from color index
    color = iso_colors[i_col]
    ; fill the rectangle
    polyfill, [ x01, x01, x12, x12 ] , [ y01, y12, y12, y01 ], color = color
  endfor ; iy
endfor ; ix 

IDL can draw arrows (or vectors). The basic command is arrow which takes the starting point and end point of a (or several) vectors. It can take coordinates for the data-, normalized- or device system.

arrow, x0,y0, x1,y1, [, color=...][, thick=...][ ,/data | /normal | /device ][, ...]

Allthough it accepts the /polar keyword it is ignored. If you have vectors in polar coordinates you will need to convert them with cv_coord to cartesian coordinates (=rectangular).
If you have a non isotropic coordinate system with different scaling and lengths of the x-axis and y-axis you need to correct for this:

• we want to have in x and y the same device coordinates for an amplitude a in u or v-direction
• to achieve this we mutliply the v component of the vector with a factor fac:
• ⇒ dx_dev = a * !X.S[1]*!D.x_size = a * fac * !Y.S[1]*!D.y_size = dy_dev
• if we want to have dx_dev = dy_dev we need
• fac = !X.S[1]*!D.x_size / (!Y.S[1]*!D.y_size)

Within factor the ratio !X.S[1]/!Y.S[1] corrects for the different scaling of the axes and the ratio !D.x_size/!D.y_size corrects for the different lengths of the axes on the paper.

Example for N vectors (u[i],v[i]) at N positions x[i], y[i]

factor = !X.S[1] / !Y.S[1] * !D.x_size / !D.y_size
arrow, x, y, x+u, y+v*factor, /data

Similar with function velovect:

plot, x,y, /nodata
factor = !X.S[1] / !Y.S[1] * !D.x_size / !D.y_size
velovect, u, v*fac, x, y, /overplot

Note that you first have to define the coordinate system with e.g. the plot command such that !X and !Y are correctly defined.

IDL knows one and two dimensional histogram.

A one dimensional histogram is e.g. generated with

histo = histogram( data , binsize=bin, locations = classes , /nan )

the classes variable will contain the lower borders of the bins with the first bin starting at min(data).
(histogramm is made via histo=lonarr((max-min)/binsize+1)& for i do ii=(data[i]-min)/binsize)& data[ii] += 1& ...

A quick and dirty plot can be generated with psym=10:

plot, classes, histo, psym=10

but as classes contains the lower borders of bins it is shifted a half binsize to the left and the plot begins and ends with an open box. Slightly better is to shift the classes and add zero values at begin and end.

N_classes = n_elements(classes)
classes_x = [ classes[0]-bin , classes, classes[N_classes-1]+bin ] + bin/2
histo_x = [ 0 , histo, 0 ]
plot, classes_x , histo_x , psym=10, $
  xtitle = 'class (unit)', $
  ytitle = 'freq.'

or you write an own plot-histo procedure … or you ask jan …

Assume you have two parallel data series data_1 and data_2 and want a joint histogram :

; set bin-size, min and max and vector with bin-borders for data set 1
bin_size_1 = 0.1 ; give an adequate value
min_data_1 = min(data_1, /nan )
max_data_1 = max(data_1, /nan )
N_bins_1 = fix((max_data_1-min_data_1)/bin_size_1)+1
bins_1 = min_data_1 + (findgen(N_bins_1) + 0.5) * bin_size_1 ; output of histo refers to lower border of bins => add 1/2binsize to center them in bin
 
; set bin-size, min and max and vector with bin-borders for data set 2
bin_size_2 = 100. ; give an adequate value
min_data_2 = min(data_2, /nan )
max_data_2 = max(data_2, /nan )
N_bins_2 = fix((max_data_2-min_data_2)/bin_size_2)+1
bins_2 = min_data_2 + (findgen(N_bins_2) + 0.5)*bin_size_2
 
; make 2d Histogramm
histo_2d = hist_2d( $
  data_1[*] , data_2[*] , $ ; using [*] allows you to put together multidimensional data as eg. images ...
  bin1 = bin_size_1, min1 = min_data_1, max1 = max_data_1, $ ; if there are nans in data idl complains bin1 is not allowed => provide min and max to circumvent
  bin2 = bin_size_2, min2 = min_data_2, max2 = max_data_2 $
  )
 
; plot it
contour, histo_2d, bins_1, bins_2, $
  Nlevels = 20, $ ; some value
  /follow, $ ; annotate iso lines
  ; /cell_fill, $ ; or use color shading
  xtitle = 'data 1', $
  xstyle = 1, $
  ytitle = 'data 2', $
  ystyle = 1

Assume you have a time series and want a 2-D-histogram with values versus hour of day

; get hour from time (given as julian date)
caldat, time, month, day, year, hour, minute, second
hour_float = hour + (minute+second*60.)*60.
; alternatively
; hour_float = ((time-0.5-long(time-0.5)))*24.0

; get data 
data = fltarr( N_time )
...

; define bin sizes
bin_data = 100.0
bin_hour = 1.0
 
; construct classes
data_min = min( data )
data_max = max( data )
N_data_class = floor((data_max-data_min)/bin_data)+1
class_data = data_min + findgen(N_data_class)*bin_data + bin_data/2

; make 2d Histogramm
histo_2d = hist_2d( data , hour_float , bin1 = bin_data, bin2 = bin_hour )

To display a 2 dim function f(x,y) as 3-dim surface as a grid mesh use Surface shade_surf:

  ; generate function data
  Nx = 100
  Ny = 100
  ; x and y = -1..1
  x = 2*findgen(Nx)/Nx-1
  y = 2*findgen(Ny)/Ny-1
  ; N*N arrays with x and y - coordinates
  xx = x  #  replicate(1,Ny) ; # is matrix product, this gives a Nx*Ny matrix with x increasing along x direction
  yy = replicate(1,Nx)  #  y ; # is matrix product, this gives a Nx*Ny matrix with y increasing along y direction
  ; 2 dim radius array
  rxy = sqrt( xx^2 + yy^2 )
  ; a nice peaked function 
  fxy = cos(rxy*2*!pi*3) * exp(-(rxy/0.3)^2)
 
  ; plot in 3D ...
  surface, fxy,x,y, xtitle = 'x', ytitle='y', ztitle='z'

IF you want to have a solid, shaded surface use shade_surf instead:

  shade_surf, fxy,x,y, xtitle = 'x', ytitle='y', ztitle='z'

To plot a 3dimensional isosurface of a 3D-field use shade_volume to get the surface and polyshade to draw it:

  ; calucalte isosurface from data = fltarr( Nx, Ny, Nz )
  shade_volume, data, iso_value, vertices, indices
 
  ; convert vertices given in index coordiantes to data coordinates
  vertices[0,*] = xcoo[ vertices[0,*] ]
  vertices[1,*] = ycoo[ vertices[1,*] ]
  vertices[2,*] = zcoo[ vertices[2,*] ]
 
  ; set up 3D transformation 
  SCALE3, XRANGE=[xmin,xmax], YRANGE=[ymin,ymax], ZRANGE=[zmin,zmax]
 
  ; draw x-y-z axes and ...
  ; ... use data from bottom of volume for basic zero plane
  contour, $
      data[ * , * , 0 ], $
      xcoo, $
      ycoo, $
      zvalue = 0.0, $ ; relative position on z-axis
      nlevels = 255, $
      /fill , $
      /cell_fill, $ ; is necessary otherwise contours are not closed or reach diagonally over the plane ...
      xtitle='x (km)', $
      ytitle='y (km)', $
      ztitle='z (km)', $
      zrange=[ zmin, zmax ], $ ; otherwise z will be rescaled to to the range of data[*,*,0]
      /t3d, $ ; use the 3dim transformation defined above
      charsize=2.5
 
  ; draw iso surface
  img = polyshade( vertices , indices , /T3D )
 
  ; save to png file
  save_png, 'iso_surf.png', img

If no device was established img will contain the image data. If the z-buffer is used (set_plot,'z') variable img will contain a dummy value and the image can be get by use of tvrd(). If you are using the z-buffer several subsequent calls to polyshade and contour (with keyword /overplot) can be made.

To display 3 dimensional data with semi transparent coloring use function voxel_proj.

There are several more functions see 3D Visualization in the Categorial list of functions.

IDL uses the system matrix !P.T to describe the 3D-projection. It can be manipulated with the functions

• surfr to define the viewing angles as rotation angles aroud the axes
• scale3 to define scaling along x,y,z and z-axis and to define the viewing angles
• t3d reset, set and manipulate the transformation matrix and its parts

There is an article Perspective projection in the english wikipedia. The basics are:

We want to map an object point in object space A=(Ax,Ay,Az) to an image point on the screen B=(Bx,By,Bz) where Bz is its position relative to the screen (behind or in front - important for hidden lines).

The observer or camera is at C=(Cx,Cy,Cz), its viewing direction is defined by the 'eulerian' angles t=(tx,ty,tz) and is position relative to the screen by E=(ex,ey,ez) (i.e. E is given in the coordinate system of the camera).

The following steps are performed:

• D' = A-C position of A relative to the camera
• D = T*D' rotate D into the camera system with rotation matrix T = f(tx,ty,tz) = [ [ Txx, Txy, Txz ] , [ Tyx , … ] , … ]
• B = (D-E)*ez/dz project D onto the screen plane

Idl system matrix !P.T is a 4*4 Matrix describing the operation B' = !P.T * (A-C)

Bx'    =  (  Txx  ,  Txy  ,  Txz    ,  -ex   )  *  (  Ax-Cx  )
By'    =  (  Tyx  ,  Tyy  ,  Tyz    ,  -ey   )  *  (  Ay-Cy  )
Bz'    =  (  Tzx  ,  Tzy  ,  Tzz    ,  -ez   )  *  (  Az-Cz  )
wz'    =  (   0   ,   0   ,  -1/ez  ,   1    )  *  (    1    )

or something like that … (there is something wrong with -ez and -1/ez)

The image coordinates are calculated then by Bx = Bx'/wz'.

The rows of rotation matrix T can be interpretad as the unit vectors describing the orientation of the screen in the object coordinate system.

PDF supports in principle 3D grafics - i.e. the PDF standard describes a format to incorporate 3D models into PDF. This format is supported by the ADOBE reader since version 8 or 9, by other readers usually not. Examples can be found in this idl thread of compgroups.
http://dx.doi.org/10.1111/j.1365-2966.2011.19900.x
http://dx.doi.org/10.1038/nature07609
http://dx.doi.org/10.1126/science.1220574
A tutorial how to produce such 3D grafics can be found at www.astrobetter.com. Three 3 steps are suggested:

1. produce a 3D representation of the data in .obj code i.e. a (ascii) text representation
2. convert the .obj files into universal 3D a XML standard for hierarchical 3D data
3. generate a PDF file with LaTeX and PDFLaTeX and embed the .u3d data by using the movie 15 package

current time as string in the standard form of your system

t = SYSTIME()

current time in UTC as string …

t = SYSTIME( /UTC )

current time in Julian Day (Days (in fractions) since November 24, 4714 BC, 12:00:00 UTC)

t = SYSTIME( /JULIAN, /UTC )

Actually /julian alone gives not UTC - you have to add the flag /UTC

current time in unix epoch (seconds since 1.1.1970 00:00:00 UTC - watch out for Jan 19, 2038)

t = SYSTIME(/seconds)

this time contains fractions of seconds, can be used to find out how long program steps take:

t0 = SYSTIME(/sec)
... do something ...
t1 = SYSTIME(/sec)
print, 'this took', t1-t0,'sec', f='(A,f0.3,A)'

current time as string formatted humand readable downto milliseconds:

t = STRING( SYSTIME(/JULIAN), FORMAT='(C(CDI02,".",CMoI02,".",CYI04,"/",CHI02,":",CMI02,":",CSF06.3))')

Actually SYSTIME(/JULIAN) makes only steps of full seconds, accordingly the fractional seconds will remain zero. But we can use the (unix-)seconds output, convert it to julian and print:

sec_per_day = 24d*60d*60d
ref_day = julday(1,1,1970,0,0,0)
t = STRING( ref_day+SYSTIME(/SEC)/sec_per_day, FORMAT='(C(CDI02,".",CMoI02,".",CYI04,"/",CHI02,":",CMI02,":",CSF06.3))')

Julian date = time in days since 1. Januar −4712 (4713 BC) 12:00 TDT ( = UTC + 1min 5.184sec )
see http://en.wikipedia.org/wiki/Julian_date or http://de.wikipedia.org/wiki/Julianisches_Datum

Julian day is a continous time variable and can easily be used to

• calulate time differences (dt = t2-t1 in days, dt=(t2-t1)*24. in hours … )
• print formated times (see below)
• annotate grafics

julian day is the IDL format for time. It is stored as a double variable counting the days.

Be aware that zero fractional part refers to 12:00 (as the start time of julian day is 12:00)

To create a time in julian date use function julday:

time = julday( month, day, year, hour, minute, second )

• You can omit the time (hour, minute, second).
• If you provide at least the hour you get a double precision variable, if not, you get an integer variable of type long.
• If one of the input variables are arrays you get an array.
• If these arrays have different size IDL uses the smallest size and gives no(!) warning.
• Zero input values for all except year are allowed
• If year is zero you get an error message “there is no year zero in the civil calendar”
• Negative input values are allowed.
• You can provide strings as input, but be aware that empty strings are interpreted as zero.

To calculate a date use procedure caldat :

caldat, time, month, day, year, hour, minute, second

• if time is an array output month, day, etc will also be arrays

print formatted date from julian day in two steps:

CALDAT, t, M,D,Y,HH,MM,SS
PRINT, D,'.',M,'.',Y,' ',HH,':',MM,':',SS,FORMAT="(i2,a,i2,a,i4,a,i2,a,i2,a,i2)"

or more elegant in one step for DD.MM.YYYY hh:mm:ss :

PRINT, t, FORMAT='(C(CDI02,".",CMOI02,".",CYI04," ",CHI02,":",CMI02,":",CSI02))'

or for YYMMDDhhmmss :

PRINT, t, FORMAT='(C(CYI02,CMOI02,CDI02,CHI02,CMI02,CSI02))'

you can also have fraction of a second for hh:mm:ss.ttt i.e. seconds downto 1/1000.sec use:

PRINT, t, FORMAT='(C(CHI02,":",CMI02,":",CSF06.3))'

the function STRING can be used as well; for DD.MM.YYYY hh:mm:ss use:

s = STRING( t, FORMAT='(C(CDI02,".",CMOI02,".",CYI04," ",CHI02,":",CMI02,":",CSI02))' )

for YYMMDDhhmmss use:

s = STRING( t, FORMAT='(C(CYI02,CMOI02,CDI02,CHI02,CMI02,CSI02))' )

The format specifiers are construced as:

C<field identifier><format specifier>[<leading zero>]<number of digits>
* <field identifier> can be D for day, MO for month, Y for year, H for hour ...
* <format specifier> can be F for float, I for integer, A for character
* <leading zero> give a zero for leading zeros
* <number of digits> 2 makes a field with a width of 2 digits/characters

(see Help for "format codes" -> "C() format codes" )

The string-formatting routine does not round the fields, but instead delivers the floor of each field. This may result in e.g. <hh-1>:59 if you ought to have hh:00 but have a 1/10000s missing. To avoid this add half of the smallest unit you do print and round times:

• for hh:mm use print, round(time*minperday+0.5)/minperday, format='(C(CHI02,":",CMI02))'
• for hh:mm:ss use print, round(time*secperday+0.5)/secperday, format='(C(CHI02,":",CMI02,":",CSI02))'
• for hh:mm:ss.ttt use print, time, format='(C(CHI02,":",CMI02,":",CSF6.3))' (fractions of a second are rounded ???)

with

minperday = 24d*60d = 1440.
secperday = 24d*60d*60d = 86400.

Julian day comes in type DOUBLE-Precision which has a mantissa of 52 bits. The smallest resolvable step is thus dt_rel=1/2^52=2.2E-16. As we have currently (2019) Julian day number t=2458695 this is dt_min=dt_rel*t*sec_per_day=4.7E-5sec=0.047msec (sec_per_day=24.*60*60=86400). You can test this with this code ():

spd=86400d&t0=julday(1,1,1970,0,0,0)&twt=1e-6&for i=0,99 do begin&t1=systime(/sec)/spd+t0&wait,twt&t2=systime(/sec)/spd+t0&print,(t2-t1)*spd&end

It will most of the time print 0.0000 and sometimes 4E-5. If you lower twt (t_wait) it will always print 0.0000

'human' date to julian date

t = JULDAY(Month, Day, Year [,hour[,minute[,second]]] )

if you omit hour, minute, second you will get type long, otherwise double

string to julian date for a string of the form DD.MM.YYYY hh:mm:ss:

; DD.MM.YYYY hh:mm:ss
; 0123456789 12345678
;   JULDAY(   month      ,   day        ,   year       ,    hour       ,    minute     ,    second      ) 
t = JULDAY( strmid(s,3,2), strmid(s,0,2), strmid(s,6,4), strmid(s,11,2), strmid(s,14,2), strmid(s,17,2) )

Time differences can be easily calculated - but be aware that julian date starts at 12:00 - if you want display time differences as time strings you must subtract a half day:

PRINT, t1-t2 - 0.5 , FORMAT='(C(CHI02,":",CMI02,":",CSI02))'

To annotate a plot axis with formatted time strings just set an appropriate format specifier in the plot command:

plot, time, data, ..., xtickformat = '(C(CHI02,":",CMI02))', ...

where time is given in julian day …

Unix time are the seconds since 1.1.1970 00:00:00 UTC
Sometimes it is also called Epoch (see http://en.wikipedia.org/wiki/Unix_time)

Originally t_unx is defined as long integer - this leads to the year 2038 Problem

Julian day (t_jul) to unix time (t_unx):

secperday = 24*60*60d
unx_start = JULDAY( 1, 1, 1970, 0, 0, 0 ) ; this is equal to 2440587.5
t_unx = (t_jul - unx_start) * secperday

and vice versa:

t_jul = t_unx/secperday + unx_start

If you have a time t_xxx as tics since t_ref (e.g. (fractional)number of 7min intervals since 3.4.1722 13:21 (just to show that it can be really weird ) then define:

tics_per_day = min_per_tic / min_per_day  or  sec_per_tic/sec_per_day  etc...
t_ref = JULDAY( mo_ref, dd_ref, yy_ref,  HH_ref, MM_ref, SS_ref )

convert julian day (t_jul) time into this exotic format (t_xxx)

t_xxx = (t_jul - t_ref) * tics_per_day

and vice versa:

t_jul = t_xxx/tics_per_day + t_ref

Examples are:

• excel time 1 = days since 1.1.1900 (nineteen hundred)
  (excel has two internal timeformats, you can switch somewhere inside excel between both)
  • t_ref = julday( 1,1,1900, 0,0,0 )
  • tics_per_day = 1.0
  • typical values (june 2023) are in the order of 45087 (45tsnd) days
  • Excel cannot deal with negative times, i.e. dates before t_ref are not valid
  • if you export excel tables as csv you might get this as time variable
• excel time 2 = days since 1.1.1904 (nineteen O four! )
  • t_ref = julday( 1,1,1904, 0,0,0 )
  • tics_per_day = 1.0
  • typical values (june 2023) are in the order of 43627 (43tsnd) days
  • if you export excel tables as csv you might get this as time variable
  • you can switch somwhere inside excel between both time systems - be careful
  • Excel cannot deal with negative times, i.e. dates before t_ref are not valid
  • Excel has a known bug in the year 1900 it thinks it errornously assumes it was a leap year ⇒ all excel dates before 1.3.1900 might be wrong (weekdays, differences etc.)
• unix time = seconds since 1.1.1970
  • t_ref = julday( 1,1,1900, 0,0,0 )
  • tics_per_day = 86400 (=24*60*60)
  • typical values (june 2023) are in the order of 1686589673 = 1.67E9 sec
  • this ist the internal time format on unix systems and derivatives (e.g. android)
  • unix time may experience an overflow or a roll over in 2038 when range of long32 is exceeded. This may lead to problems in the computer world.
• days since 1.1.2000 ( time format in RPG data ? )
  • t_ref = julday(1,1,2000,0,0,0)
  • tics_per_day = 1.
  • typical values (june 2023) are in the order of 8563 days
• seconds since 1.1.2000 ( time format in Juelich radiation data )
  • t_ref = julday(1,1,2000,0,0,0)
  • tics_per_day = 86400.
  • typical values (june 2023) are in the order of 739905199 = 7.3990510e+008 = 740Mio seconds

The day of the year (doy) - sometimes incorrectly named 'Julian day' - can be calulated from a Julian date by:

CALDAT, t_julian, month, day, year
t_doy = floor( t_julian - JULDAY( 1,1,year ) ) + 1

Note: Day of year starts with 1 on January first ⇒ add one. The simplified formula doy = t_julian - floor(t_julian/days_per_year)*days_per_year - with days_per_year=365.25 does not work - why ?

The hour of the day (hod) is the fractional part of the date as julian day shifted by half a day/twelve hours:

t_hod = ( t_julian-0.5 - floor( t_julian-0.5 ) ) * 24.

or alternatively and more intuitive, but slower

CALDAT, t_julian, month, day, year, hh, mm, ss
t_hod = hh+( mm + ss/60.)/60.

IDL has some predefined functions for interpolation:

• interpol interpolates a 1dim vector to new abscissa values - either on a regular grid or with given abscissas
• interpolate can do bilinear interpolation on two-dimensional regular or irregular grids
• bilinear can do bilinear interpolation on regular grids
• krig2d does 2 dimensional Kriging i.e. a distant weighted interpolation from neighboring points.

The interpolate routine takes only floating point index values for the interpolation. If you want to interpolate to certain nodes x_ipol and y_ipol you have to calculate the respective indices. This can be done by the interpol routine as follows. Assume we have a data array z_org = fltarr(Nx_org,Ny_org) with coordinates x_org=fltarr(Nx_org) and y_org=fltarr(Ny_org) and want interpolate it to z_ipol = fltarr(Nx_ipol,Ny_ipol) with coordinates x_ipol=fltarr(Nx_ipol) and y_ipol=fltarr(Ny_ipol):

; bilinear interpolation 
ix_ipol = interpol( findgen(Nx_org) , x_org , x_ipol )
iy_ipol = interpol( findgen(Ny_org) , y_org , y_ipol )
z_ipol = interpolate( z_org , ix_ipol , iy_ipol , /grid  )

If you want to interpolate only over certain dimensions they must be the last ones. You may use transpose to rearrange. Lets assume you have an array depending on x, y and time, i.e. z=fltarr(Nx,Ny,Nt) and you want to interpolate to certain x and y coordinates you can do this as follows:

z_ipol = transpose( interpolate( transpose(z_org,[2,0,1]) , ix_ipol , iy_ipol , /grid  ), [1,2,0] )

IDL can transform geographical coordinates into cartesian x,y coordinates and vice versa. To do that you have to define your coordinatesystem with MAP_PROJ_INIT. To calculate your projected coordinates (X,Y) from longitude and latitude use MAP_PROJ_FORWARD and to get longitude and latitude for a (X,Y) pair use MAP_PROJ_INVERSE.

Example for the universal transverse mercator projection (UTM) for West Germany (UTM zone 32 with central meridian at 9°East):

Init the transformation:

  map_struct = MAP_PROJ_INIT( $
    101 , $ ; UTM
    CENTER_LONGITUDE = longitude , $
    )

The result map_struct contains information of the projection and must be handed to the functions for the caluclation of coordinates. The parameter CENTER_LONGITUDE must be given to define which UTM grid to be used (read wikipedia if you want to understand the details). It should be a longitude close to the locations you are interested in. Following the UTM definition IDL will round it to a multiple of 6° relative to 3° (i.e. CENTER_LONGITUDE=3+n*6 with n an integer). Calcualted UTM coordinates will be conform to the UTM specifaications but the direction of UTM-northing may deviate by several degrees from true north. If you need a coordinate system directed exactly to north at your location use instead:

  map_struct = MAP_PROJ_INIT( $
    109 , $ ; transverse mercator
    DATUM = 8 , $ ; WGS84
    CENTER_LONGITUDE = longitude , $
    FALSE_EASTING = 500000. $
    )

Calculate some UTM values from longitude long and latitude latt:

  UTM = MAP_PROJ_FORWARD( long , latt , map_structure=map_struct  )

where long and latt can be arrays (fltarr(N)) and UTM will be of type fltarr(2,N) containing the calculated UTM coordinates.

The inverse: calculate some longitude and latitude values from UTM=fltarr(2,N)

  LonLat = MAP_PROJ_inverse( UTM, map_structure=map_struct )

A list of the ID's of the projections can be found at MAP_PROJ_INIT help article.

You can also generate a map and draw into this map by using MAP_SET :

  center = [ center_lat , center_lon ]
  d_ll = 30.
  map_box = [ center - [1,1]*d_ll , center + [1,1]*d_ll ]
 
  MAP_SET, /Azimuthal, /isotropic, $
    center_lat , center_lon , 0 , $
    limit = [ lat_a, lon_a, lat_b, lon_b ]
 
  ; project an image to the map with img=bytarr(Nx,Ny)
  read_jpeg, '/data/sfb1211/A01/srtm/other/images/Envisat_global_land_cover_map.jpg', img
  latmin_jpg = -94. &  latmax_jpg = +106. & lonmin_jpg = -180. & lonmax_jpg = +180.
  ; position to put image
  start_x = lonmin_map
  start_y = latmin_map
  ; wrap image - MAP_image can deal only with 2dim array - we need to do it channel by channel
  img_mapped_red = MAP_IMAGE( reform(img[0,*,*]) , startX, startY, lonmin=lonmin_jpg, lonmax=lonmax_jpg, latmin=latmin_jpg, latmax=latmax_jpg )
  img_mapped_grn = MAP_IMAGE( reform(img[1,*,*]) , startX, startY, lonmin=lonmin_jpg, lonmax=lonmax_jpg, latmin=latmin_jpg, latmax=latmax_jpg )
  img_mapped_blu = MAP_IMAGE( reform(img[2,*,*]) , startX, startY, lonmin=lonmin_jpg, lonmax=lonmax_jpg, latmin=latmin_jpg, latmax=latmax_jpg )
  ; put channels togehter to one image
  img_mapped_rgb = transpose( [ [[img_mapped_red]], [[img_mapped_grn]], [[img_mapped_blu]]], [2,0,1] )
  ; put image on canvas
  TV, img_mapped_rgb, startX, startY, true=1
 
 
  ; draw coastlines, continents and Country borders in high resolution 
  MAP_CONTINENTS, /COASTS, /CONTINENTS, /COUNTRIES, /HIRES
 
  ; rivers - idl draw allways all rivers - an hierarchy as in GMT is not there (https://docs.generic-mapping-tools.org/6.3/pscoast.html#i)
  MAP_CONTINENTS, /RIVERS, /HIRES, color = 'ff0000'x
 
  ; draw a grid, make nice black and white axes around it
  MAP_GRID , /box_axes
 
  ; plot and annotate some data (lat=fltarr(N), lon=fltarr(N), name=strarr(N))
  oplot, lon, lat, psym=4
  for i = 0, n_elements(lat)-1 do xyouts, lon[i], lat[i], name[i]
 

IDL's atan function can be fed with x and y component of a vector to get the direction angle from -pi to +pi = -180° to +180°.

phi = atan( x, y )

It performs the complete angle determination based on the signs of x and y, i.e. something like

function atan_2, y, x
  ; determination of azimuth angle from x and y components
  ; this is principally done if you call atan( y,x )
if  x eq 0. then begin
  if y gt 0. then phi =  !pi/2. else $
  if y lt 0. then phi = -!pi/2. $
             else phi =  0.0
endif $; x eq 0
else begin ; x ne 0
  phi0 = atan( y/x )
  if (y gt 0. ) and (x gt 0.) then phi = phi0       else $ ; first  quadrant
  if (y le 0. ) and (x gt 0.) then phi = phi0       else $ ; second quadrant
  if (y lt 0. ) and (x lt 0.) then phi = phi0 - !pi else $ ; third  quadrant
  if (y gt 0. ) and (x lt 0.) then phi = phi0 + !pi        ; fourth quadrant
endelse 
return, phi
end ; function 

As the wind directon azimuth angle is differently defined than mathematical azimuth angle i give here the ultimative conversion laws between wind direction azi, wind speed spd and vector components u (W-E) and v (S-N):

; spd and azi  to  u and v
u = - spd * sin(azi*!dtor)
v = - spd * cos(azi*!dtor)

; u and v  to  spd and azi
azi = atan(  -u , -v ) / !dtor
if azi lt 0 then azi += 360.
spd = sqrt( u^2 + v^2 )

Factor !dtor converts degrees to radians, it is of type FLOAT (=single pecision) for DOUBLE precisions use !CONST.DtoR

Examples:

• North wind: azi = 0° ⇒ u = 0 (no W-E component), v = -spd (wind is blowing from north to south)
• East wind: azi = 90° ⇒ u = -spd (wind is blowing from East to West), v = 0 (no N-S component)
• West wind: azi =270° ⇒ u = +spd (wind is blowing from West to east), v = 0 (no N-S component)

The error can be estimated with gaussian error propagation:

<m>

\phi = arctan(-u/-v)
     = arctan(x)  \ ; \  
x = {-u}/{-v}

</m>

<m>

\Delta \phi^2 = ({d\phi}/du \Delta u)^2 + ({d\phi}/dv \Delta v)^2
             = ({d\phi}/dx dx/du \Delta u)^2 + ({d\phi}/dx dx/dv \Delta v)^2

</m>

if we consider <m>{d\phi}/dx=1/{1+x^2}</m> we get

<m>

\Delta \phi^2 = 1/(u^2+v^2)^2 ( (v \Delta u)^2 + (u \Delta v)^2 )

</m>

(or quantiles)

IDL can give you the median but not any other arbitrary percentile. So do it yourself:

; define 5 percentiles
N_perc = 5
p1 = 0.05 ; 5% and 95%
p2 = 0.25 ; 25% and 75% = lower and upper quartiles
; add the median and mirror p1 and p2:
p_fac = [ p1 , p2 , 0.5, 1-p2 , 1-p1 ]

; determine the percentiles
percentiles = (data[ sort[data] ])[ p_fac*(n_elements(data)-1) ]

This way +/-infinity and Nan values are included. IDL sorts them in the following way:

data[sort[data]] = [[-inf.values], [finite.values], [+inf.values], [Nan.values] ]

You have to count the respective non-finite values and include it in the index calculation:

N_neg_inf = long( total(finite(data,/infinity,sign=-1)) )
N_pos_inf = long( total(finite(data,/infinity,sign=+1)) )
N_nan     = long( total(finite(data,/nan             )) )
percentiles = (data[ sort[data] ])[ N_neg_inf + p_fac*(n_elements(data)-1-(N_neg_inf+N_pos_inf+N_nan)) ]

An implementation as function can be found in /home/hatpro/idl/lib/f_percentile.pro

There is no exact method to determine percentiles. Imagine a dataset with 50% of the values lower than a vlaue y1 and the other 50% larger than a value y2 with a large gap between y1 and y2.

• Where ist the median in this case? Is it y1, y2 or somewhere in between?
• Is interpolation between y1 and y2 allowed?
• In the same way any percentile may not be exactly determined
• wikipedia gives a method for interpolation between these values…

Random numbers can be generated with

x = RANDOMU( seed )
x = RANDOMU( seed, N )

for a single number or an array with N elements of Uniform distrinbuted random numbers in [0,1)
Expectation value is thus μ=0.5 and standard deviation σ=3.46.
(Uniform disributed radnom numbers have an expectation value μ=(a+b)/2 and standard deviation σ=sqrt(12)*(b-a) if they lie in the range [a,b))

x = RANDOMN( seed )
x = RANDOMN( seed, N )

for as single number or an array with N elements of Normal distrinbuted random numbers with expectation value μ=0 and standard deviation σ=1.0. There is theoretically no upper and lower border for these random numbers but |x-µ| will probably never be much larger than 4σ.

BUT there is no Zufall in computers !
All these routines produce pseudo random numbers
The random number is calculated from seed !
With every random number a new seed is returned for the next random number …
The following code will allways produce the same sequence of ten “radnom” numbers:

seed = 0
for i =0, 10 do PRINT, RANDOMU(seed)

(you may also omit seed = 0 as idl sets seed than internally to 0)

To generate some randomness you can use the time:

seed = SYSTIME(/seconds)
for i =0, 10 do PRINT, RANDOMU(seed)

See Numerical recipies chapter 7 or wikipedia.

every number can be displayed in many different ways. In IDL this canbe done with formatted output.
Let e.g. x = 16(decimal):

• binary: print, x, format='(B)' → 10000
• octal: print, x, format='(O)' → 20
• decimal: print, x, format='(I)' → 16
• hexadecimal: print, x, format='(Z)' → 10
• roman: print, ääääää … ok that does not work - but it should look like XVI

The principle of all these number representations (except for the roman) is positional notation (wikipedia):
A number written as a sequence of digits A_i:
x = A₃A₂A₁A₀
is mathemtically the sum:
x = A₃*P³ + A₂*P² + A₁*P¹ + A₀*P⁰

• Where P is the base of the used number representation.
  • in the decimal system P=10
  • in the binary system P=2
  • in the octal system P=8
  • in the hexadecimalsystem P=16
  • …
• Digits A_i can lie in the interval 0 to P-1.
  • in the decimal system (P=10) valid digits are 0..9
  • in the binary system (P=2) valid digits are 0..1
  • in the octal system (P=8) valid digits are 0..7
  • in the hexadecimalsystem (P=16) valid digits have the decimal values 0..15 - but since digits should be only one character wide the values 10 to 15 are replaced by the first six letters of the alphabet: 10→A, 11→B, 12→C, 13→D, 14→E, 15→F
• Every position i in the written number is related to the base value Pⁱ
  • the rightmost position in a decimal number has the base value 1 = 10⁰,
    the position to the left has the base value 10 = 10¹
    the position to the left has the base value 100 = 10²
    …
  • the rightmost position in a binary number has the base value 1 = 2⁰,
    the position to the left has the base value 2 = 2¹
    the position to the left has the base value 4 = 2²
    …
  • the rightmost position in a octal number has the base value 1 = 8⁰,
    the position to the left has the base value 8 = 8¹
    the position to the left has the base value 64 = 8²
    …

Since computers do internally work with binary numbers it migth be sometimes useful to write binary numbers - but this creates easily very long numbers. That is the reason to use octal or hexadecimal numbers. Their base is a power of 2 and thus compress some binary digits into one digit. In the octal system 3 binary postions are compressed into one (P_oct=8=2³) and in the hexadecimal system 4 binary postions are compressed into one (P_hex=16=2⁴):

bin   oct      bin  oct ...   bin  hex(dec)    bin   hex(dec)
000    0       100   4       1000   8  (8)    1100    C  (12)
001    1       101   5       1001   9  (9)    1101    D  (13)
010    2       110   6       1010   A (10)    1110    E  (14)
011    3       111   7       1011   B (11)    1111    F  (15)

To convert a number x into any number system on the base P 'manually' (by mental arithmetic ?) divide the number by P and note the divisons rest (a number between 0 and P-1) as the rightmost digit. Continue the same with the result of the division and write the divisons rest to the left of your previous digit until you end up with zero …

To set IDL variables to certain binary, octal or hexadecimal values write the number in the desired representation as string followed by a letter inidcating the used number system:

x = '10000'b
x = '20'O
x = '10'x

To convert strings given in a number system to integer values use reads with an apropriate format specifier:

; binary
s = '10000'
x = 0L
reads, s , x , format='(B)'
print, x

; octal
s = '20'
x = 0L
reads, s , x , format='(O)'
print, x

; hexadecimal
s = '10'
x = 0L
reads, s , x , format='(Z)'
print, x

In many cases the status of a system is delivered as flags in a binary number. Every position (i.e. bit OR flag) in the binary number is related to a certain status of a subsystem. A 0 represents false or off or not active etc. A 1 (or a set bit) represents vice versa a true or on or active etc. To check these flags it is convenient to use 'binary' operators (here are operators meant dealing with binary numbers - the same word is used for operators between two numbers…):

• c = a AND b : if bit at position i is set in a and b it will be set in c too and otherwise not:
  a=1010_bin b=0001_bin ⇒ c=0000_bin=0_dec
  a=1010_bin b=0010_bin ⇒ c=0010_bin=2_dec
  …
• c = a OR b : if bit at position i is set in a or b it will be set in c and otherwise not:
  a=1010_bin b=0001_bin ⇒ c=1011_bin=11_dec
  a=1010_bin b=0010_bin ⇒ c=1010_bin=10_dec
  …
• c = a XOR b : if bit at position i is set in a only or b only it will be set in c and otherwise not (exclusive or):
  a=1010_bin b=0001_bin ⇒ c=1011_bin=11_dez
  a=1010_bin b=0010_bin ⇒ c=1000_bin=8_dec
  …
• c = NOT a : if bit at position i is not set in a it will be set in c and otherwise not:
  a=1010_bin ⇒ c=0101_bin=5_dez
  …

Dont mix these operators with logical operators

To check whether in a status variable stat several flags are set define a mask and check by comparing stat with mask using of the AND operator:

; define a mask to check for bits 2, 6, and 7
mask = '00000100'b OR '01000000'b OR '10000000'b
; check
IF (stat AND mask) ne 0 THEN BEGIN
  print, "some of the flags are set ..."
  ...
ENDIF

in rare cases it might be necessary to see the representation of a variable in memory. This can be done in the following way:

• copy content of variable into a long or long64 integer variable
• print it formatted in hex, octal or binary

see https://www.nv5geospatialsoftware.com/Support/Self-Help-Tools/Help-Articles/Help-Articles-Detail/ArtMID/10220/ArticleID/19447/2723

For single precision variables we need a 4Byte=32Bit integer variable i.e. type long. We print in Hexadecimal (formatcode Z), Octal (O) and Binary (B)

x = 1.234
xl = long(x,0)
print, x, xl, xl, xl, f='(f," ",Z," ",O," ",B)'

For double precision variables we need an 8byte = 8*8bit = 64bit integer variable, i.e. Long64:

x = 1.234d
xl = long64(x,0)
print, x, xl, xl, xl, f='(f," ",Z," ",O," ",B)'

This way we can see how NAN's and infinity are represented:

 print, long( !values.f_nan,0),      long64( !values.d_nan,0),      f='(" fnan=",Z,"  dnan=",Z)'
 print, long(-!values.f_nan,0),      long64(-!values.d_nan,0),      f='("-fnan=",Z," -dnan=",Z)'
 print, long( !values.f_infinity,0), long64( !values.d_infinity,0), f='("+finf=",Z," +dinf=",Z)'
 print, long(-!values.f_infinity,0), long64(-!values.d_infinity,0), f='("-finf=",Z," -dinf=",Z)'

will give:

 fnan=    7FC00000  dnan=      7FF8000000000000
-fnan=    FFC00000 -dnan=      FFF8000000000000
+finf=    7F800000 +dinf=      7FF0000000000000
-finf=    FF800000 -dinf=      FFF0000000000000

The -NaN values is a bit strange as NaN is Not-a-Number but IDL also produces NaNs where the sign flag is set.

A fourier transform is done with function FFT. to plot the spectrum of a time series y running over a time interval of length T do the following:

FT_y = FFT( y )

N = n_elements( y )
N2 = N/2
i2 = 1+findgen(N2-1)
freq = i2/N * 1/T
FTy_cos = real_part( FT_y[i2] )
FTy_sin = imaginary( FT_y[i2] )
FTy_pow = FT_y[i2] * conj(FT_y[i2])

plot, freq, FTy_cos, color='0000ff'x
plot, freq, FTy_sin, color='ff0000'x
plot, freq, FTy_pow

The power spectrum eventually needs some adjustment - depending on how you define it.

You can use FFT to calculate the auto-covariance and cross-covariance of two time series x and y:

ft_x = FFT( x )
ft_y = FFT( y )
auto_cov_x = FFT( ft_x * CONJ(ft_x) , /INVERSE )
cross_cov_xy = FFT( ft_x * CONJ(ft_y) , /INVERSE )

For long time series this is much faster than a classical sum( xi'*yi' ).

And you can do this also for two dimensional data sets x and y …

Wavelet analysis (WV) splits similar to [fast]fouriertransform (FFT) a dataseries into components at different scales (or wavelengths or frequencies). In difference to the FFT it gives also a localization of analysed scales, i.e. it gives information when or where in a dataseries a certain frequency had a a large amlitude.

Idl has a Wavelet Toolkit consisting of a number of graphical user interfaces and functions. It is described in the IDL help under Introduction to the IDL Wavelet Toolkit.

But you need to have a license for the toolkit to use the graphical user interface. We do not have the license!

Nevertheless the basic functions are accessible without the license. In the help there is a list of the functions . Documentation is rather poor sometimes it is helpful to study the source code of the functions you find in the IDL-lib/wavelet/source/ directory (that is in our implementation in /opt/itt/idl/lib/wavelet/source/). It seems as if all the functions were written by Chris Torrence in 2000. He refers in his documentation to

Torrence, C. and G.P. Compo, 1998: A practical guide to wavelet analysis. Bull. Amer. Meteor. Soc., v. 79, pp. 61-78. DOI:10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2

Basically you can do two kinds of wavelet transformations,
the so called discrete wavelet transform (DWT) and the continous wavelet transform (CWT):

• DWT is based on orthognal wavelets like the Haar or the Daubechie Wavelets. It splits the scale/space domain of the wavelet transform in separate regions with no overlapping information. It will deliver wavelet coefficients in the so called pyramidical order: two coefficents for the lowest scale, four for the second, … and N coefficent for the N'th scale.
  • call with
    wv = WV_DWT( data, Scaling, Wavelet, Ioff, Joff [, ... ] )
    • data = fltarray(N_data) input dataseries
    • scaling = fltarr(N_coeff) scaling coefficients of the wavelet, sometimes called the 'father wavelet', see wavelet functions below
    • wavelet = fltarr(N_coeff) wavelet coefficients, sometimes called the 'mother wavelet', see wavelet functions below
    • Ioff = scalar integer, offset of the supprot for the scaling function, set to -N/2+2 to center it over each point in data
    • Joff = scalar integer, offset of the supprot for the wavelet function, set to -N/2+2 to center it over each point in data
    • further details on WV_DWT
• CWT can be based on virtually any function (of course there are some restrictions) the best known are the morlet wave package or the mexican hat. The information received by the transformation has some overlap in the scale/space domain. It will deliver the wavelet coefficients as fltarr(N_data,N_scale)
  • call with
    wv = WV_CWT( data, Family, Order [, ... ] )
    • data = fltarray(N_data), input dataseries
    • Family = string, name of the wavelet function. Is expanded to WV_FN_family.PRO, see wavelet functions below
    • Order = float, parameter forwarded to the wavelet function, e.g. the number of zero intercepts of the morlet wavelet.
    • furhter details on WV_CWT

IDL provides seven wavelet functions in directory /opt/itt/idl/lib/wavelet/source/ with names

wv_fn_coiflet.pro
wv_fn_daubechies.pro
wv_fn_gaussian.pro
wv_fn_haar.pro
wv_fn_morlet.pro
wv_fn_paul.pro
wv_fn_symlet.pro

These functions return a structure with information about the wavelet. The values of the wavelet funcion itself come in parameter wavelet. Despite their similar names there are two different types of wavelet functions in IDL: DISCRETE and CONTINOUS:

• discrete are the coiflet, daubechies, haar and symlet wavelet
  • they can be used in the discrete wavelet transform WV_DWT only !
  • they are called by
    wv_info = WV_FN_<name>( [Order, Scaling, Wavelet[, Ioff, Joff]] )
  • the 'mother wavelet' in array wavelet is returned in the third parameter without keyname
  • Scaling is output parameter containing the 'father wavelet' as array
  • wavelet is the wavelet function
• continous are the gaussian, morlet and paul wavelet
  • they can be used in the continous wavelet transform WV_CWT only !
  • they are called by
    wv_info = WV_FN_<name>( [Order, Scale, N ][, /DOUBLE] [, FREQUENCY=variable] [, /SPATIAL] [, WAVELET=variable])
  • the 'mother wavelt' is returned to the variable specified with keyword WAVELET
  • scale is a scalar input parameter containing the scale of the function (its width if you want)
  • in wavelet is returned the fourier transform of the wavelet function, except when keyword /SPATIAL ist set
  • gaussian is a derivative of the gaussian bell where parameter 'order' is the order of the derivative: 1 is the first derivative, 2 is the second derivative (sometimes called the 'mexican hat wavelet')
  • morlet is a sine with damped amplitude parameter 'order' is called 'nondimensioanl frequency parameter' in the comments of the function. It is the ratio between amplitude damping constant and frequency. It is a measure for the number of full oscialltions.

In principle you can construct your own wavelet function for the continous wavelet transform by writing a function following the name and call conventions:

• the filename should be WV_FN_<name>.PRO where <name> is then used as paramter family when calling WV_CWT
• the source code can stand anywhere in the idl path, e.g. in the directory of your current project.
• for a continous wavelet you must provide parameters Order, Scale, N, and named paramters WAVELET, SPATIAL, DOUBLE and FREQUENCY.
• the function must return a structure containing information about the function - allthough it is not evaluated in WV_DWT.PRO
• if flag SPATIAL is set it must return the wavelet
• if flag SPATIAL is not set it must return the fourier transform of the wavelet and the related frequencies

But so far i get only weird results trying to implement a continous Haar wavelet …

This is a lot of code - mainly to reconstruct the values of scale and time for the wavelet coefficients.

; generate wavelet function
; info record for wavelet function
info_rec = WV_FN_HAAR( 1, scaling, wavelet, ioff, joff)
 
; ioff and joff need to be set !
ioff = 0 ; alternativley -N/2+2
joff = 0 ; alternativley -N/2+2
 
; Compute the discrete wavelet transform
wave_disc = WV_DWT( data, scaling, wavelet,   ioff, joff )
 
; power spectrum if desired
; wave_disc = sqrt(abs(wave_disc))
 
; construct scale and time arrays for the wavelet coefficients
di = N/2
; first two coefficients in wv_disc contain the 'residuums' S0 and S1
time_disc = [ 0 , di ] + di/2
scale_disc = [ di , di ]
; the following two coefficients belong to the largest scale
time_disc  = [ time_disc  , time_disc  ]
scale_disc = [ scale_disc , scale_disc ]
; the following coefficents belong to scales 2^i
while di gt 2 do begin
  di /= 2
  for j = 0L, n-1, di do begin
    time_disc = [ time_disc , di/2+j ]
  endfor ; j
  scale_disc = [ scale_disc , replicate(di,n/di) ]
endwhile
; convert to units of time by multiplying with timestep dt
time_disc = time_disc * dt + time[0]
scale_disc *= dt
 
; make a contour plot 
CONTOUR, $
  wave_disc[2:n-1] , $ ; drop the first two coefficients as they contain the residuums
  time_disc[2:n-1], scale_disc[2:n-1], $ ; reconstructed time and scale from above
   TITLE = 'discrete wavelet transform', $
   XTITLE = 'time', $
   YTITLE = 'period', $  
   /IRREGULAR, $ ; idl will triangulate the data before drawing in a regular grid
   /YLOG , $ ; log scaling on scale axis
   YRANGE = [ max(scale_disc) , min(scale_disc) ], $ ; invert direction of axis
   NLEVELS = 100, $ ; number of isoline levels
   /FILL, $  ; fill isolines with color
   XSTYLE = 1, $
   YSTYLE = 1

; continous wavelet transform
; wave_cont = WV_CWT(data, 'morlet',   8, /PAD, SCALE=scale_cont )
 
; power spectrum if desired
; wave_cont = sqrt(abs(wave_cont))
 
; convert to units of time by mutliplying with timestep dt
scale_cont *= dt
 
; Contour plot
CONTOUR, wave_cont, time, scale_cont, $
   TITLE = 'continous wavelet transform', $
   XTITLE = 'time', $
   YTITLE = 'period', $  
   /YLOG , $ ; log scaling on scale axis
   YRANGE = [ max(scale_cont) , min(scale_cont) ], $ ; invert direction of axis
   NLEVELS = 100, $ ; number of isoline levels
   /FILL, $  ; fill isolines with color
   XSTYLE = 1, $
   YSTYLE = 1

IDL can handle most of arithmetic errors without stopping.

You just get a message at the end of the program like:

% Program caused arithmetic error: Floating divide by 0

or:

% Program caused arithmetic error: Floating illegal operand

The latter might be something like assigning values outside the int-range to an array which has been defined as a=intarr(N).

To debug errors you can set the system variable

!EXCEPT = 2

this forces IDL to print an error message including the line number after each the line where an error occurs.
This should only used for debugging purposes - it slows down IDL by about 5%.

The default is !EXCEPT=1 meaning that error messages are printed at the end of the program. With !EXCEPT=0 no error messages are printed at all.

native idl routines and examples can be found in the lib directory of the idl distribution. IT can be found by investigating the !PATH variable. Currently it is /opt/itt/idl8.2/idl82/lib/

For plotting on the screen and for building pathes your code may need to know which operating system your code is running on. You can find this out by reading the !VERSION system variable:

on one of our linux machine it looks like:

IDL> help, !version , /stru
** Structure !VERSION, 8 tags, length=104, data length=100:
   ARCH            STRING    'x86_64'
   OS              STRING    'linux'
   OS_FAMILY       STRING    'unix'
   OS_NAME         STRING    'linux'
   RELEASE         STRING    '8.2.1'
   BUILD_DATE      STRING    'Aug 20 2012'
   MEMORY_BITS     INT             64
   FILE_OFFSET_BITS
                   INT             64

On a windows system it may look like

IDL> Help, !Version, /Structure
 ** Structure !VERSION, 8 tags, length=104, data length=100:
    ARCH            STRING    'x86_64'       OS              STRING    'Win32'
    OS_FAMILY       STRING    'Windows'
    OS_NAME         STRING    'Microsoft Windows'
    RELEASE         STRING    '8.0'
    BUILD_DATE      STRING    'Jun 17 2010'
    MEMORY_BITS     INT             64       FILE_OFFSET_BITS
                    INT             64 

For the correct device To plot on the screen (i.e. open a window to plot in ) you may use:

if !VERSION.OS_FAMILY eq 'Windows' then set_plot, 'Win' $
                                   else set_plot, 'X'

And similar for pathes:

if (!version.os_family eq 'linux')  then $
  srtm_dbase = '/home/jschween/exp/ideen/idl/srtm/srtm/' $
else $ 
  srtm_dbase = 'C:\Users\jschween\Downloads\srtm\'

Since November 2012 we have idl 8.2.1

Simply type idl and idl 8.2 will be executed. if no license for idl 8.2 is available 6.4 will be used. if you want to you use 6.4 explicitly, you need to type idl64. idlhelp and idlde will use the 8.2 versions. for 6.4 use idlhelp64 or idlde64, respectively.

To directly start version 6.4. type /opt/itt/idl6.4/idl64/bin/idl

To directly start version 8.2.1. type /opt/itt/idl6.4/idl64/bin/idl

If your code needs to know under which idl version it is running you may use the !VERSION system variable:

if !VERSION.RELEASE ge '8.2.1' then begin
  ...
endif else begin
  ...
endelse

To find out how many users are currently in IDL excute:

/opt/itt/idl8.2/idl82/bin/bin.linux.x86_64/lmstat -c /opt/itt/license/license.dat -f idl

in the listing should appear beside some other ten lines something like:

gop: license server UP (MASTER) v11.6
...
Users of idl:  (Total of 60 licenses issued;  Total of 18 licenses in use)  
...
 jschween qeb qeb/localhost:11 (v8.2) (gop/1700 522), start Wed 11/11 10:48, 6 licenses
...

The line starting with “Users of idl: …” says that there are 18 license in use. As one allways needs six licences when running idl this means that 18/6=3 (i.e. three) instances of idl are running. The line below indicates that user jschween on machine gop is running one instance of idl (using 6 licenses) v8.2 since Nov. 11.

If this is missing and/or you get a long strange output the license manager is probably not running.

To restart the license manager:

• log in to gop
• check whether it is running:

ps -elf | grep idl

• there should appear the licences manager daemon (lmgrd)
• if not:

sudo /etc/init.d/sys5_idl_lmgrd start

password is your password …

If licencses are obviously blocked by dead idl processes (running for several days, see section about license manager above). You may want to kill them (be careful not to destroy hard work of your colleagues !):

• Login on the respective machine
• see what idl processes are running by executing:

ps -elf | grep idl

• output may look like:

...
0 S hatpro   18936     1  0  80   0 -  3513 wait    2019 ?        00:00:00 /bin/bash /usr/local/bin/idl /home/dturner/aerioe/tmp/.run_aerioe.pro
0 R hatpro   18950 18936 98  80   0 - 85711 -       2019 ?        151-17:57:02 /opt/itt/idl8.2/idl82//bin/bin.linux.x86_64/idl /home/dturner/aerioe/tmp/.run_aerioe.pro
...

• indicating that user hatpro has an idl job ”run_aerioe.pro“ which was started in 2019 and is running for 151 days and ~18hours
• its process ID is 18950 and it has been started by process 18936
• eventually check with the (probable) owner whether this is intended
• to terminate this process do kill 18950 (this is the process ID in this case !)