Add IRF write methods

[thomasarmstrong]

Related to issue #1157. Added features to EffectiveAreaTable2D, EnergyDispersion2D and Background3D which enable a Bintable HDU to be constructed in the current CTA IRF format (or at least related to the format specified in https://github.com/open-gamma-ray-astro/gamma-astro-data-formats). It is rather simple at the moment and could be improved with versioning, as mentioned by @kosack, along with any other suggestions. Additionally I need to work on the PSF classes a little bit as the data format isn't quite compatible. just need to have theta_lo and theta_hi instead of just offset, and the possibility of writing out additional PSF dimensions? ( CTA IRF currently has the tables SCALE, SIGMA_1, AMPL_2, SIGMA_2, AMPLE_3, SIGMA_3). I will work on this in the next couple of days.

[cdeil]

@thomasarmstrong - How do you feel about re-exposing the axes as properties on the IRF class?

My suggestion would be to remove them, and instead to show in the class-level docstring an example how to obtain the axes via the data attribute (the NDDataArray). My reasoning is that it's better to start minimal, because it's better to have little code and adding is always possible in the future, but removing causes anger after people (or e.g. code in ctapipe) starts using it. I know at the moment some of our IRF classes have them re-exposed and some don't, so if you agree with this suggestion, please remove them from all IRF classes here.

If you think it's useful to re-expose them: please add a one-line docstring to each property, so that it shows up in a good way in the API docs. I.e. something like this:

Energy axis (`~gammapy.utils.nddate.BinnedDataAxis`).

Maybe @kosack also has an opinion here?

[thomasarmstrong]

@cdeil - so the idea would be to remove all the @propery definitions for energy, offset etc for all the IRF classes, and replace things like

energy = self.energy.bins

with

energy = self.axis('energy').bins

For example in EffectiveAreaTable2D.to_effective_area_table?

[cdeil]

Yes. But I don't have a strong preference here, so you can also go with what you prefer here.

[kosack]

that would be beneficial later if we have add more axes than just energy, offset, etc.

[cdeil]

I've had a quick look and left some inline comments.

This looks very good already and is a great addition!

But we first have to discuss where to use astropy.io.fits.BinTable or astropy.table.Table, and whether to use astropy.io.fits.Header or collections.OrderedDict for metadata (i.e. what Table uses).

There's pros and cons. My suggestion would be to use Table and OrderedDict for meta as much as possible, and to only go to fits.BinTable at the very last minute on FITS read / write (or explicitly via http://docs.astropy.org/en/stable/io/fits/api/tables.html#table-to-hdu).

This has the advantage of a cleaner and better code (Table is much better) and we can then also add e.g. HDF5 or other IRF serialisation in addition to FITS.
But for metadata it's a lossy conversion from astropy.io.fits.Header to collections.OrderedDict, one looses e.g. the comment for the header keys.

@thomasarmstrong @kosack - Do you have an opinion here or already a strategy / policy in CTApipe on this topic?

[cdeil]

This TestEnergyDispersion2D.test_write is failing:
https://travis-ci.org/gammapy/gammapy/jobs/302936895#L2666
Please have a look and let me know if it's not clear why.
Maybe it's about lower and upper case mismatch?

[cdeil]

There's formatting issues with your docstrings.

Locally I ran this:

make clean
python setup.py build_docs # very slow, 5 min
make doc-show

and see something like this:

Can you please fix up your docstrings?

I think in the examples section you have to do :: to declare a code block and indent:

Here's my code::

    print('thank you for this pull request!')

or you can use >>> without indentation if it's just a few lines of code or you want to mix prose and code:

Here's some code:

>>> print('spam!')

[thomasarmstrong]

I did not know that would slow it down so much! I will clean it up.

[cdeil]

This is an example of the code I'd like to avoid. astropy.table.Table computes the FITS header keywords for format and dim automatically. unfortunately, unit isn't handled automatically, but for now we can do it via a utility function (http://docs.gammapy.org/en/latest/api/gammapy.utils.fits.table_to_fits_table.html) and then later we can contribute this to the Table <-> BinTable converter in Astropy core.

[thomasarmstrong]

Table is indeed more simple, I think I had switched to BinTable just because I couldn't get Table to take the units for each column, but at least it goes from

c1 = fits.Column(name='ENERGY_LO', array=np.asarray([self.energy.lo]),
                 format='{}E'.format(self.energy.nbins), unit='{}'.format(self.energy.unit))
c2 = fits.Column(name='ENERGY_HI', array=np.asarray([self.energy.hi]),
                 format='{}E'.format(self.energy.nbins), unit='{}'.format(self.energy.unit))
c3 = fits.Column(name='OFFSET_LO', array=np.asarray([self.offset.lo]),
                 format='{}E'.format(self.offset.nbins), unit='{}'.format(self.offset.unit))
c4 = fits.Column(name='OFFSET_HI', array=np.asarray([self.offset.hi]),
                 format='{}E'.format(self.offset.nbins), unit='{}'.format(self.offset.unit))
c5 = fits.Column(name='EFFAREA', array=np.asarray([self.data.data.T]),
                 format='{}E'.format(self.energy.nbins * self.offset.nbins),
                 dim='({},{})'.format(self.energy.nbins, self.offset.nbins),
                 unit='{}'.format(self.data.data.unit))

to

table = Table([[self.energy.lo], [self.energy.hi], [self.offset.lo], [self.offset.hi], [self.data.data.T]],
              names=('ENERGY_LO', 'ENERGY_HI', 'OFFSET_LO', 'OFFSET_HI', 'EFFAREA'), meta=self.meta)

I believe table_to_fits_table already has the option to get the hdu name from the meta file, I suppose something similar to this could be done with the units for now?

[cdeil]

I think I would write it like this:

table = Table(meta=self.meta)

and then start filling one column at a time.
(Table is basically a dict of column name strings mapped to Column objects that hold the data, so filling column by column is efficient).

In [11]: table['ENERG_LO'] = aeff.data.axis('energy').lo

Assigning a Quantity object like that will create a column object with the unit filled from the quantity, which is what we want here:

In [12]: table.info()
<Table length=42>
  name    dtype  unit
-------- ------- ----
ENERG_LO float32  TeV

In [13]: table['ENERG_LO']
Out[13]: 
<Column name='ENERG_LO' dtype='float32' unit='TeV' length=42>
0.0125893
0.0158489
0.0199526
0.0251189
0.0316228
0.0398107
0.0501187
0.0630957
0.0794328
      0.1
 0.125893
 0.158489
 0.199526
 0.251189
 0.316228
 0.398107
 0.501187
      ...
  5.01187
  6.30957
  7.94328
     10.0
  12.5893
  15.8489
  19.9526
  25.1189
  31.6228
  39.8107
  50.1187
  63.0957
  79.4328
    100.0
  125.893
  158.489

In [14]: type(table['ENERG_LO'])
Out[14]: astropy.table.column.Column

So to_table would return a Table.
And then read and write would call read and write on that Table object, or if that looses the units or other metadata, we should go via gammapy.utils.fits.table_to_fits_table and gammapy.utils.fits.fits_table_to_table (and if that still doesn't work properly, fix up those utility functions, but only do that work of proper Table to BinTable conversion once in a central place in those helper functions. If it's unclear or you run into any issues (not unheard of with FITS), comment here and I'll take a closer look.

[cdeil]

Additionally I need to work on the PSF classes a little bit as the data format isn't quite compatible.

I would suggest to do that in a separate PR and keep this one limited to AEFF, EDISP and BACKGROUND.

Note that http://docs.gammapy.org/en/latest/api/gammapy.irf.PSF3D.html (the PSF class that's most important and probably the one to work on first) hasn't been changed to use NDDataArray yet, so doing that change will likely lead to quite a bit of work of fixing up callers and high-level anlaysis code, i.e. something best done separate from a PR with the main focus on FITS I/O.

[cdeil]

@thomasarmstrong - I've left one last round of inline comments. It's quite a few changes, but I hope they should all relatively straightforward to address. Let me know if you have any questions or disagree with something!

Finally, there should be at least a 1-line docstring for the to_table methods, like

def to_table(self):
    """Convert to `~astropy.table.Table`."""

And how about adding one-line methods like this?

def to_fits(self):
   """Convert to `~astropy.io.fits.BinTable`."""
    return table_to_fits_table(self.to_table())

I know its a bit of boilerplate, but it's convenient to have from calling code doing FITS I/O like in your example_write_irf.py, while still letting us keep the "clean" serialisation format to_table and from_table which could be the basis for other serialisation formats like already directly to HDF5, which Astropy Table supports.

[cdeil]

Please follow the spec (see here, i.e. use ENERG_LO and ENERG_HI.

We can discuss and make changes to the formats of course, or first introduce / prototype new formats in Gammapy / ctapipe and then write a spec for them after. But let's add support for the spec we have and have been using for current IACTs for a few years first.

[thomasarmstrong]

Will do (same as for the background), I am not sure why CTA has gone for slightly different definitions in their IRFs?

[cdeil]

I am not sure why CTA has gone for slightly different definitions in their IRFs?

My 2 cents: CTA hasn't committed to any formats. What we have is preliminary and variants are due to missing communication in the past. Concerning the specific issue, the discussion for CTA should continue here: https://forge.in2p3.fr/boards/236/topics/1824?r=2314#message-2314 gamma-astro-data-formats.readthedocs.io is just an attempt to write up the formats that we use, to have a format description at all and to avoid duplicating it (in this case in the ctapipe, Gammapy and Gammalib docs).

There is a lot of discussion and complaining, but I think basically we just need more communication and collaboration between pipeline and science tool developers and things can be OK from here on out for CTA.

[cdeil]

I see, you were commenting on reco vs true energy axis names.
This has been discussed previously, e.g. here: open-gamma-ray-astro/gamma-astro-data-formats#63
I don't remember what I or others said, in principle I don't care much about FITS serialisation details. The problem with changing something is that it needs an active communication channel and collaboration, and as Karl suggested versioning. I hope that this will come soon, but for now just sticking with what we have in the spec (and used e.g. in HESS for years) seems simplest.

[cdeil]

Please follow the spec (see here, i.e. use BKG.

[cdeil]

I think it would be best if the "to_table" method didn't modify self. It's supposed to be a serialisation, so the self object should remain unchanged IMO.

Would this work?

meta = self.meta.copy()
table = Table(meta=meta)

Is adding the name key really needed or useful? If yes, you could add it to this copy of the meta dict here, or move it to some other place, e.g. __init__.

This comment applies to all to_table methods.

[thomasarmstrong]

It could be taken out all together, but the user would have to make sure to include it when adding the provenance data. The function table_to_fits_table extracts this in order to assign the name to the fits table. An alternative/additional solution could be:

def table_to_fits_table(table, name=None):
    ...
    if name is None:
         name = table.meta.pop('name', None)

which could allow this to be placed earlier (i.e. to_fits could have the same format). Either way, should I still include the

meta = self.meta.copy()
table = Table(meta=meta)

so that when additional information is included (such as in the example) it doesn't modify self?

[cdeil]

Hmm ... OK, now I see that handling of name is difficult because:

• BinTableHDU has name separate from header (but presumably they have some auto-interaction e.g. on read / write with header keywords).
• Table doesn't have a name attribute, so the urge to combine header and name into meta is strong, but then the question is what all the places are where to do this.

I see that in Table.read and Table.write the HDU name round-trips through the EXTNAME key in meta:

In [1]: from astropy.table import Table

In [2]: t = Table.read('http://gamma-astro-data-formats.readthedocs.io/en/latest/_downloads/aeff_2d_full_example.fits')
Downloading http://gamma-astro-data-formats.readthedocs.io/en/latest/_downloads/aeff_2d_full_example.fits
|=======================================================================================================================|  14k/ 14k (100.00%)         0s
WARNING: hdu= was not specified but multiple tables are present, reading in first available table (hdu=1) [astropy.io.fits.connect]

In [3]: t.meta
Out[3]: 
OrderedDict([('OBS_ID', 31415),
             ('LO_THRES', 0.1),
             ('HI_THRES', 50),
             ('HDUDOC',
              'https://github.com/open-gamma-ray-astro/gamma-astro-data-formats'),
             ('HDUVERS', '0.2'),
             ('HDUCLASS', 'GADF'),
             ('HDUCLAS1', 'RESPONSE'),
             ('HDUCLAS2', 'EFF_AREA'),
             ('HDUCLAS3', 'FULL-ENCLOSURE'),
             ('HDUCLAS4', 'AEFF_2D'),
             ('EXTNAME', 'EFFECTIVE AREA')])

In [4]: t.write('/tmp/test.fits')
WARNING: UnitsWarning: 'None' did not parse as fits unit: At col 0, Unit 'None' not supported by the FITS standard.  [astropy.units.core]

In [5]: from astropy.io import fits

In [6]: hdu_list = fits.open('/tmp/test.fits')

In [7]: hdu_list.info()
Filename: /tmp/test.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU       4   ()      
  1  EFFECTIVE AREA    1 BinTableHDU     35   1R x 5C   [20D, 20D, 5D, 5D, 100D]   

But as far as I know they have no support for writing multiple HDUs to file, this would clobber, but by default fails if the file exists:

t.meta['EXTNAME'] = 'AEFF V2'   
t.write('/tmp/test.fits') 

So basically I think sane FITS extension name handling is an unsolved issue in Gammapy (and probably ctapipe)?
@thomasarmstrong - I think maybe what you suggest could be the best way to proceed: add a name option for the functions and methods that convert to BinTableHDU. But probably it should default to the key EXTNAME in meta, not name like it is now?

It's a bit perverse, but this is the hardest part of this PR. So if you can add tests that actually write to FITS files to establish the current behaviour and show that it works, i.e. a caller can have control over HDU names that get written, that would be useful.

[cdeil]

copy & paste error on the line filling table['DETX_LO']: should be self.data.axis('detx').lo, you have self.data.axis('detx').hi

[cdeil]

You currently have this code to insert a new length 1 dimension, to fit the array into a single FITS table row as in the format that we use:

table['ENERG_LO'] = [self.data.axis('energy').lo]*self.data.axis('energy').unit

I would suggest to write it in this style instead everywhere:

table['ENERG_LO'] = self.data.axis('energy').lo[np.newaxis]

The result is the same, I just find the code a little shorter and simpler.

In the from_table methods I think we should use this to again drop this "fit into single row" dimension that is a curiosity of the FITS BINTABLE serialisation format:

energy_lo = table['ENERG_LO'].quantity[0]

Then the code is pretty symmetrical and easy to read.

@thomasarmstrong - Does that sound OK to you?

If yes, would you mind also fixing up these lines to use that coding pattern?
https://gist.github.com/cdeil/40c787d6188bee1c2f214246a84d3ca0
(previously we used squeeze to remove the extra first length 1 dimension. But that will actually also squeeze out other axes if they have length 1, which can happen, so really that is a bug waiting to happen)

[ZiweiOu]

Dear Doctor Deil, My name is Ziwei Ou. I am a PhD student from IPN Orsay, working for CTA simulations. One of my main propose to study GammaPy is plotting IRFs of CTA. Recently I check the Python script from GammaPy call "example_fermi_psf" in the examples file. There is a problem shows " TypeError: kernel() got an unexpected keyword argument 'pixel_size' " for " kernel = psf.kernel(pixel_size=pixel_size, offset_max=offset_max) " Do you have any suggestion or advice for that? Thank you very much. Best wishes, Ziwei

…

On Tue, Nov 21, 2017 at 9:47 AM, Christoph Deil ***@***.***> wrote: ***@***.**** requested changes on this pull request. @thomasarmstrong <https://github.com/thomasarmstrong> - I've left one last round of inline comments. It's quite a few changes, but I hope they should all relatively straightforward to address. Let me know if you have any questions or disagree with something! Finally, there should be at least a 1-line docstring for the to_table methods, like def to_table(self): """Convert to `~astropy.table.Table`.""" And how about adding one-line methods like this? def to_fits(self): """Convert to `~astropy.io.fits.BinTable`.""" return table_to_fits_table(self.to_table()) I know its a bit of boilerplate, but it's convenient to have from calling code doing FITS I/O like in your example_write_irf.py, while still letting us keep the "clean" serialisation format to_table and from_table which could be the basis for other serialisation formats like already directly to HDF5, which Astropy Table supports. ------------------------------ In gammapy/irf/energy_dispersion.py <#1211 (comment)>: > @@ -924,3 +910,17 @@ def peek(self, figsize=(15, 5)): edisp.plot_matrix(ax=axes[2]) plt.tight_layout() + + def to_table(self): + #TODO : ETRUE or ENERG? See https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/edisp/index.html Please follow the spec (see here <https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/edisp/index.html>, i.e. use ENERG_LO and ENERG_HI. We can discuss and make changes to the formats of course, or first introduce / prototype new formats in Gammapy / ctapipe and then write a spec for them after. But let's add support for the spec we have and have been using for current IACTs for a few years first. ------------------------------ In gammapy/irf/background.py <#1211 (comment)>: > @@ -114,3 +115,17 @@ def read(cls, filename, hdu='BACKGROUND'): filename = make_path(filename) hdulist = fits.open(str(filename)) return cls.from_hdulist(hdulist, hdu=hdu) + + def to_table(self): + #TODO: BGD or BKG - see https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/bkg/index.html Please follow the spec (see here <https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/bkg/index.html>, i.e. use BKG. ------------------------------ In gammapy/irf/background.py <#1211 (comment)>: > @@ -114,3 +115,17 @@ def read(cls, filename, hdu='BACKGROUND'): filename = make_path(filename) hdulist = fits.open(str(filename)) return cls.from_hdulist(hdulist, hdu=hdu) + + def to_table(self): + #TODO: BGD or BKG - see https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/bkg/index.html + table = Table() + table['DETX_LO'] = [self.data.axis('detx').hi]*self.data.axis('detx').unit + table['DETX_HI'] = [self.data.axis('detx').hi]*self.data.axis('detx').unit + table['DETY_LO'] = [self.data.axis('dety').lo]*self.data.axis('dety').unit + table['DETY_HI'] = [self.data.axis('dety').hi]*self.data.axis('dety').unit + table['ENERG_LO'] = [self.data.axis('energy').lo]*self.data.axis('energy').unit + table['ENERG_HI'] = [self.data.axis('energy').hi]*self.data.axis('energy').unit + table['BGD'] = [self.data.data]*self.data.data.unit + self.meta.update({'name':'BACKGROUND'}) I think it would be best if the "to_table" method didn't modify self. It's supposed to be a serialisation, so the self object should remain unchanged IMO. Would this work? meta = self.meta.copy() table = Table(meta=meta) Is adding the name key really needed or useful? If yes, you could add it to this copy of the meta dict here, or move it to some other place, e.g. __init__. This comment applies to all to_table methods. ------------------------------ In gammapy/irf/background.py <#1211 (comment)>: > @@ -114,3 +115,17 @@ def read(cls, filename, hdu='BACKGROUND'): filename = make_path(filename) hdulist = fits.open(str(filename)) return cls.from_hdulist(hdulist, hdu=hdu) + + def to_table(self): + #TODO: BGD or BKG - see https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/bkg/index.html + table = Table() + table['DETX_LO'] = [self.data.axis('detx').hi]*self.data.axis('detx').unit copy & paste error on the line filling table['DETX_LO']: should be self.data.axis('detx').lo, you have self.data.axis('detx').hi ------------------------------ In gammapy/irf/effective_area.py <#1211 (comment)>: > @@ -570,3 +563,15 @@ def peek(self, figsize=(15, 5)): self.plot_energy_dependence(ax=axes[0]) self.plot_offset_dependence(ax=axes[1]) plt.tight_layout() + + def to_table(self): + table = Table() + table['ENERG_LO'] = [self.data.axis('energy').lo]*self.data.axis('energy').unit You currently have this code to insert a new length 1 dimension, to fit the array into a single FITS table row as in the format that we use: table['ENERG_LO'] = [self.data.axis('energy').lo]*self.data.axis('energy').unit I would suggest to write it in this style instead everywhere: table['ENERG_LO'] = self.data.axis('energy').lo[np.newaxis] The result is the same, I just find the code a little shorter and simpler. In the from_table methods I think we should use this to again drop this "fit into single row" dimension that is a curiosity of the FITS BINTABLE serialisation format: energy_lo = table['ENERG_LO'].quantity[0] Then the code is pretty symmetrical and easy to read. @thomasarmstrong <https://github.com/thomasarmstrong> - Does that sound OK to you? If yes, would you mind also fixing up these lines to use that coding pattern? https://gist.github.com/cdeil/40c787d6188bee1c2f214246a84d3ca0 (previously we used squeeze to remove the extra first length 1 dimension. But that will actually also squeeze out other axes if they have length 1, which can happen, so really that is a bug waiting to happen) — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#1211 (review)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/Ad9AcZ6pA6mPQPjkqD1UoLD4PNEF4x2Kks5s4o4MgaJpZM4QgOj5> .

[cdeil]

@ZiweiOu - Please call me Christoph (or @cdeil on Github). Also - please always start a new thread on the Gammapy mailing list (https://groups.google.com/forum/#!forum/gammapy) or the Gammapy issue tracker (https://github.com/gammapy/gammapy/issues/new), unless a given thread or issue is exactly on point or very much related to your issue.

In this case, it's not related to the work that @thomasarmstrong is doing here, so I've gone ahead and filed a separate issue #1216 concerning PSF kernels. @ZiweiOu - You can subscribe to that issue if you want to follow the discussion on the PSF kernel (see https://help.github.com/articles/subscribing-to-and-unsubscribing-from-notifications/#subscribing-to-issues-and-pull-requests).

[cdeil]

@thomasarmstrong - You have this:

gamma = Quantity(table['GAMMA'].quantity[0], table['GAMMA'].unit)

This should be the same, just passing through the Quantity constructor one times less:

gamma = table['GAMMA'].quantity[0]

Can you please change here and review the diff and also do it in the other cases, at least where it's clearly the same. (changing e.g. from Angle to Quantity might break callers, altough for uniformity we probably should change to Quantity everywhere. (Angle is a Quantity subclass with a little sugar.

[cdeil]

@thomasarmstrong - I think the EnergyDispersion2D class should be changed to remove this transpose in the from_table and to_table, and probably then you also need to update the order of axes as listed in __init__ and then it should work. That would change the numpy array axis order in memory, but code everywhere in Gammapy should evaluate axes by name and all calling code should still work.

That would still be a temp workaround, I think the final solution would be to parse the axis order from the CREF header key. That would make Gammapy independent of the axis order in the FITS file, we could then read any axis order. For generating new FITS files, a class-level attribute that specifies the "preferred" or "default" axis order would still be needed, i.e. for code that creates IRFs from scratch like you want to in ctapipe.

Previous discussion on IRF axis order is at open-gamma-ray-astro/gamma-astro-data-formats#28 and probably other issues.

My suggestion would be that you try the small change I mentioned and see if python setup.py test -V still passes, but leave the larger change involving CREF to a future PR. (or directly postpone this to a future PR).

[cdeil]

@thomasarmstrong - Why do you call .all() here? I think what that does is to convert the data arrays to bool and check if all entries are 0 or non-0. This is not a useful assert, no? Did you mean this?

assert hdu.data['DETX_LO'][0] == bkg_3d.data.axis('detx').lo.value

I think that would be a useful assert and work as long as the array dtypes are the same, i.e. no float64 to float32 conversions occur. I.e. bit-by-bit comparison should work here, but if it doesn't, then you should call numpy.testing.assert_allclose.

[thomasarmstrong]

There seems to be some problems elsewhere in gammapy from changing the data axis, I will work on this and add a commit shortly.

[cdeil]

@thomasarmstrong - Thanks! If it's confusing, I can also take a look and fix up other parts of the Gammapy code or tests by adding a commit here. Let me know.

[thomasarmstrong]

@cdeil - I am starting to wonder if it might be easier to leave the transpose of the effective area and energy dispersion data arrays... There are many places throughout the code that call for the data/axis to be in the original format, which seems like a lot of work to do for little gain? I have little opinion on what is better so I defer to your insight.

[cdeil]

@thomasarmstrong - I would suggest to defer the EDISP transpose change to later.

Otherwise - is this PR ready from your side?

[thomasarmstrong]

@cdeil - I have reverted some of the changes and it should be ready now. I guess the next step is to open a new issue/PR for the PSF classes to use NDDataArray and to be compatible with the current CTA IRFs?

[cdeil]

@thomasarmstrong - Thanks! I don't have time today, but I will merge this in tomorrow.

I think the next steps are:

• Improve code and tests for gammapy.irf.PSF3D, especially use NDDataArray to store the IRF and make FITS I/O work.
• Improve code and tests for gammapy.irf.EnergyDependentMultiGaussPSF. There I think we also want to use NDDataArray, but currently I think we have no case of analytical IRFs using it, and it will need some thought what the interpolation behaviour should be. I think at the moment it is interpolating the model parameters of the analytical model, and that's acceptable behaviour and should be possible to go to NDDataArray and FITS I/O like for the other IRF classes.
• I also think a reminder issue that the data array transpose on IRF FITS read / write should be looked at and probably best to clean up Gammapy code to avoid it. (not sure, needs some thought / discussion)
• And more generally we might or not might want to change NDDataArray to support arbitrary axis orders, and in FITS I/O use the CREF header key to determine axis order. (not sure, needs some thought / discussion)
• We probably also want some code to smooth and fits IRF histograms to derive smooth or anlytical IRFs. For that I'm not sure where it should live. It could partly be methods on the IRF classes, or separate functions or classes in gammapy.irf or in ctapipe.
• Finally, we have several old very preliminary formats and classes for IRF code, that comes from first attempts to export to FITS in HESS from a few years ago. There I have the overview and will clean that up, i.e. mostly remove them, in the coming weeks.

Of course pull requests (probably one per class is best, keep them small) are welcome any time, but it's also clear that this is a big task for the coming months and we should try to discuss and distribute the work a bit. From the Gammapy side, we now have bi-weekly calls on eZuce on Fridays at 11 am, the next ones are Dec 1 and 15. We are also starting to plan our next coding sprint, which will be early next year in Paris. If you're interested, please fill https://doodle.com/poll/2q8rw8n8edpm6fn5 and for the weekly calls, let @bkhelifi or me know if you have time to join and talk about IRFs.