Small plot straight from GAIA

[enzodbspace]

Hello,

I prepared this quick example as a proof of concept. It can probably be improved. The script downloads a predefined circular area (centered on RA and Dec, and with a defined radius) of basic star data from the GAIA catalog, builds the Starplot catalog, and generates the map of the area. It requires astroquery to download the catalog.

The catalog query are async (launch_job_async) to avoid being limited to a maximum of 2000 records, and is currently responding very slowly. So you will likely have to wait several minutes. If you want to modify the query (more fields, data filters, etc), please see the Astroquery GAIA documentation.

With the below code, I generated this map centered on the open cluster NGC 3293 (Gem Cluster, in Carina), with 2° radius.

Please take into account the directory settings for the files, the coordinate format, and the diameter (I don't know what happens if it's too large, spacetime rupture?). All this in the config section.

Tip: If you have already download and build the catalog, and only want to modify the map, please comment the first two functions at the end query() and build().

from astroquery.gaia import Gaia
from astropy import coordinates
from astropy.coordinates import SkyCoord
import astropy.units as u
from astroquery import gaia

import csv
import time

from shapely import Point

from starplot import Star, MapPlot, Path, Mercator, Field, callables, _
from starplot.data import Catalog
from starplot.styles import PlotStyle, extensions


#
# Config
#

data_directory      = "data"                    # Directory for catalog download and catalog final buid
plot_directory      = "output"                  # Directory for plot output
catalog_output_file = "GAIA_temp_catalog.csv"   # GAIA Stars download filename
catalog_build_file  = "GAIA_build_catalog"      # Build catalog name

coor_ra_center      = "10h35m48.0s"             # RA Center
coor_dec_center     = "-58d13m00s"              # Dec Center
circle_radius       = 2                         # Circle area radius (deg) - ! WARNING: be nice, or you'll open a wormhole
magnitude_min       = -2                        # Magnitude Min
magnitude_max       = 14                        # Magnitude Max

plot_output         = "gaia_map.png"            # Starplot Output


# Coordinates in degrees (for Starplot)
coords = SkyCoord(ra=coor_ra_center, dec=coor_dec_center)
coor_ra_center_deg  = coords.ra.degree
coor_dec_center_deg = coords.dec.degree


#
# GAIA Stars Download
#

def query():

	print("GAIA Stars Download")
	
	gaia.ROW_LIMIT = None
	
	# Center Coordinates RA Dec
	c = coordinates.SkyCoord(coor_ra_center + " " + coor_dec_center, frame='icrs')
	
	
	# GAIA Query by radius
	query = f"""
	SELECT source_id, ra, dec, phot_g_mean_mag, bp_rp 
	FROM gaiaedr3.gaia_source 
	WHERE CONTAINS(POINT('ICRS', gaiaedr3.gaia_source.ra, gaiaedr3.gaia_source.dec), 
	               CIRCLE('ICRS', {c.ra.deg}, {c.dec.deg}, {circle_radius})) = 1
	AND gaiaedr3.gaia_source.phot_g_mean_mag >= {magnitude_min}
	AND gaiaedr3.gaia_source.phot_g_mean_mag <= {magnitude_max}
	"""
	
	print("Launching async job: query database. Please wait. This may take several minutes...")
	
	job = Gaia.launch_job_async(query)
	
	result = job.get_results()
	
	row_count = 0
	
	with open(data_directory + "/" + catalog_output_file, 'w') as catalog:
		# First Line
		catalog_row = "GAIA_ID,RA,DEC,MAG_G,BP_RP"
		catalog.write(catalog_row + '\n')
	
		# Results iteration
		for row in result:
			gaia_id     = row['source_id']
			ra          = row['ra']
			dec         = row['dec']
			mag_g       = row['phot_g_mean_mag']
			color_bp_rp = row['bp_rp']
	
			catalog_row = str(gaia_id) + "," + str(ra) + "," + str(dec) + "," + str(mag_g) + "," + str(color_bp_rp)
	
			row_count += 1
			print("# " + str(row_count) + " " + str(gaia_id))
	
			catalog.write(catalog_row + '\n')
	
	print("Entries: " + str(row_count))


#
# STARPLOT Catalog Build
#

def parse_int(value):
	return int(value) if value else None


def parse_float(value):
	return float(value) if value else None


def stars():

	print(data_directory + "/" + catalog_output_file)

	with open(data_directory + "/" + catalog_output_file, "r") as starscat:
		reader = csv.DictReader(starscat)

		for row in reader:

			ra = parse_float(row.get("RA"))
			dec = parse_float(row.get("DEC"))
			geometry = Point(ra, dec)
			name = row.get("GAIA_ID")

			yield Star(
				pk=parse_int(row.get("GAIA_ID")),
				hip=None,
				ra=ra,
				dec=parse_float(row.get("DEC")),
				constellation_id=None,
				magnitude=parse_float(row.get("MAG_G")),
				ra_mas_per_year=0,
				dec_mas_per_year=0,
				bv=None,
				geometry=Point(ra, dec),
				epoch_year=2000,
				bayer=None,
				flamsteed=None,
				name=name,
			)

def build():

	print("Building Star Catalog...")
	time_start = time.time()

	catalog = Catalog(path=data_directory + "/" + catalog_build_file + ".parquet")
	catalog.build(
		objects=stars(),
		chunk_size=200_000,
		columns=[
			"pk",
			"hip",
			"ra",
			"dec",
			"magnitude",
			"bv",
			"ra_mas_per_year",
			"dec_mas_per_year",
			"constellation_id",
			"geometry",
			"ccdm",
			"epoch_year",
			"name",
			"bayer",
			"flamsteed",
		],
		partition_columns=[],
		sorting_columns=["magnitude"],
		compression="snappy",
		row_group_size=100_000,
	)

	duration = time.time() - time_start
	print(f"Done. {duration:.0f}s")


#
# STARPLOT Map Generation
#

def plot():

	print("Map Plot...")

	style = PlotStyle().extend(
		extensions.BLUE_DARK,
		extensions.MAP,
	)

	# Stars
	style.star.label.font_size = 8
	style.star.label.font_weight = 400

	p = MapPlot(
		projection=Mercator(center_ra=coor_ra_center_deg),
		ra_min=coor_ra_center_deg - (circle_radius * 2),
		ra_max=coor_ra_center_deg + (circle_radius * 2),
		dec_min=coor_dec_center_deg - (circle_radius),
		dec_max=coor_dec_center_deg + (circle_radius),
		style=style,
		resolution=3000,
		scale=1.1,
	)

	gaia = Catalog(
		path=Path(data_directory + "/" + catalog_build_file + ".parquet"),
		hive_partitioning=True,
		spatial_query_method="healpix",
	)

	p.stars(
		catalog=gaia,
		show_names=False,
		bayer_labels=False,
		flamsteed_labels=False,
		sql_labels=False,
		where_labels=[_.magnitude <= 4]
	)

	p.gridlines(
		ra_locations=[d for d in range(0 * 15, 24 * 15, 1)],
		dec_locations=[d for d in range(-90, 90, 1)],
		ra_formatter_fn=lambda ra: f"{round(ra * 15)}\u00B0",
	)

	p.export(plot_directory + "/" + plot_output, padding=0.5)


query() # Query GAIA database and download the stars
build() # Build the catalog
plot()  # Plot the map

Enjoy!

[steveberardi]

coool, thanks for sharing!! I'll have to investigate astroquery some more and maybe add a few examples plotting data from it in the examples section.

btw, I've been working on building a Starplot catalog builder for Gaia DR3 -- the initial version is available at https://github.com/steveberardi/starplot-gaia-dr3/releases/tag/v0.1.0. Most queries with it should take a few seconds or less.

also, quick tip: if you put python after the three backticks of the code block in markdown it'll add nice syntax highlighting.

[enzodbspace]

Nice tip!

Yes! I've seen it. I think it's great! I thought the direct download option was for those who want to do something specific without downloading all of GAIA (maybe there could be a reduced version, like mag <= 15 or something similar?) Last year it was very fast, almost no waiting. I think they're making some big changes for DR4.

One idea, but perhaps for another thread, is to have a catalog-agnostic build function: instead of adapting the script to the catalog, the user adapts the catalog to one more generic script (sometimes it's easier to manipulate the data in a spreadsheet) Especially for DSO.

[steveberardi]

maybe there could be a reduced version, like mag <= 15 or something similar?

Yeah, I'm planning to build multiple versions with different magnitude ranges like that, once I get more confirmation that the magnitude/BV calculations are correct and that there aren't any bugs/etc

One idea, but perhaps for another thread, is to have a catalog-agnostic build function: instead of adapting the script to the catalog, the user adapts the catalog to one more generic script (sometimes it's easier to manipulate the data in a spreadsheet) Especially for DSO.

Not sure I follow this. Do you mean have a function something like build_catalog_from_csv ? I've considered something like that, but handling the geometry field can get tricky. I also thought using an iterable of the sky object models you want in the catalog was a good way to ensure objects get cataloged with the correct schema.

[enzodbspace]

Do you mean have a function something like build_catalog_from_csv?

Yes, something like that. I've done it many times with custom catalogs for various applications, like the good old SkyMap Pro. But I understand that the user should do things correctly, otherwise it won't work properly.

(English is not my first language, sorry for any grammatical errors 🙃)