View source for Two-dimensional histograms (xfig example)



[[Category:SLxfig]]

<p style="font-size:15px"><b>Plotting two-dimensional histograms using the png and xfig modules</b></p style="font-size:15px"><br/>
<i>example added by Manfred</i>

Similarly to the [[Example: color-coded maps with 2d histograms, png, and xfig|first example]] for creating [[Example: color-coded maps with 2d histograms, png, and xfig|color-coded maps]] with the [[Plotting_with_the_S-Lang_Xfig_module|SLxfig module]], the goal of this one is to produce a color-coded map,
but this time, the color-scale represents a density of (a possibly large number of) data points.

=== ISIS' histogram2d function ===

==== Usage ====
The <code>histogram2d</code> function is used in the following way (see also [http://pulsar.sternwarte.uni-erlangen.de/hanke/science/X-ray/ISIS/cgi/help.cgi?subj=histogram2d|help("histogram2d");]):
  num = histogram2d(A, B, Agrid, Bgrid);
''A'' and ''B'' have to be ''Double_Type'' arrays of the same length,
''Agrid'' and ''Bgrid'' have to be ''Double_Type'' arrays with monotonically increasing elements.

Now, <code>num[a,b]</code> contains the <strong>number of points</strong> (<code>A</code>, <code>B</code>)
falling into each 2d-bin defined by the lower bin boundaries <code>Agrid[a]</code> and <code>Bgrid[b]</code>, i.e.:
  num[a,b]  ==  howmany( (Agrid[a] <= A < Agrid[a+1]) and (Bgrid[b] <= B < Bgrid[b+1]) );

==== Overflow bins ====
The last bin in each row or column is an //**overflow bin**// whose upper limit is at infinity, i.e.,
  num[-1,b]  ==  howmany( (Agrid[-1] <= A) and (Bgrid[b] <= B < Bgrid[b+1]) );
  num[a,-1]  ==  howmany( (Agrid[a] <= A < Agrid[a+1]) and (Bgrid[-1] <= B) );
and
  num[-1,-1]  ==  howmany( (Agrid[-1] <= A) and (Bgrid[-1] <= B) );

If one wants to remove these overflow bins,
one can simply apply the index-array ''[:-2]'' (which means "take all elements up to the second last one")
in both dimensions in order to create a smaller two-dimensional array:
  num_ = num[ [:-2], [:-2] ];

The last bins corresponding to the reduced 2d-array ''num_'' //end// now at the same position where the grids //end//, too:
  num_[-1,-1]  ==  num[-2,-2]  ==  howmany( (Agrid[-2] <= A < Agrid[-1]) and (Bgrid[-2] <= B < Bgrid[-1]) ); 
Before, the grids specified always lower bin boundaries,
but now, the grids have one element more than the correpsonding array-dimensions of ''num_''.
  
Note that ISIS (one-dimensional) ''histogram'' function does not produce an overflow bin.

==== Order of the indices ====
The first index of ''num'' (called ''a'' above) will always correspond
to the first [and third] argument of the ''histogram2d'' function (here called ''A'' [and ''Agrid'']),
and the second index (''b'') to the second [and fourth] argument (''B'' [and ''Bgrid'']).

As most functions that operate on two dimensional arrays (like ''plot_image'', ''png_write'', ''ds9_view'', ...)
expect that the //first index// of the 2d-array to correspond to the y-axis,
whereas only the //second index// to the x-axis,
one will have to pass the y-data [and -grid] before the x-data [and -grid]:
  num = histogram2d(Y, X, Xgrid, Ygrid);


=== Using two-dimensional histograms ===

Just as in the [[isis:slxfig:ccmap|previous example]],
these two dimensional histograms can be transformed into a color-coded image,
which can then be included in an [[isis:slxfig|xfig plot]] which adds back the original coordinate system.

==== A "banana-plot" ====

In an X-ray gratings spectrometer (//Chandra//-HETGS, //XMM//-RGS),
the dispersion is not only proportional to the photon's wavelength, but also to the dispersion order.
Different orders, however, can be separated with the internal energy resolution of the CCD detector.

The following example plots the density of events in a position--energy diagram:

===== Writing the png image =====

<pre>
variable evts = fits_read_table("evt2.fits");          % reading a Chandra level-2 event file

variable x_grid = [2300:5700:#320];                    % x position on the detector (320 values)
variable energy_grid = [0:10000:#240];                 % energies in eV (240 values from 0 to 10 keV)

variable num = histogram2d(evts.energy, evts.x,        % the data (Y, X)
                           energy_grid, x_grid         % the grids (Ygrid, Xgrid)
                          )[[:-2],[:-2]];              % removing overflow bins

require("png");                                        % load png module
variable RGB = png_gray_to_rgb(log(1+num), "ds9sls");  % transform num values logarithmically to RGB values, using the ds9sls color map
png_write_flipped("banana.png", RGB);                  % write RGB values to a png file
                                                       % png_write_flipped places RGB[0,0] in the *lower* left corner
</pre>
This produces just the bare image "banana.png":

{{:isis:slxfig:banana.png|bare image "banana.png"}}

With the [[isis:slxfig|S-Lang xfig module]], the coordinate system can, however, be easily added (back again).

===== Plotting the png into an SLxfig plot =====

It is important to set the correct coordinates for the SLxfig plot object,
as this information is no more present in the bare image.
As discussed above (in the section on "overflow bins"),
the last bins of the reduced 2d-array end where the original grids end
-- and obviously the first bins start where the grids start.
This is all the information that needs to be passed to the ''world'' method of an SLxfig plot object:
<pre>
require("xfig");
variable p = xfig_plot_new();
p.world(x_grid[0], x_grid[-1], energy_grid[0]/1000, energy_grid[-1]/1000);  % defining the coordinate system; energy in keV instead of eV
p.plot_png("banana.png");
p.xlabel("Sky X Position");
p.ylabel("Energy [keV]");
p.render("banana-plot.pdf");
</pre>
The resulting plot looks like this:

{{:isis:slxfig:banana-plot.png|banana-plot with coordinate system}}