Improvements to reading / writing FITS images #2246

jzuhone · 2019-04-17T15:36:57Z

This PR introduces improvements to FITS file writing and reading. Specifically:

Specifies a new standard for writing metadata to FITS files using the FITSImageData class and its subclasses. This standard uses the header to write unit information and other metadata to the FITS file so that when it is read back into yt as a dataset it is more self-describing.
Defines a new class, YTFITSDataset, which is a subclass of FITSDataset which looks specifically for this header metadata.
Simplifies the logic in FITSHierarchy for determining units.
Makes the "beam" unit part of SkyDataFITSDataset.
Adds a convolve method to FITSImageData instances to convolve them with any kernel created using AstroPy or a simple Gaussian kernel.
Specifies that the default image resolution for FITSImageData subclasses should be 512. The previous behavior was inconsistent across the various subclasses.

The new yt-based standard for writing FITS datasets with self-describing information and reading them back in as datasets is set out in a YTEP:

yt-project/ytep#5

Things left to do:

add built-in method to YTCoveringGrid to write these grids to FITS files
change docs
fix / add more tests

jzuhone · 2019-04-23T01:38:44Z

Assuming the tests pass, this is now finally ready to go. I know it's huge and that it's not a very familiar part of the code, for which I apologize.

… them. Move the beam stuff to the sky fits dataset because that's what it is for.

…t exists) to the file headers

…it parsers which should in theory be more robust

…ch means that all of its derived arrays and quantities will also be. Some other minor stylistic nits.

…eader in the FITSDataset class

…t spectrum as 3rd dimension

matthewturk

lgtm!

matthewturk · 2019-05-14T17:25:33Z

doc/source/examining/loading_data.rst

 * ``BUNIT``: The units of the image

 FITS header keywords can easily be updated using AstroPy. For example,
 to set the ``BTYPE`` and ``BUNIT`` keywords:

 .. code-block:: python

-   import astropy.io.fits as pyfits
-   f = pyfits.open("xray_flux_image.fits", mode="update")
+   from astropy.io import fits


how modern! 😁

matthewturk · 2019-05-14T17:27:07Z

yt/frontends/fits/data_structures.py

@@ -69,11 +67,13 @@
 sky_prefixes = list(sky_prefixes)
 spec_prefixes = list(spec_names.keys())

-field_from_unit = {"Jy":"intensity",
-                   "K":"temperature"}
+field_from_unit = {"Jy": "intensity",


Do we ever expect these to be imported outside of this file?

Nope, I'll fix that.

matthewturk · 2019-05-14T17:27:21Z

yt/frontends/fits/data_structures.py

@@ -92,7 +93,7 @@ def __init__(self,ds,dataset_type='fits'):
        self.dataset_type = dataset_type
        self.field_indexes = {}
        self.dataset = weakref.proxy(ds)
-        # for now, the index file is the dataset!
+        # for now, the index file is the dataset


:inception_noise:

matthewturk · 2019-05-14T17:27:56Z

yt/frontends/fits/data_structures.py

-            if field_units[0] == "/":
-                field_units = "1%s" % field_units
-            return field_units
+            try:


Is there any way we could flatten this?

If possible, flattening this would make it a bit more readable!

matthewturk · 2019-05-14T17:28:40Z

yt/frontends/fits/data_structures.py

@@ -410,21 +402,22 @@ def _set_code_unit_attributes(self):
            length_factor = self.wcs.wcs.cdelt[0]
            length_unit = str(file_units[0])
            mylog.info("Found length units of %s." % (length_unit))
-        else:
+        elif "length_unit" not in self.units_override:


Any chance we'll need a final else block?

I'll add one. Thanks for catching this.

munkm

Overall, this looks great! I noticed the YTEP still has the old PR (not the one you rebased with 4.0) linked to it for management. Maybe include this PR in that too?

munkm · 2019-05-15T22:08:44Z

doc/source/visualizing/FITSImageData.ipynb

   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
-    "A 3D FITS cube can also be created from a covering grid:"
+    "A 3D FITS cube can also be created from regularly gridded 3D data. In yt, there are covering grids and \"arbitrary grids\". The easiest way to make a grid is using `ds.r`:"


I think that the transition from mentioning yt's covering grids and arbitrary grids, and then saying that a user should make a grid with ds.r() is a little bit confusing. Maybe a transition sentence that the ds.r() is creating an ArbitraryGrid object based on the boundaries you're passing it would help?

Sounds good

munkm · 2019-05-15T22:19:18Z

doc/source/examining/loading_data.rst

@@ -992,7 +989,7 @@ zero, and the left and right edges of the domain along this axis are
 Examples of Using FITS Data
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^

-The following IPython notebooks show examples of working with FITS data in yt,
+The following Jupyter notebooks show examples of working with FITS data in yt,
 which we recommend you look at in the following order:

 * :ref:`radio_cubes`


Do you also want to add a link to the FITSImageData notebook in /visualizing/ here?

Sounds good

munkm · 2019-05-15T22:27:53Z

doc/source/examining/loading_data.rst

-   import astropy.io.fits as pyfits
-   f = pyfits.open("xray_flux_image.fits", mode="update")
+   from astropy.io import fits
+   f = fits.open("xray_flux_image.fits", mode="update")
   f[0].header["BUNIT"] = "cts/s/pixel"


In the notebook you showed updating header values with the update_header() method. Can that be used to do the same thing as the code block here too?

This is actually a different situation where we show how to make FITS files more self-descripting so yt can read them properly.

Ah ok! That makes sense! Is there any threshold of requirements that yt needs to read FITS files properly?

Not really a threshold, it's more that the more self-descripting your data is, the better it is supported. Even data that's not very self-descripting can still be loaded.

munkm · 2019-05-15T22:32:19Z

doc/source/examining/loading_data.rst

@@ -833,7 +847,7 @@ The third way is if auxiliary files are included along with the main file, like

 .. code-block:: python

-   ds = load("flux.fits", auxiliary_files=["temp.fits","metal.fits"])
+   ds = yt.load("flux.fits", auxiliary_files=["temp.fits","metal.fits"])


Just checking: this is an example of a case where there might be auxiliary files. We don't have any datasets in the FITS frontend that have these, right?

@munkm we don't, no. I don't think I have any on hand at the moment, actually.

That's ok! I just wanted to check that I shouldn't be able to run it with any files that I currently have locally.

munkm · 2019-05-15T22:47:55Z

yt/visualization/tests/test_fits_image.py

@@ -87,10 +91,12 @@ def test_fits_image():
    cut = ds.cutting([0.1, 0.2, -0.9], [0.5, 0.42, 0.6])
    cut_frb = cut.to_frb((0.5, "unitary"), 128)

-    fid3 = FITSImageData(cut_frb, fields=[("gas","density"), ds.fields.gas.temperature], units="cm")
+    fid3 = cut_frb.to_fits_data(fields=[("gas","density"),


Maybe change fields=[("gas","density"), to add a space between the fields, like fields=[("gas", "density"), here and to the fields in the fits_cut three lines down?

brittonsmith

One totally minor comment, respond or ignore as you like. I'm happy for this to be merged.

brittonsmith · 2019-05-21T10:02:49Z

yt/frontends/fits/fields.py

+
+        # We need to have a bunch of species fields here, too
+        ("metal_density",   ("code_mass/code_length**3", ["metal_density"], None)),
+        ("hi_density",      ("code_mass/code_length**3", ["hi_density"], None)),


Do you want to also add the yt standard alias names (as per YTEP 3) for these species fields, so for example, 'H_p0_density`?

I talked to @jzuhone about this and he's going to address it in a different PR.

yt/visualization/fits_image.py

jzuhone force-pushed the fits branch 2 times, most recently from a0d51dd to a475b1c Compare April 22, 2019 19:01

jzuhone changed the title ~~[WIP] Improvements to reading / writing FITS images~~ Improvements to reading / writing FITS images Apr 23, 2019

jzuhone force-pushed the fits branch from 9592d59 to 2c04416 Compare April 25, 2019 02:30

jzuhone and others added 25 commits May 10, 2019 14:46

Allow the user to set the image center

d399c23

shorter

a7f3d4d

Separate out the yt-generated FITS files into a new class which reads…

4a82a53

… them. Move the beam stuff to the sky fits dataset because that's what it is for.

Refactor of this to write all of the units and the current time (if i…

8685e5d

…t exists) to the file headers

Add an attribute for the header's comments

6efba9a

Loop over and get the units from the header

68296dc

Get the current_time from the header if it's there

286703d

Special handling for velocity and magnetic units

c38c18f

A method to fix the input current time

be76386

Bugfix

cae102d

I obviously didn't understand how that worked

192789c

Bugfix

cf24b34

Bugfixes

5c071b9

Try to determine the units of the image using AstroPy and yt-based un…

af20d62

…it parsers which should in theory be more robust

Bugfix

4501dda

Oops

b603797

More refining on this

f551ed8

Removing stuff we don't need anymore

31dd2ff

Bugfix

6192ea7

Make the left and right edges in the length units of the dataset, whi…

7ba6d64

…ch means that all of its derived arrays and quantities will also be. Some other minor stylistic nits.

Bugfixes

5677dac

Fix comment

63e0ce5

Drop the fits suffix from the repr

c5f8a03

Fix bad header keywords

4924370

Fix magnetic field units

9b9f946

jzuhone added 22 commits May 10, 2019 14:46

Restore some old code to allow length units to be detected from the h…

a8f3284

…eader in the FITSDataset class

This fixes a regression--we want to include cubes with linear axes bu…

e07f132

…t spectrum as 3rd dimension

Make sure we can handle code units

143dbc7

Fix deprecation warnings and whitespace

65cb0fc

A function to change the name of a fits image

93810af

Some doc fixes

fa5e244

Some doc work

52bd6e1

Updating docs for FITS data

b4b7b00

Forgot this

5b71a06

bugfix

171e10b

bugfix

fef67ba

Make sure we can only do one image, and make sure that we make a copy

90672db

Fix bugs

60a6c50

Pass extra arguments to convolve

339d1ed

This only works for 2D data

6e67659

Update FITS tests

dc03347

Update docs

d452e45

These should be set by the dataset

5df175f

Failsafe against using code units here

3a2c01e

Fix corner case of length units which are not pure

b3ffe21

account for potentially duplicate names

e614199

Don't fail here but convert to cgs instead

6b98e5d

jzuhone force-pushed the fits branch from 573d53e to 6b98e5d Compare May 10, 2019 18:46

matthewturk approved these changes May 14, 2019

View reviewed changes

munkm approved these changes May 15, 2019

View reviewed changes

Responding to comments by matthewturk and munkm

0d2c955

brittonsmith approved these changes May 21, 2019

View reviewed changes

brittonsmith reviewed May 21, 2019

View reviewed changes

yt/visualization/fits_image.py Show resolved Hide resolved

munkm merged commit 099ac4a into yt-project:yt-4.0 May 29, 2019

jzuhone deleted the fits branch May 29, 2019 19:36

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Improvements to reading / writing FITS images #2246

Improvements to reading / writing FITS images #2246

jzuhone commented Apr 17, 2019 •

edited

Loading

jzuhone commented Apr 23, 2019

matthewturk left a comment

matthewturk May 14, 2019

matthewturk May 14, 2019

jzuhone May 19, 2019

matthewturk May 14, 2019

matthewturk May 14, 2019

munkm May 15, 2019

matthewturk May 14, 2019

jzuhone May 19, 2019

munkm left a comment

munkm May 15, 2019

jzuhone May 19, 2019

munkm May 15, 2019

jzuhone May 19, 2019

munkm May 15, 2019

jzuhone May 19, 2019

munkm May 19, 2019

jzuhone May 20, 2019

munkm May 15, 2019

jzuhone May 19, 2019

munkm May 19, 2019

munkm May 15, 2019

jzuhone May 19, 2019

brittonsmith left a comment •

edited

Loading

brittonsmith May 21, 2019

munkm May 29, 2019

Improvements to reading / writing FITS images #2246

Improvements to reading / writing FITS images #2246

Conversation

jzuhone commented Apr 17, 2019 • edited Loading

jzuhone commented Apr 23, 2019

matthewturk left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

munkm left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

brittonsmith left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jzuhone commented Apr 17, 2019 •

edited

Loading

brittonsmith left a comment •

edited

Loading