Plotting example with Matplotlib - efficiency?

[ursomniac]

So, I've been following the awesome examples of how to use Skyfield with Matplotlib and have managed to create some things rather quickly:

Finder chart: (note that Mercury is in front of M75 AND the proximity of Venus --- of course it was cloudy...)

Views of planets with moons:

... and last night I did my version of a star chart using the show-and-tell from @bathoorn in another discussion as a guide.

(Note, the Moon shown is just the Unicode symbol: it's not trying to represent the actual phase here.)

SO: I want to improve things. What I'm finding is that the code I grabbed/borrowed from the "Drawing a finder chart for comet NEOWISE" example seems to grab the entire sky BEFORE filtering them for limiting magnitude (and that's the only filter). As a result, running several finder charts can take quite a long time and in my case, for Uranus and Neptune, I only have a FOV of 8° so plotting ~8.000 stars from Hipparcos (or the BSC) when you'll only see ~100 is a lot of wasted computation.

BUT I can't think of a way to improve "pre-filtering" that isn't just "put another layer of expensive computations" in place of "plot everything and just let the scaling window show what's there"; viz.: "only get the stars that are within a FOV / 2. radius of the center (ra, dec) of the projection... and a simple/naïve filter of "ra ± FOV", "dec ± FOV" will fail when you get closer to the poles (and it doesn't matter if you transform to another system), and any other "calculate the distance to the center of the plot" method means you have to do it for every star on the list. I suppose I could test the relative speed of Matplotlib vs. Python doing spherical geometry (I suppose it's just N * 4 or 5 trig function calls - so it MIGHT not be too bad?).

Anyway I'm curious if anyone has any insights?

P.S. Background/CSB:

I'm doing this as part of a Django site I built that has a target DSO model, and a model for potential observing sites. I just purchased an SE8, and have spent the ensuing 7 weeks (!) of cloudy weather, coding. :-)

[brandon-rhodes]

It would be interesting to see a profile showing how much time is spent in each routine, to measure what fraction of your computation time is spent rotating the 8,000 stars. While rotating the entire Hipparcos catalog would be expensive, if you have already pre-filtered the catalog to only 8,000 stars then I would expect the plot coordinates to be generated very quickly, when compared to more expensive operations like the precession and nutation of the Earth's poles. (Or maybe those aren't involved here, because you're using J2000 coordinates instead of the RA and dec of that particular date?)

I think you are correct that, when dealing with a plot that involves such a huge expanse of sky, there might not be a useful way to put the 8,000 stars into groups that make it easy to eliminate the ones that will be off the chart. If I'm reading the chart correctly, it includes more than half of the stars in the sky?

I'm not experienced with optimizing matplotlib operations, but a profile might suggest whether matplotlib itself is part of the slowdown here. I wonder what the speed difference would be if the Python code simply produced an SVG file, and another program was then called to render that to a raster? Or if a non-raster output format were chosen, like delivering SVG as the output in the browser?

If you had a sample snippet of code that draws a sample plot, then I could try profiling it and seeing where the expense goes.

In any case, I like how your plot looks, it uses color effectively to draw my eye to the planets and to the deep space objects!

[ursomniac]

Well, it's a work in progress, but the code is at https://github.com/ursomniac/skytour

I'll try to think up a way to do the performance testing.

[ursomniac]

Ooooh - it might get a little hairy, but there might be a way to use Django F() or Func() expressions to run a transformation from ra/dec to distance, and then do an .exclude() or .filter on that... Of course that'd require putting the Hipparcos data into the database (which isn't ideal).
BUT there might be a way to do it in numpy as well (I'm really new to numpy).

[ursomniac]

Xo I did find something that does speed things up.

I haven't tried it on the Hipparcos dataset yet, but I've added it to my todos.

Basically, you can use scipy's KTree to do a nearest-neighbor tree on a series of object from a set of coordinates.
This is further accelerated by using the chord length between (ra, dec)1 and (ra,dec)*. (Basically, drill a tunnel through the sphere and calculate the distance = sqrt(∂x^2 + ∂y^2 + ∂z^2) for each object. The distances of course won't be the angular distance but the rankings are the same, and you CAN convert from chord length to arc length easily ex post facto.

So seems to me that if add X, Y, Z, columns (on a unit sphere) to anything with (ra, dec), you can quickly slice through them. For example, for the skychart above, on a unit sphere is distance > sqrt(2), then the object is below the horizon. For a field/finding chart you have to work out the FOV size (and here for my 8° FOV chart above, I had to do some fiddling, but it was doable), and it could go through my 1500+ DSO list (with X,Y,Z, stored as a model property) in a fraction of a second, so that I'm now only plotting say 10 things on the chart instead of 1500 with 1590 "out of bounds".

I'll try to write up something more coherent soon.

[brandon-rhodes]

Neat! Though I can't quite imagine exactly what the code looks like, the approach sounds like a good one. I'll have to read about KTree sometime soon, in case it can help more folks using Skyfield.

[ursomniac]

• Yeah - I had to do some fudging, and I still can't figure out why. The chord lengths (somehow) end up being short by a factor of 90-100 for an 8° FOV; if I multiply them by 100, then they're "right" in terms of "which objects are on the plot". So my code is a little kludgy; if I can sort that out I can refactor the code to something more coherent. If I had to guess, the spheres (map vs. chords) have different radii, so I need to confirm that somehow.