cbyrohl / scida

scida is an out-of-the-box analysis tool for large scientific datasets. It primarily supports the astrophysics community, focusing on cosmological and galaxy formation simulations using particles or unstructured meshes, as well as large observational datasets. This tool uses dask, allowing analysis to scale.
https://scida.io
MIT License
26 stars 4 forks source link

plotting ra,dec from eFEDS FITS file (without units?) #150

Open dnelson86 opened 7 months ago

dnelson86 commented 7 months ago

I copied the example code for the SDSS FITS file, for a different file:

import matplotlib.pyplot as plt
import numpy as np
import scida

path = '/virgotng/mpia/obs/eROSITA/eFEDS/'
filename = 'fm00_merged_020_EventList_c001.fits'

ds = scida.load(path + filename)

fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(111, projection="aitoff")
sc = ax.scatter(ds['RA'], ds['DEC'], s=0.05, c=ds['PI'], rasterized=True)
fig.colorbar(sc, label="Photon Energy [eV]")
ax.set_xticklabels(['14h','16h','18h','20h','22h','0h','2h','4h','6h','8h','10h'])
ax.set_xlabel("RA")
ax.set_ylabel("DEC")
ax.grid(True)
plt.savefig("eFEDS.png", dpi=150)

which produces the following image:

eFEDS

but the distribution of RA,DEC in this file (in theory) are localized to a small sky region like this:

Footprint-of-eFEDS-given-by-the-black-solid-line-Cluster-candidates-present-in-the

Is this possibly because the fields in this FITS file have no units?

WARNING:scida.io._base:FITS file support is work-in-progress.
WARNING:scida.interfaces.mixins.units:Cannot determine units from neither unit file nor metadata for '/DEC'.
WARNING:scida.interfaces.mixins.units:Cannot determine units from neither unit file nor metadata for '/PI'.
WARNING:scida.interfaces.mixins.units:Cannot determine units from neither unit file nor metadata for '/RA'.

or possibly something else?

dnelson86 commented 7 months ago

Documentation for the structure of such files:

https://erosita.mpe.mpg.de/dr1/eSASS4DR1/eSASS4DR1_ProductsDescription/extensions_description_dr1.html

dnelson86 commented 7 months ago

Same result with:

sc = ax.scatter(ds['RA']*ds.ureg['deg'], ds['DEC']*ds.ureg['deg'], s=0.05, c=ds['PI'], rasterized=True)

WARNING:pint.util:Calling the getitem method from a UnitRegistry is deprecated. use `parse_expression` method or use the registry as a callable.
WARNING:pint.util:Calling the getitem method from a UnitRegistry is deprecated. use `parse_expression` method or use the registry as a callable.
/u/dnelson/.local/envs/py3/lib/python3.9/site-packages/numpy/ma/core.py:2820: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  _data = np.array(data, dtype=dtype, copy=copy,
/u/dnelson/.local/envs/py3/lib/python3.9/site-packages/matplotlib/axes/_axes.py:4353: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
  c = np.asanyarray(c, dtype=float)
cbyrohl commented 7 months ago

One can test this without units:

ds_nounits = scida.load(path + filename, units=False)
ds = scida.load(path + filename)

which gives the same result except for carrying "unknown units" as a pint.Quantity for

assert np.all(ds.data["RA"].magnitude.compute() == ds.data["RA"].compute())

I would expect this problem to be unrelated to the unit attachment?

dnelson86 commented 7 months ago
In [5]: RA.min()
Out[5]: 126.20363814596136

In [6]: RA.max()
Out[6]: 145.82010938548257

In [7]: DEC.min()
Out[7]: -3.223283315519501

In [8]: DEC.max()
Out[8]: 6.099576083087285
dnelson86 commented 7 months ago

The result is the same with ds = scida.load(path + filename, units=False).

cbyrohl commented 7 months ago

Isn't the computed data range correct?

dnelson86 commented 7 months ago

Yes, just the visual plot is wrong. Turning off the aitoff projection with ax = fig.add_subplot(111) fixes everything, and produces this image (7.6M photons in eFEDS):

eFEDS4

Is this just some matplotlib bug, or some strange interaction between scida and matplotlib?

cbyrohl commented 7 months ago

It is the same when you pass these arrays cast to numpy, so I do not think this has anything to do with scida.

Regarding units for eFEDS: If there is no metadata, we do not make an assumption for the units here. We can provide a minimal unit file.

dnelson86 commented 7 months ago

With correct x-axis, and sorting by photon energy:

eFEDSd

dnelson86 commented 7 months ago

All 172M photons of the full eRASS1, unfortunately I think (vs Merloni+24 Fig 6) it suffers from a similar bug, i.e. the data should not fill the sky:

eRASS1

cbyrohl commented 7 months ago

Did you check the matplotlib docs? I could imagine it uses radians rather than degrees.