Message Boards Message Boards

Issues with AstronomicalData and SunPosition

Posted 10 years ago

TL;DR

I have three issues with getting sun positions in Mathematica V10:

  • It is extremely slow compared to V9
  • It doesn't seem to yield correct results
  • a) It needs an Internet connection, which b) is not assured to always return results and c) if not used optimally can easily consume your monthly allowance of W|A calls

Long story

With the advent of V10 AstronomicalData has been deprecated, as shown on its documentation page: enter image description here.

I'm not an astronomer, so my main use of this function has been restricted to its capability to get the sun position using functions calls like this:

{
  AstronomicalData["Sun", {"Azimuth", {2013, 3, 1, #, 0, 0}, {52.37`, 4.89`}}, TimeZone -> 1], 
  AstronomicalData["Sun", {"Altitude", {2013, 3, 1, #, 0, 0}, {52.37`, 4.89`}},  TimeZone -> 1]
} & /@ Range[0, 23]

{{341.47732, -43.93930}, {2.33417, -45.21747}, {22.94232, -43.19133}, {41.52167, -38.28400}, {57.55208, -31.30276}, {71.47253, -23.03159}, {84.02194, -14.08294}, {95.91940, -4.92723}, {107.80166, 4.03418}, {120.23770, 12.40167}, {133.72303, 19.72563}, {148.59433, 25.48377}, {164.84179, 29.12209}, {181.93103, 30.19428}, {198.91284, 28.55117}, {214.89602, 24.41962}, {229.46400, 18.28613}, {242.70064, 10.70703}, {254.98998, 2.19073}, {266.84760, -6.82566}, {278.85594, -15.94505}, {291.66723, -24.75464}, {306.00911, -32.75641}, {322.57956, -39.29877}}

So, this gets me the sun positions in steps of an hour during a particular day in Amsterdam (TZ 1).

The same call still works in V10, though it now returns numbers with units; degrees in this case. On its first call, it reads some paclet information from a Wolfram server, but on any successive call no Internet connection is needed. I will be going into detail about timing further on, but I'll say here that the V10 function takes about three times longer than its V9 namesake. I blame the addition of units for that.

With AstronomicalData apparently deprecated we are supposed to use its successors. In this case I need SunPosition. A direct translation of the above would be:

SunPosition[GeoPosition[{52.37`, 4.89`}], DateObject[{2013, 3, 1, #, 0, 0}, TimeZone -> 1]] & /@ Range[0, 23]

{{95.7, -5.1}, {107.6, 3.9}, {120.0, 12.3}, {133.5, 19.6}, {148.3, 25.4}, {164.5, 29.1}, {181.6, 30.2}, {198.6, 28.6}, {214.6, 24.5}, {229.2, 18.4}, {242.5, 10.9}, {254.8, 2.4}, {266.6, -6.7}, {278.6, -15.8}, {291.4, -24.6}, {305.7, -32.6}, {322.2, -39.2}, {341.3, -43.5}, {2.0, -44.8}, {22.5, -42.9}, {41.1, -38.0}, {57.1, -31.1}, {71.0, -22.9}, {83.6, -13.9}}

As with the new AstronomicalData the output is actually in degrees which I have removed in the above output for the sake of clarity. There are a few things to note:

  • SunPosition uses position and date objects, the latter being new in V10
  • SunPosition does not have a TimeZone option, but you can set it in DateObject
  • SunPosition can use the old lat/long list position indication. It also can use a date list to enter the date instead of a DateObject. In the latter case you are out of options with respect to time zones and you have to add the appropriate amount of time offset
  • It is extremely slow, and it may even time-out:

enter image description here

  • Last but not least: the results seem to be plain wrong. It suggests that sunrise is somewhat before 1 am, which is -of course- incorrect. I assume that this has something to do with a $GeoLocation setting for the observer of the sun positions, but I haven't managed to sort out what I am supposed to enter to get the correct sun positions for the location provided in the same call.

As to timing: I noticed very inconsistent timings for SunPosition compared to AstronomicalData, so I used the following code to collect a somewhat more statistical sound sample:

SetAttributes[timingTest, HoldFirst];
timingTest[code_, repeats_Integer] :=
   Table[
      ClearSystemCache[];
      code // AbsoluteTiming // First,
      {repeats}
    ]

Using this, I collected timing of 20 calls to the following code snippets:

  • AstronomicalData V9 and V10: As above
  • SunPosition: As above
  • SunPosition without GeoPosition, just the lat/long list.
  • SunPosition without GeoPosition, and also without the DateObject date, just a classical date list (with the hour set to +1 to accommodate TZ 1)
  • SunPosition V10 without GeoPosition and with the Map (/@) gone and replaced by a DateRange inside the call.

In the last case, the returned value is a TimeSeries object from which I extract the positions using the "Paths" method:

 SunPosition[{52.37`, 4.89`}, DateRange[{2013, 3, 1, 1, 0, 0}, {2013, 3, 1, 24, 0, 0}, "Hour"]]["Paths"][[1, All, 2]]

The results were as follows:

enter image description here

Clearly, the SunPosition results are very disappointing. Getting the sun positions with SunPosition is almost 40 times slower than using the old V9 method (which, I should add, wasn't particularly quick either. I have an implementation in Mathematica code which is faster). The V10 implementation of AstronomicalData is also more than three times slower than the V9 version. The DateRange version of the call saves a lot of communication overhead. Still, it is almost five times slower than in V9.

The cause of all this slowness is that SunPosition simply does a call to Wolfram|Alpha. Sniffing the communication one sees the following string passed to the server:

"1:eJxTTMoPSuNgYGAoZgESPpnFJcHcQEZwaV5AfnFmSWZ+XhoTsmxR/6GvGjH9wg4Qhr6XQxobsnzmXXYGhkxmIC+TEUSIgwggZihigIJgoAIGj/yizKr8PJigA5yBZtubwB1yrdxeDkXVIuvcH1aJOBRzAqUcS0vycxNLMpMBSAArww=="

which can be turned into readable form using Uncompress:

{"SunPosition", {4.89, 52.37}, {2013, 3, 1, 23, 0, 0.}, "Horizon", 2., 2., {52.09, 5.12}, Automatic}

Here, we can recognize the lat/long of the position I used (but with lat/long reversed - Is this somehow significant?). At the end is my own $GeoLocation, but I don't believe it is used at all (and it shouldn't: I'm asking for the sun position over Amsterdam, not where I live). Changing it with Block I get the same results:

Block[{$GeoLocation = GeoPosition[{52.37`, 40.89`}]}, SunPosition[{52.37`, 4.89`}, {2013, 3, 1, 1, 0, 0}]]

Apart from the slowness, there's the issue of the necessary Internet connectivity (Want to give a demonstration and you don't have Internet? Sorry, you're out of luck).

And what of the use of W|A calls? Each of the SunPosition tests (except the last one) took me 20 * 24 = 480 calls. So this part of my testing only already took 1440 calls, and one should be reminded that a typical Home Use license allows for only 3,000 calls per month. Things like this can go pretty fast. In fact, I once wrote an application that calculates the impact of building changes on shadows around your house throughout the year. It does in the order of 17,000 AstronomicalData calls. I couldn't implement that naively using SunPosition and have it actually work. Clearly, one should now use the DateRange version of the call as much as possible.


To wrap up: I have one real question, i.e., how to get SunPosition to return the same values as AstronomicalData, and a request to the WRI team: please put SunPosition in the kernel and don't use W|A calls, because the situation as it is now is rather annoying and IMHO a real step backwards.

POSTED BY: Sjoerd de Vries
27 Replies

I cannot comment on the accuracy of the model used as I was not involved in that implementation. I can only confirm that NASA and the US Naval observatory were not the sources used. This does not mean the results are wrong, just that it uses different algorithms. There are many different algorithms for such computations and some are more accurate than others, may be slower than others, and are more/less accurate over longer date ranges. I believe most of the code used in AstronomicalData/SunPosition came from Meeus and/or VSOP87 and are listed in the sources used by AstronomicalData:

http://reference.wolfram.com/language/note/AstronomicalDataSourceInformation.html

As for atmospheric refraction, as documented, it is not taken into account. Its possible this could be made available in a future release, but its not currently available. Wikipedia and other sources mention several correction factors that can be used. You should be able to apply one of those of your choosing to modify the results of SunPosition in the mean time. As with all such algorithms, they also vary in accuracy as modeling the dynamic effects of the atmosphere is highly non-trivial and location/temperature/pressure dependent. http://en.wikipedia.org/wiki/Atmospheric_refraction

If you are looking for a way to make AstronomicalData/.SunPosition match NASA/Naval Observatory, I doubt that is possible without a total reimplementation to use whatever algorithm they use and there are no current plans to do that.

POSTED BY: Jeffrey Bryant
Posted 10 years ago

Thank you. And thanks for the link. Interestingly one of the datasources was the NASA one.

POSTED BY: John Drummond
Posted 10 years ago

Jeff, thanks for that. => Would you know what the accuracy of the SunPosition Calculations would or should be?

Here's the US Naval Laboratory for another comparison http://aa.usno.navy.mil/ they only have per minute on the site, but the numbers agree with the NASA 17:13 20.2 279.0 17:14 20.0 279.2 17:15 19.8 279.4 17:16 19.6 279.5

As Sjoerd kindly pointed out there's Mathematica did not take Atmospheric Refraction into consideration which will be something like 10 arc minutes at 5 degrees of altitude 1 at 45 and 35 at the horizon, or at 20 degrees of altitude something like 3 arc minutes (using Bennett's 1982 calculation) which is about 0.05 degrees. However this should increase the perceived altitude. So for the SunPosition Calculation this correction is big enough in magnitude, nor is it in the right direction to help with the reconciliation, as the SunPosition is supposed to exclude it.

=> Any thoughts how I reconcile these differences or am I missing something obvious?

=> Would you know what is supposed to be the most accurate estimate (I understand some approximation would be made for atmospheric refraction)?

POSTED BY: John Drummond

One reason for the minor difference between AstronomicalData and SunPosition is that the older AstronomicalData code base uses a simplistic conversion to Julian date for its internal computations. SunPosition/Alpha use a more modern/less naive conversion. The result being that you get a slightly different answer with the newer result less affected by errors in this particular component.

POSTED BY: Jeffrey Bryant
Posted 10 years ago

Are you sure the output is ok? I've refreshed with RebuildPacletData[]

If http : // ssd.jpl.nasa.gov/horizons.cgi?s_loc = 1 #top is the gold standard Then for the sample point I looked at {43.953074,4.472485} and at 2014-07-22 17:15:00 UTC AstronomicalData["Sun",.... does not agree with SunPosition[GeoPosition and both are different from the nasa calculation e.g. for the altitude for this point respectively I get 19.85892, 20.06311334010335, and 19.7609

Unless I've done some thing stupid I can't see that 0.3 degrees out is much of a success


copied from mathematica 10 notebook

In[26]:= {AstronomicalData["Sun",{"Azimuth",{2014,7,22,18,15,0},{43.953074,4.472485}},TimeZone->1], AstronomicalData["Sun",{"Altitude",{2014,7,22,18,15,0},{43.953074,4.472485}},TimeZone->1]} Out[26]= {279.28408[Degree],19.85892[Degree]} In[27]:= SunPosition[GeoPosition[{43.953074,4.472485}], DateObject[{2014,07,22,19,15, 0}, TimeZone -> IntegerPart[LocalTimeZone[Entity["City",{"Paris","IleDeFrance","France"}], DateObject[{2014,07,22,19,15, 0},TimeZone-> 0]]]] ]//FullForm Out[27]//FullForm= List[Quantity[279.096879879350734008354638938468698465083.6876323489782243,"AngularDegrees"],Quantity[20.063113340103351466195182502734913517392.5442756941572635,"AngularDegrees"]]

Result of values in the jpl.nasa site 279.3737 19.7609

POSTED BY: John Drummond

Two things: Firstly, I do not really consider NASA to be the gold standard; I used these words only in my reply to Nasser as he seemed to regard them as such.

Secondly, Mathematica delivers the actual sun position, not its apparent position. The latter would imply taking atmospheric refraction into account. At sunrise that amounts to about 0.5 degree, and is at a minimum when the sun is at zenith position. So, 0.3 deg off isn't really much of a problem, or so it seems.

POSTED BY: Sjoerd de Vries
Posted 10 years ago

I suppose I'm still not clear why the two different methods in Mathematica (AstronomicalData and SunPostion) don't agree. If I look at the Azimuth for which refraction presumably doesn't make much difference then I get

AstronomicalData 279.28408
SunPosition 279.09687
and NASA 279.3737

which there's roughly 0.1 degrees between AstronomicalData and NASA and 0.2 degrees between AstronomicalData and SunPosition

=> What would you consider to be the gold standard?

I'm interested as I was wondering how accurate one could plausibly make a sundial for for example the next 100 years and thought it would be fairly easy to set up on Mathematica, but am now not sure. the timezone correction is 4 minutes a degree, i.e. 0.1 degrees is 24 seconds 0.1 degrees is something like 3.5 mm at 2 metres

POSTED BY: John Drummond

Output seems to be OK now, great! Now, if the speed could be improved somewhat...

POSTED BY: Sjoerd de Vries

I spoke too soon. It hasn't pushed to the production server just yet. Should hopefully be on production tomorrow. Sorry if I caused additional confusion.

POSTED BY: Jeffrey Bryant

This issue should be resolved now that the most recent production push has happened. This was a time zone conversion problem. Results should now more closely match the results from AstronomicalData. This issue affected several functions mentioned in this thread and should all be addressed now. Thank you for letting us know!

POSTED BY: Jeffrey Bryant

I have found that, in addition to SunPosition returning the wrong result, DaylightQ and SiderealTime also produce the wrong result. SiderealTime code is almost certainly called by the SunPosition code, so this might be a good place to start looking.

POSTED BY: Andrew Mattingly

I'm investigating. I will let you know if I have found the problem.

POSTED BY: Jeffrey Bryant

Users here get the same results I do. I'm investigating for any issues I can think of that might cause such behavior.

POSTED BY: Jeffrey Bryant

I have now asked two users that are in my time zone, they get the same result as me. Another user is on GMT and his results differ from both yours and mine.

It really seems the calculation take your own locale into account, which it shouldn't.

POSTED BY: Sjoerd de Vries

I checked by setting my time zone to Champaign's

Block[{$TimeZone = -5}, 
  SunPosition[GeoPosition[{52.37`, 4.89`}], DateObject[{2013, 3, 1, #, 0, 0}, TimeZone -> 1]] & /@   Range[0, 23]
]

I now get your results.

So, you get the correct results because part of the calculation somehow assumes your time zone, which yields incorrect results for anyone not in your time zone. You are among the few who can't see the incorrect results ;-)

POSTED BY: Sjoerd de Vries

BTW I get the impression that daylight savings also causes problems. When I ask:

In[22]:= Sunrise[GeoPosition[{52.37`, 4.89`}],  DateObject[{2013, 3, 1}, TimeZone -> 1]]

I get

Out[22]= DateObject[{2013, 3, 1}, TimeObject[{8, 28}], TimeZone -> 2.]

whereas I feel it should be

DateObject[{2013, 3, 1}, TimeObject[{7, 28}], TimeZone -> 1.]

as the location is GMT + 1 and DST was not in effect on March,1 (it started March 31).

The documentation says: "Daylight Saving Time corrections must be included in the time zone", which means that fixed locations change from one time zone to another depending on DST.

POSTED BY: Sjoerd de Vries

I saw your analysis. The API call includes a fair amount of information that may or may not be needed for a specific input. I see nothing obviously wrong with the API arguments being sent assuming for you:

$TimeZone -> 2
$GeoLocation -> GeoPosition[{52.09, 5.12}]

You only gave a single API call, but did not specify which value of date was used as input for that call.

Your input date specifies that you are trying to feed in a different timezone (1) from where you are (2). TimeZone conversions need this information to be correct or you could get confusing results.

POSTED BY: Jeffrey Bryant

Amsterdam is in TimeZone 1. The way Mathematica is set-up for daylight savings time entails it adds an hour if DST is in effect. On march 1, however, DST is not in effect and TZ 1 is correct for both me and Amsterdam. Anyway, setting TimeZone to two still has the sun rise in the middle of the night.

But I think that's beside the matter. With SunPosition you're asking for the position of the sun at a given day and time above a specific location. All these are the same in my call and in yours and should therefore return the same result. They don't.

In the meantime I've asked another user for his results and he comes up with exactly the same results I got.

POSTED BY: Sjoerd de Vries
POSTED BY: Jeffrey Bryant

No, I'm on a standard V10.0.0

I just executed the code again and still get the same result:

In[6]:= SunPosition[GeoPosition[{52.37`, 4.89`}], 
   DateObject[{2013, 3, 1, #, 0, 0}, TimeZone -> 1]] & /@ Range[0, 23]

Out[6]= {{Quantity[95.7, "AngularDegrees"], 
  Quantity[-5.1, "AngularDegrees"]}, {Quantity[107.6, 
   "AngularDegrees"], 
  Quantity[3.9, "AngularDegrees"]}, {Quantity[120.0, 
   "AngularDegrees"], 
  Quantity[12.3, "AngularDegrees"]}, {Quantity[133.5, 
   "AngularDegrees"], 
  Quantity[19.6, "AngularDegrees"]}, {Quantity[148.3, 
   "AngularDegrees"], 
  Quantity[25.4, "AngularDegrees"]}, {Quantity[164.5, 
   "AngularDegrees"], 
  Quantity[29.1, "AngularDegrees"]}, {Quantity[181.6, 
   "AngularDegrees"], 
  Quantity[30.2, "AngularDegrees"]}, {Quantity[198.6, 
   "AngularDegrees"], 
  Quantity[28.6, "AngularDegrees"]}, {Quantity[214.6, 
   "AngularDegrees"], 
  Quantity[24.5, "AngularDegrees"]}, {Quantity[229.2, 
   "AngularDegrees"], 
  Quantity[18.4, "AngularDegrees"]}, {Quantity[242.5, 
   "AngularDegrees"], 
  Quantity[10.9, "AngularDegrees"]}, {Quantity[254.8, 
   "AngularDegrees"], 
  Quantity[2.4, "AngularDegrees"]}, {Quantity[266.6, 
   "AngularDegrees"], 
  Quantity[-6.7, "AngularDegrees"]}, {Quantity[278.6, 
   "AngularDegrees"], 
  Quantity[-15.8, "AngularDegrees"]}, {Quantity[291.4, 
   "AngularDegrees"], 
  Quantity[-24.6, "AngularDegrees"]}, {Quantity[305.7, 
   "AngularDegrees"], 
  Quantity[-32.6, "AngularDegrees"]}, {Quantity[322.2, 
   "AngularDegrees"], 
  Quantity[-39.2, "AngularDegrees"]}, {Quantity[341.3, 
   "AngularDegrees"], 
  Quantity[-43.5, "AngularDegrees"]}, {Quantity[2.0, 
   "AngularDegrees"], 
  Quantity[-44.8, "AngularDegrees"]}, {Quantity[22.5, 
   "AngularDegrees"], 
  Quantity[-42.9, "AngularDegrees"]}, {Quantity[41.1, 
   "AngularDegrees"], 
  Quantity[-38.0, "AngularDegrees"]}, {Quantity[57.1, 
   "AngularDegrees"], 
  Quantity[-31.1, "AngularDegrees"]}, {Quantity[71.0, 
   "AngularDegrees"], 
  Quantity[-22.9, "AngularDegrees"]}, {Quantity[83.6, 
   "AngularDegrees"], Quantity[-13.9, "AngularDegrees"]}}

Did you see my analysis above? It showed that the actual W|A call also includes the caller's location. If that is somehow used in the calculation that would explain why we see different results. However, the results should not depend on that and it seems likely this is a bug.

POSTED BY: Sjoerd de Vries

The speed overhead is a known issue that we are actively planning to address in a future release.

POSTED BY: Jeffrey Bryant

Good to hear Jeff! Any suggestions as to how to get SunPosition return the same results as AstronomicalData?

POSTED BY: Sjoerd de Vries

Very good posts, Sjoerd. For my projects, I've started doing the calculations in other languages, using well-known datasets and algorithms, and then finally feeding the various chunks of calculated data to Mathematica for higher-level assembly and presentation. It's not the most desirable solution, but as Mathematica becomes a higher and higher level language, it's a pragmatic approach, at the moment. Of course, the temptation is to try and do it in just one environment - but then, which one?! :)

For one-off questions or exploring a new problem space, firing off a few dozen WolframAlpha queries is fine, but, like elsewhere in computing, relying on the cloud isn't always going to be the ideal solution. There are some obvious disadvantages, too: you don't know where it comes from or how it's been calculated, how reliable it is, or when or whether it's been updated or modified, missing and unexplained data points can't be tracked down and repaired, and so on. Unless the cloud service is quick and very reliable, and well documented, it's easy to end up with problems that aren't easily fixed.

POSTED BY: C ormullion

Thanks, Cormullion. If you have a use for my converted sun position code just let me know.

POSTED BY: Sjoerd de Vries

I can't answer the other issues since I do not use these functions. But on this point you raise

It needs an Internet connection

How would propose that Mathematica get the current ephemeris data for solar system planets and the sun without using an updated database somewhere on the internet?

I use JPL Horizon website myself for this information. You might want to see if there is a way to use this site via an API call to obtain this information. It is free.

POSTED BY: Nasser M. Abbasi
POSTED BY: Sjoerd de Vries

I converted to Mathematica a VBS macro that you can find in this Excel sheet. It is written by Greg Pelletier and based on a NOAA Javascript program which was used to drive their sun calculation page. I performed a direct and ugly conversion without any optimization whatsoever. Testing it with the method above I found that it is about 30 times faster than V9 AstronomicalData and about 1100 times faster than V10's SunPosition.

Testing the output and taking the NASA website you refer to as the golden standard I get for the sun positions starting at 1 am (using AstronomicalData for V10):

Azimuth:

enter image description here

Elevation:

enter image description here

I declare winner the one with the smallest absolute deviation from the NASA data.

As to the need for ephemeris: The following plot compares the results for Pelletier's and V10's AstronomicalData in the year 3013, thousand years in the future (CircleDot is Mathematica):

enter image description here

POSTED BY: Sjoerd de Vries
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract