# Run with jystilts (java -jar jystilts.jar)
# Jar file: http://www.starlink.ac.uk/stilts/jystilts.jar

from stilts import *

# Assembles monolithic TGAS FITS file from 16 TGAS chunks at ESAC.
basename = ("http://cdn.gea.esac.esa.int/Gaia/tgas_source/fits/" +
            "TgasSource_000-000-%03d.fits")
chunks = [tread(basename % (i)) for i in range(16)]
tgas = tcat(chunks)

# Plot all points on all-sky projection - looks like a density map
# if you don't zoom in.
plot2sky(projection="aitoff",
         layer1="mark", in1=tgas, lon1="ra", lat1="dec", color1="cyan")

# Plot weighted density map, with quantitative colour axis,
# giving the mean number of photometric observations over the
# whole sky.  Use an ecliptic projection.
plot2sky(projection="aitoff", datasys="equatorial", viewsys="ecliptic",
         layer1="mark", in1=tgas, lon1="ra", lat1="dec",
         shading1="weighted", weight1="phot_g_n_obs", combine1="mean")

# .. think up your own variations.