TypeError due to astropy

[adonini]

Lunching the first script write_dl1.py with the files from the new divergent simulations I get the following error due to astropy when creating reco_result:

Traceback (most recent call last):
  File "i/protopipe/protopipe/scripts/write_dl1.py", line 442, in <module>
    main()
  File "/protopipe/protopipe/scripts/write_dl1.py", line 222, in main
    ) in preper.prepare_event(source, save_images=args.save_images):
  File "/protopipe/protopipe/pipeline/event_preparer.py", line 1163, in prepare_event
    for tel_id in point_altitude_dict.keys()
  File "/protopipe/protopipe/pipeline/event_preparer.py", line 1163, in <dictcomp>
    for tel_id in point_altitude_dict.keys()
  File "/miniconda3/envs/protopipe/lib/python3.7/site-packages/astropy/coordinates/sky_coordinate.py", line 257, in __init__
    frame_cls(**frame_kwargs), args, kwargs)
  File "/miniconda3/envs/protopipe/lib/python3.7/site-packages/astropy/coordinates/sky_coordinate_parsers.py", line 244, in _parse_coordinate_data
    valid_components.update(_get_representation_attrs(frame, units, kwargs))
  File "/miniconda3/envs/protopipe/lib/python3.7/site-packages/astropy/coordinates/sky_coordinate_parsers.py", line 589, in _get_representation_attrs
    valid_kwargs[frame_attr_name] = repr_attr_class(value, unit=unit)
  File "/miniconda3/envs/protopipe/lib/python3.7/site-packages/astropy/coordinates/angles.py", line 618, in __new__
    self.wrap_angle = wrap_angle
  File "/miniconda3/envs/protopipe/lib/python3.7/site-packages/astropy/coordinates/angles.py", line 654, in wrap_angle
    self._wrap_internal()
  File "/miniconda3/envs/protopipe/lib/python3.7/site-packages/astropy/coordinates/angles.py", line 645, in _wrap_internal
    super().__setitem__((), value)
  File "/miniconda3/envs/protopipe/lib/python3.7/site-packages/astropy/units/quantity.py", line 1060, in __setitem__
    self.view(np.ndarray).__setitem__(i, self._to_own_unit(value))
  File "/miniconda3/envs/protopipe/lib/python3.7/site-packages/astropy/units/quantity.py", line 1372, in _to_own_unit
    raise TypeError("cannot convert value type to array type "
TypeError: cannot convert value type to array type without precision loss

If I use a file from the old simulations the script run without error. Both the simulations were done using the same CORSIKA and SIMTEL version, only difference is the add of a TELESCOPE 0 in the simtel conflict file that correspond to the array pointing. I checked and this should not interfere, I have 19 telescope with ID from 1 to 19.
The value of point_azimuth_dict[tel_id] and point_altitude_dict[tel_id] are both float32. One difference with the old simulation is that I have some negative value for the azimuth.

point_altitude_dict[tel_id]:  1.2257447242736816 rad float32
point_azimuth_dict[tel_id]:  -3.0990066528320312 rad float32

Maybe this can be the problem for astropy? @vuillaut @maxnoe
The paths to both new and old simulation on the GRID can be find on redmine here https://forge.in2p3.fr/issues/36268

[adonini]

An update: I have produced new divergent simulation with only positive values for azimuth angles and I didn't receive any error and the dl1 generation was successful. For some reasons astropy is not happy with the negative values.

[maxnoe]

That's strange, as it is normally quite happy with negative azimuth, you can even decide yourself if you want to have it between (-180, 180) or (0, 360):

In [5]: s = SkyCoord(alt=[70, 70, 70], az=[-90, 10, 270], unit=u.deg, frame='altaz')

In [6]: s
Out[6]: 
<SkyCoord (AltAz: obstime=None, location=None, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
    [(270., 70.), ( 10., 70.), (270., 70.)]>

In [7]: s.az.wrap_at(180 * u.deg)
Out[7]: <Angle [-90.,  10., -90.] deg>

In [8]: s.az.wrap_at(360 * u.deg)
Out[8]: <Angle [270.,  10., 270.] deg>

[maxnoe]

I also cannot make out which part of the code raises the error, you traceback does not match the current master.

Which branch are you using?

[adonini]

The branch is the same as the master, no modification, and for the new simulations I simply added 360deg to the negative values. The problem appears at the line 1154 in event_preparer.py. I added some lines for debugging that's why they don't match in the output.

[maxnoe]

What is:

type(point_altitude_dict[tel_id]) ?

If it is numpy.ndarray or astropy.Quantity, what is

point_altitude_dict[tel_id].dtype?

[adonini]

The values are all float32. And the first is class astropy.units.quantity.Quantity

[maxnoe]

Ok, then I think you should convert these to float64 before trying to convert, best already when reading in.

With .astype(np.float64)

[adonini]

ok that solved the problem I analysed one file successfully. I don't see why only the negative values need so much precision. Thanks anyway.

[HealthyPear]

If you want you can open a PR to fix this, then I can keep it open until I update the master with the API change that is currenty in the draft PRs

[adonini]

If we think about real data, it's possible to have negative values, but maybe it's best to transform them to positive ones if the negative ones need more memory.

[maxnoe]

No, I don't think it has anything to do with the sign. It is just generally true that not all float64 values are exactly representable as float32, so this will be always an issue when reading real world floating point data as float32 and dowcasting with the rule safe:

In [10]: np.array([3.2137982173819273192739127319]).astype('float32', casting='safe')
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-10-f6231f54925a> in <module>
----> 1 np.array([3.2137982173819273192739127319]).astype('float32', casting='safe')

TypeError: Cannot cast array data from dtype('float64') to dtype('float32') according to the rule 'safe'

So the sign is not the problem here.

[maxnoe]

I fail to reproduce this, what is your astropy version?

[adonini]

It's actually 3.2.3. I didn't update the environment installation since the master is still using ctapipe 0.7. So I stayed with the version I got when installing protopipe.

[maxnoe]

seems to be fixed in astropy 4

[HealthyPear]

I'll keep this issue open until the next minor release.

There we should use astropy 4 by default because, if I am not wrong, ctapipe right now requires the latest astropy up to 5 excluded.

[maxnoe]

I could not even reproduce this error with simple examples using astropy 3.

[HealthyPear]

ok that solved the problem I analysed one file successfully. I don't see why only the negative values need so much precision. Thanks anyway.

@adonini could you open a PR with what you did from Max's advice if it solved this problem?

Also, I did not understand if you did it only for the negative values or to all coordinates.
The second way seems safer I think

[adonini]

At the end for my results I used the "positive" simulations, that didn't need any modification of the code. Anyway I was able to produce DL1 files converting only the negative values for the older simulations. I can do some tests and see if also DL2 files are produced. But I would wait for the minor release and see how it goes with the newer version of astropy for a PR. No need to take up more memory than necessary.

[adonini]

I could not even reproduce this error with simple examples using astropy 3.

I can give you one of my files if you want.

[maxnoe]

Don't worry please about the memory here. If the simplest solution is converting to float64, do that.

[HealthyPear]

@adonini Could you post here the LFN of a file that causes this crash?

The test file that I added for ctapipe was one of these 2 (sorry I do not remember exactly the run number)
and it came from the page hat you linked

gamma_20deg_180deg_run1___cta-prod3-demo-2147m-LaPalma-baseline.simtel.gz
gamma_20deg_180deg_run9___cta-prod3-demo-2147m-LaPalma-baseline.simtel.gz

[adonini]

If you took them before August, then is the very first sample and it was ok. Every pointing had positive angles. You can retrieve the files you want from here /vo.cta.in2p3.fr/MC/PROD3/LaPalma/gamma/simtel/2233/Data/000xxx . It's the simulation for 1.5x the nominal FoV.