Add RaDec to AltAz transformation tests

[cdeil]

I would like to add tests for conversion to the HorizonSystem.

@dsberry Are there any existing reference results?

If not we could simply use a few observatory locations and a grid of dates?

[dsberry]

On 23 November 2012 09:57, Christoph Deil notifications@github.com wrote:

I would like to add tests for conversion to the HorizonSystem.

@dsberry https://github.com/dsberry Are there any existing reference results?

If not we could simply use a few observatory locations and a grid of dates?

I'm not aware of any reference results myself (that's not to say they do not exist though). Probably best to do as you suggest - use the various packages to convert the table of standard FK5 J2000 positions to (az,el) values assuming some arbitrary geographical position and time, and compare the results. Maybe do this for several different geographical positions and times.

David

[astrofrog]

Since this is not implemented in astropy.coordinates yet, I would say this is lower priority, but it would indeed be good.

[cdeil]

I need this test anyways for another non-public code, so I want to do it now. It won't hurt to have reference results available when someone implements the transformation to HorizontalCoordinates in the future.

I'm thinking of taking the existing input.txt as (lon, lat) positions on Earth (maybe a subset) and for each (lon, lat) checking a small grid of altitudes and observation times, e.g.:

altitudes = [1e3, 6e3, 10e3, 20e3] # unit: km
times = [1900, 1950, 1990, 2000, 2010, 2020]

@dsberry Is there anything I have to be aware of, e.g. that certain coordinates are only defined for certain time ranges, or that certain Earth ellipsoid models are not defined for very low or high altitudes or ... or is this straightforward?

[dsberry]

On 23 November 2012 11:31, Christoph Deil notifications@github.com wrote:

I need this test anyways for another non-public code, so I want to do it now. It won't hurt to have reference results available when someone implements the transformation to HorizontalCoordinates in the future.

I'm thinking of taking the existing input.txt as (lon, lat) positions on Earth (maybe a subset) and for each (lon, lat) checking a small grid of altitudes and observation times, e.g.:

altitudes = [1e3, 6e3, 10e3, 20e3] # unit: km times = [1900, 1950, 1990, 2000, 2010, 2020]

@dsberry https://github.com/dsberry Is there anything I have to be aware of, e.g. that certain coordinates are only defined for certain time ranges, or that certain Earth ellipsoid models are not defined for very low or high altitudes or ... or is this straightforward?

We've only used horizon coords for ground-based work, so I'm not sure about the altitude limits. The notion of an horizon becomes a bit abstract when you are 20,000 km high!

Make sure you specify the time correctly - for instance pyast expects the epoch to be given in the TDB timescale, although using the TT timescale (=TAI + 32.184) instead should introduce errors of no more than 0.03 arc-sec. And pyast also needs the dut1 value at the epoch to be supplied. This is the UT1-UTC correction, in seconds. See http://www.starlink.rl.ac.uk/docs/sun211.htx/node486.html.

For each earth position, how many sky positions are you going to test?

Information about AST's (az,el) coordinate system is available at http://www.starlink.rl.ac.uk/docs/sun211.htx/node612.html#System

(search for "AZEL").

David

[astrofrog]

@cdeil - FYI, you now have push access to this repository, since you've been doing so much work on it!

[cdeil]

In the attached generate_observers.py I take 100 random (lon, lat) and for each one these altitude and time values:

# Sample of altitudes (distances from the Earth's center) in km
altitude = np.array([6e3, 10e3], dtype='int')

# Sample of UTC times
time = ['1900-01-01', '1950-01-01', '1990-01-01', '2000-01-01', '2010-01-01', '2020-01-01']

When running the script I get this warning:

WARNING: iauUtctai: dubious year for UTC (before 1960.0 or 5 years beyond last known leap second) [astropy.time.core]

So there might be problems for pre-1960 or post-2012 times. Should we simply not test those?

[dsberry]

On 23 November 2012 13:00, Christoph Deil notifications@github.com wrote:

In the attached generate_observers.py I take 100 random (lon, lat) and for each one these altitude and time values:

Sample of altitudes (distances from the Earth's center) in km

altitude = np.array([6e3, 10e3], dtype='int')

Sample of UTC times

time = ['1900-01-01', '1950-01-01', '1990-01-01', '2000-01-01', '2010-01-01', '2020-01-01']

When running the script I get this warning:

WARNING: iauUtctai: dubious year for UTC (before 1960.0 or 5 years beyond last known leap second) [astropy.time.core]

So there might be problems for pre-1960 or post-2012 times. Should we simply not test those?

I think so. The current formulation in AST is based on UTC and leap seconds, so it's problematic to use times before UTC began (in 1960) or too far into the future (when other leap seconds may have been introduced).

David

[cdeil]

I've added a few horizontal coordinates computed with pyephem to get started. I'll look tomorrow how to compute them with AST and TPM.

[dsberry]

On 23 November 2012 19:46, Christoph Deil notifications@github.com wrote:

I've added a few horizontal coordinates computed with pyephem to get started. I'll look tomorrow how to compute them with AST and TPM

For AST, it's just a case of setting the ObsLon, ObsLat, ObsAlt and DUT1 attributes in the SkyFrame (see http://starlink.jach.hawaii.edu/docs/sun211.htx/node651.html), and then setting the SkyFrame's System attribute to "AZEL" . The rest is just the same as for the other systems.

David

[cdeil]

@dsberry Is there some python package that let's me compute DUT1 for a given UTC time?

[cdeil]

Note to self: Before doing this full test, first try to get approximately equal results from pyast, pytpm and pyephem for this one case: http://phn.github.com/pytpm/conversions.html?highlight=azimuth#fk5-equinox-and-epoch-j2000-0-to-topocentric-observed

[astrofrog]

@cdeil - about DUT1, if astropy.time (http://docs.astropy.org/en/latest/time/index.html) doesn't do it, should it be implemented?

[dsberry]

On 26 November 2012 13:13, Christoph Deil notifications@github.com wrote:

@dsberry https://github.com/dsberry Is there some python package that let's me compute DUT1 for a given UTC time?

— Reply to this email directly or view it on GitHubhttps://github.com/astropy/coordinates-benchmark/pull/10#issuecomment-10714782.

DUT1 is determined experimentally. Historical values are available at

http://data.iers.org/products/214/14443/orig/eopc04_08_IAU2000.62-now

AST does not incorporate this table since it gets extended (in principle) every day. Instead, AST expects the DUT1 value to be supplied with the data, in the same way that epoch, equinox, etc, is supplied with the data.

David

[cdeil]

If I understand correctly, DUT1 = UT1 - UTC is not implemented in astropy.time at the moment. If it's needed to compute Horizontal coordinates, we will want it sooner or later in astropy, right? @taldcroft and others: should I open an issue?

I see this: https://github.com/astropy/astropy/issues/352

[taldcroft]

@cdeil - yes, you should go ahead and open an issue. IIRC this is nominally assigned to @iguananaut since he had someone working on an implementation that would handle getting updated tables automatically. I'm not sure about the status of that work.

[embray]

I'm not so sure either--it was an intern who I think may be slowly, painfully flaking out. I should have assigned an easier task but it's hard for me to judge what "easy" is--I thought it would be. They did come up with some working regular expressions that I might be able to salvage, but otherwise if you want to move ahead with that we can do it.

[eteq]

I believe astropy/astropy#351 is the issue for this (I re-assigned from astropy/astropy#502, as I believe this is a duplicate).

[cdeil]

I rebased this PR for sky to horizon coordinate checks and plan to work on this as soon as sky to horizon coordinate transformations become available in Astropy.

@eteq is there a branch or issue I can watch for sky to horizon coordinate transformations in Astropy?

[eteq]

@cdeil not yet because it requires the ERFA vectorizing PR to be merged. But I will be sure to @ mention you as soon as there is something to test.

[cdeil]

@eteq I had a quick look at the ERFA vectorizing PR ... it's huge and almost impossible to review on Github. Thanks for all that work!

[eteq]

Yeah, it's a bit annoying, but we had to check some auto-generated code because it wasn't clear how to integrate it into the setup. The hope is to get that in for 1.1 or something so that future reviews were be easier. I can't take all the credit, though - @jwoillez did a lot of the nitty-gritty templating and numpy stuff!

[astrofrog]

@cdeil - it would be great if you could manage to get pyast and pyephem (or any other two codes) to at least roughly agree before we do it in astropy. That is, we want to make sure we have independent reference values that are validated by at least two packages. So if possible, we should try and get this PR to work and merge it before @eteq finishes implementing it in astropy.

[eteq]

@cdeil - one thing that will be in astropy (possibly in v1.0, possibly later, depending on time constraints) will be at least an option for atmospheric refraction corrections. So that requires a few more observer parameters:

1. pressure
2. temperature
3. humidity
4. wavelength of observation

Of course if pyast and pyephem don't support these it's not that critical to get them in here for comparison, but it might be useful to at least have an option for accuracy-checking purposes.

[dsberry]

pyast corrects for diurnal aberattion but not for atmospheric refraction or polar motion.

[cdeil]

@eteq In this PR I only want to do obs.pressure = 0 ... as a reminder I've opened a new issue for the extensions you suggest: #40

[astrofrog]

@cdeil - can you add the changes from summary so we can see what the level of agreement is?

[eteq]

@cdeil - ok, that's fine too. It's possible the default will be to have preesure = 1 atm (or whatever the altitude indicates), but that can be dealt with once it's in the code. So this is fine by me for now.

[astrofrog]

@cdeil - will you have time to get this to work today or tomorrow, so that we can test @eteq's implementation once it's ready?

[cdeil]

Yes, I can work on this tomorrow.

One issue I'm having is that I don't have a verified test case from the literature or another software package and I've been getting conflicting results from different packages. But I don't remember the details ... I'll have another look and give this another try tomorrow.

[astrofrog]

@cdeil - thanks!

[cdeil]

Note to self: Astropy now has AltAz - RaDec transformations ... add precision tests here using https://github.com/astropy/astropy/pull/3217

[astrofrog]

@cdeil - regarding your previous comment, I wouldn't worry too much about having literature 'reference' values for now - if we can get astropy to agree with both codes already here then it's already a good indication we are doing something right. I think for the initial astropy merge 1" accuracy would be plenty.

[astrofrog]

@cdeil - it looks like the values from the two tools here are still very different? Do you have any ideas what might be going on?

[cdeil]

No ... I didn't have time to look at this and figure it out yet. I want to look at composite models now and review that, but I might have time to come back to this in a few hours.

[cdeil]

@astrofrog - Working on this now

[cdeil]

Astropy now gives roughly consistent results with PyEphem.

I did not update the pyast tests, because I'd first like to re-structure coordinates-benchmarks into a Python package, so that we can have a utils.py and import shared functions / constants from there.

@astrofrog - OK to merge this now and continue in new PRs? Do you want to review this?

[astrofrog]

@cdeil - thanks! Would it be possible to commit the changes to the summary files before we merge this PR?

[cdeil]

@astrofrog - The summary and ./run_benchmark.py --tasks horizontal isn't implemented yet. Fo now you have to do this:

$ cd tools/astropy
$ python convert_horizontal.py
$ cd tools/pyephem
$ python convert_horizontal.py

I will continue working on coordinated-benchmarks this week, but I'd like to do it in separate PRs.

[astrofrog]

Ok, then feel free to merge this for now - thanks!