Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix method camera_to_sky for real events #421

Merged
merged 7 commits into from
Jun 4, 2020

Conversation

JouvinLea
Copy link
Collaborator

@JouvinLea JouvinLea commented May 26, 2020

This PR fix the problem in the method computing camera_to_sky for event to take properly into account the time of the observation in the computation.

For the moment, just use time if provided and the default obs time in utils if not printing a warning.

I think the solution to use HA/DEC for MC event instead of RA/DEC should be done in another method when you compute ra/dec for those SkyCoords. For example, if I read correctly to compute the HA from RA, we need to extract the sideral time from the obs time and the location.

scale_time = sky_coords.obstime.scale
longitude = sky_coords.location.longitude
latitude = sky_coords.location.latitude
time=Time(obstime, scale, location=(longitude, latitude)
HA = sky_coords.icrs.ra.to("hourangle") -time.sidereal_time("apparent")
dec = sky_coords.icrs.dec

Also I need to improve the test to add a check on a ra/dec value when we add the time.

@codecov
Copy link

codecov bot commented May 26, 2020

Codecov Report

Merging #421 into master will increase coverage by 0.07%.
The diff coverage is 100.00%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #421      +/-   ##
==========================================
+ Coverage   41.46%   41.53%   +0.07%     
==========================================
  Files          76       76              
  Lines        6179     6187       +8     
==========================================
+ Hits         2562     2570       +8     
  Misses       3617     3617              
Impacted Files Coverage Δ
lstchain/reco/tests/test_utils.py 100.00% <100.00%> (ø)
lstchain/reco/utils.py 68.93% <100.00%> (+0.97%) ⬆️
lstchain/tests/test_lstchain.py 100.00% <100.00%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 3995ebb...0668c90. Read the comment docs.

@@ -203,16 +203,24 @@ def reco_source_position_sky(cog_x, cog_y, disp_dx, disp_dy, focal_length, point
return camera_to_sky(src_x, src_y, focal_length, pointing_alt, pointing_az)


def camera_to_sky(pos_x, pos_y, focal, pointing_alt, pointing_az):
def camera_to_sky(pos_x, pos_y, focal, pointing_alt, pointing_az, time = None):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe we should change the name of the function. As @maxnoe pointed out in #420 , this method is converting to a HorizonFrame

maxnoe
maxnoe previously requested changes May 27, 2020
horizon_frame = AltAz(location=location, obstime=time)
else:
horizon_frame = AltAz(location=location, obstime=obstime)
print("No time given, use default observation time. To be use only for MC data.")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do not use print here, if at all, use logging. But It is unnecessary here.



def camera_to_sky(pos_x, pos_y, focal, pointing_alt, pointing_az, time = None):
def camera_to_altaz_frame(pos_x, pos_y, focal, pointing_alt, pointing_az, time = None):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The frame in the end is a little strange (Either its camera_frame_to_altaz_frame or camera_to_altaz).
I would just leave out the frame.

time should be obstime to be consistent with the astropy API, and then you also don't need the if statement below.

Copy link
Collaborator Author

@JouvinLea JouvinLea May 27, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @maxnoe for your usefull input.
Here I disagree with you but may I misunderstood. I think for MC events we need to use the obstime_mc in order to be able to create a posteriori a ra/dec for MC event in HA/dec as proposed @moralejo . To perform that, I understood we will need:

longitude = sky_coords.location.longitude
latitude = sky_coords.location.latitude
time=Time(obstime, scale, location=(longitude, latitude)
HA = sky_coords.icrs.ra.to("hourangle") -time.sidereal_time("apparent")
dec = sky_coords.icrs.dec

So we need to have a false obstime to be able to then transform to HA.

Copy link
Member

@maxnoe maxnoe May 27, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think @moralejo was proposing to actually calculate HA, just mentioning the possibility.

At the moment, there is no usecase I know of for having simulations in another frame than horizon.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that the monte carlo processing should be fully agnostic of any obstime. Introducing a fake obstime is a hack (that should be unnecessary).

If it at some position is necessary to introduce a fake time, I would be glad to know, because that means we nead it fixed.


"""

if time:
Copy link
Member

@maxnoe maxnoe May 27, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This if is not needed. Just use obstime=obstime and it works fine also in the case where obstime is None


Returns
-------
sky frame: `astropy.coordinates.sky_coordinate.SkyCoord`

sky frame: `astropy.coordinates.AltAz`
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually the return value is a astropy.coordinates.SkyCoordinate that has a frame of AltAz.

So maybe write `astropy.coordinates.SkyCoord` in AltAz frame`

pointing altitude in angle unit
pointing_az: `~astropy.units.Quantity`
pointing altitude in angle unit
time: `~astropy.time.Time`
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Call it obstime

@@ -43,8 +44,8 @@

# position of the LST1
location = EarthLocation.from_geodetic(-17.89139 * u.deg, 28.76139 * u.deg, 2184 * u.m)
obstime = Time('2018-11-01T02:00')
horizon_frame = AltAz(location=location, obstime=obstime)
obstime_mc = Time('2018-11-01T02:00')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just remove this.

Why was this added anyway? What is so special about 2018-11-01T02:00 that it shoud be used as time fo mc?

if obstime:
horizon_frame = AltAz(location=location, obstime=obstime)
else:
horizon_frame = AltAz(location=location, obstime=obstime_mc)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again: just don't use an obstime for mc please. It has no meaning.

horizon_frame = AltAz(location=location, obstime=obstime)
else:
horizon_frame = AltAz(location=location, obstime=obstime_mc)
if not obstime:
logging.warning("No time given, use default observation time. To be use only for MC data.")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Now, you don't use any obstime (obstime is None, you are not using the global)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes this is what you ask no?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But the logging message says something about a default obstime, you not using any obstime

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And I would say, as this is the expected situation for every simulated event, this should not be a "warning".

Warnings should be for unexpected conditions that require the user to do something. Meaningless warnings are dangerous because then users start ignoring them and potentially start ignoring also the important warnings.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is because you are too fast, I already modify the print in the time you did your comment

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree for the warning

@maxnoe
Copy link
Member

maxnoe commented May 27, 2020

Btw. the frame to get local hourangle and declination would be the Celestial Intermediate Reference System:
https://docs.astropy.org/en/stable/api/astropy.coordinates.CIRS.html

EDIT: now, the origin of the CIRS is the CIP, so not yet a local coordinate frame.

@maxnoe
Copy link
Member

maxnoe commented May 27, 2020

There is an open issue for requesting the local equatorial frame in astropy here:
astropy/astropy#10239

Copy link
Contributor

@rlopezcoto rlopezcoto left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JouvinLea is this PR ready (all comments taken into account, except this small one I just spotted in one comment)?
Are you also planning to implement the writing of RA/Dec coordinates in the DL1/2 dataframe?

@@ -227,9 +236,13 @@ def camera_to_sky(pos_x, pos_y, focal, pointing_alt, pointing_az):
focal = 28*u.m
pointing_alt = np.array([1.0, 1.0]) * u.rad
pointing_az = np.array([0.2, 0.5]) * u.rad
sky_coords = utils.camera_to_sky(pos_x, pos_y, focal, pointing_alt, pointing_az)
sky_coords = utils.camera_to_altaz_frame(pos_x, pos_y, focal, pointing_alt, pointing_az)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

camera_to_altaz_frame -> camera_to_altaz

@maxnoe
Copy link
Member

maxnoe commented Jun 2, 2020

Is the obstime defined globally in the utils.py file used somewhere else? It should be removed.

@JouvinLea
Copy link
Collaborator Author

if you think it can be merge, I am personally done.
otherwise tell me what I should change.

Obstime at the beginning of the files is used to define the horizon_frame used in other methods . And I will not do this clean up in this PR. maybe one PR can be done to clean the utils.py in the future.

Once this is merge, I will do another PR for implementing ra/dec in the DL2 dataframe but other PR is better because this is another subject.

Copy link
Contributor

@rlopezcoto rlopezcoto left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, I agree with @JouvinLea and we should move on with this. @maxnoe could you open an issue on the cleanup of utils.py according to what we discussed here?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants