diff --git a/README.md b/README.md index 8ddb9bb..6f005f9 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,67 @@ # gaia-astrometry-net -Code to generate astrometry.net index files from the GAIA catalog +*Code to generate astrometry.net index files from the GAIA DR2 catalog* + +## Prerequisites +- numpy +- sqlalchemy +- astrometry +- astropy +- lcogt_logging +- sep +- pyyaml + + +## Installation +All of the dependencies below install automatically when you run python setup.py +install (or via pip) except *astrometry*. The *astrometry* package is not +currently on PyPi and is provided directly by astrometry.net. This can be installed +by cloning the astrometry.net repo at https://github.com/dstndstn/astrometry.net +and running python setup.py install in the top directory. + +Once the *astrometry* package is installed, the code can be installed in the +usual way by running python setup.py install. + +## Usage +There are two main purposes for the code in this repo. The first is to generate +astrometry.net index files from the GAIA DR2 catalog. The second is wrap +astrometry.net to use the new index files to solve for the astrometry of an +image. + +The index file creation logic is mostly in index_file_creator.py. To create the index files +go into an empty directory and run +``` +create_gaia_index_files +``` +from the command line. +This will then use the network mounted of the GAIA DR2 csv files to create astrometry.net +index files. The first issue is that csv files do not cover consistent patches of the sky. +As such, the first thing that we do is build an SQLite database with all of the positions +and fluxes of the GAIA sources. This DB includes two tables: sources and positions. +The positions table is indexed using an r-tree index so that cone searches are +fast (even for the 1.5 billion GAIA sources). Generating the DB can take a few days +of cpu time to complete. The final DB will be ~200 GB. After we have the queryable +catalog (the DB), we split the catalog into healpixels (see https://healpix.jpl.nasa.gov/). +The sources in each healpixel are stored in a binary table extension of a fits file +that astrometry.net can read. Finally, we use astrometry.net's internal tools +to create index files that astrometry.net can use for solves. We create a range +of index files that have quads for many different image scales +(large field of view to small field of view images). The meta data for the index +files is stored in an ascii file called `gaia-dr2-index-files.dat`. + +Once the index files are created, we need to wrap astrometry.net to use them. +The code to do this is mostly in `solve.py`. + +## Examples +To solve an image using the GAIA astrometry.net index files, you can run +``` +solve-astrometry-gaia --filename {filename} --index-file-path {path} --ra {ra} --dec {dec} --radius {radius} +``` +`filename` is the name of the image to solve. +The `index-file-path` argument should point to the directory with the GAIA astrometry index files. +By default, the code will try to extra RA and DEC from the header, but the +user can manually override those values using `--ra` and `--dec`. The default +search radius is 2 degrees but this can also be overridden by the user. + +## Support +[API documentation]() +[Create an issue](https://issues.lco.global/) diff --git a/gaia_astrometry_index_files/csv.py b/gaia_astrometry_index_files/csv.py new file mode 100644 index 0000000..c1bb991 --- /dev/null +++ b/gaia_astrometry_index_files/csv.py @@ -0,0 +1,18 @@ +import gzip +import logging +logger = logging.getLogger(__name__) + + +def parse_csv(catalog_file): + with gzip.open(catalog_file, 'rt') as csv_file: + columns = csv_file.readline().split(',') + indices = [columns.index(i) for i in ['ra', 'ra_error', 'dec', 'dec_error', 'phot_g_mean_flux', 'phot_g_mean_flux_error']] + logger.info('Making rows variable') + rows = [parse_row(line, indices) for line in csv_file] + return rows + + +def parse_row(line, indicies): + # This only works for GAIA DR2. We could get the headers from + values = line.split(',') + return (values[index] for index in indicies) diff --git a/gaia_astrometry_index_files/dbs.py b/gaia_astrometry_index_files/dbs.py index d5df8e8..abb59c2 100644 --- a/gaia_astrometry_index_files/dbs.py +++ b/gaia_astrometry_index_files/dbs.py @@ -1,9 +1,9 @@ -import gzip import sqlite3 from sqlalchemy import Column, Integer, Float, create_engine, pool from sqlalchemy.ext.declarative import declarative_base from sqlalchemy.orm import sessionmaker import logging +from gaia_astrometry_index_files.csv import parse_csv logger = logging.getLogger(__name__) @@ -13,8 +13,8 @@ class Star(Base): __tablename__ = 'sources' id = Column(Integer, primary_key=True, autoincrement=True) - ra = Column(Float) #, index=True - dec = Column(Float) # , index=True + ra = Column(Float) + dec = Column(Float) ra_error = Column(Float) dec_error = Column(Float) g_flux = Column(Float) @@ -74,18 +74,10 @@ def make_gaia_db(catalogs): for catalog_file in catalogs: logger.info('Adding {filename} to db.'.format(filename=catalog_file)) cursor.execute('begin transaction') - with gzip.open(catalog_file, 'rt') as csv_file: - columns = csv_file.readline() - logger.info('Making rows variable') - rows = [parse_row(line) for line in csv_file] + rows = parse_csv(catalog_file) logger.info('Adding records to db') cursor.executemany('INSERT INTO sources (ra, ra_error, dec, dec_error, g_flux, g_flux_error) values(?, ?, ?, ?, ?, ?)', rows) cursor.execute('COMMIT') cursor.execute('CREATE VIRTUAL TABLE if not exists positions using rtree(id, ramin, ramax, decmin, decmax);') cursor.execute('insert into positions (id, ramin, ramax, decmin, decmax) select id, ra, ra, dec, dec from sources;') cursor.close() - - -def parse_row(line): - values = line.split(',') - return values[5], values[6], values[7], values[8], values[47], values[48] \ No newline at end of file diff --git a/gaia_astrometry_index_files/main.py b/gaia_astrometry_index_files/index_file_creator.py similarity index 85% rename from gaia_astrometry_index_files/main.py rename to gaia_astrometry_index_files/index_file_creator.py index 53209c7..ae307be 100644 --- a/gaia_astrometry_index_files/main.py +++ b/gaia_astrometry_index_files/index_file_creator.py @@ -108,8 +108,8 @@ def parse_catalog(filename): return catalog_data -def query_sources(db_address, ramin, ramax, decmin, decmax): - connection = sqlite3.connect('gaia.db') +def query_sources(ramin, ramax, decmin, decmax, db_address='gaia.db'): + connection = sqlite3.connect(db_address) cursor = connection.cursor() sql_command = 'select sources.g_flux, sources.ra, sources.dec from sources, positions ' \ 'where positions.ramin >= {ramin} and positions.ramax <= {ramax} ' \ @@ -143,25 +143,25 @@ def get_sources_in_healpixel(healpixel, db_address): if decmin < -90.0 or decmax > 90.0: ramin = 0.0 ramax = 360.0 - sources_in_healpix = query_sources(db_address, ramin, ramax, decmin, decmax) + sources_in_healpix = query_sources(ramin, ramax, decmin, decmax) elif ramin < 0.0: - sources_in_healpix = query_sources(db_address, ramin, ramax, decmin, decmax) + sources_in_healpix = query_sources(ramin, ramax, decmin, decmax) ramin = 360.0 - delta_ra ramax = 360.0 - sources_from_wrap = query_sources(db_address, ramin, ramax, decmin, decmax) + sources_from_wrap = query_sources(ramin, ramax, decmin, decmax) for key in sources_in_healpix: sources_in_healpix[key] += sources_from_wrap[key] elif ramax > 360.0: - sources_in_healpix = query_sources(db_address, ramin, ramax, decmin, decmax) + sources_in_healpix = query_sources(ramin, ramax, decmin, decmax) ramin = 0.0 ramax = delta_ra - sources_from_wrap = query_sources(db_address, ramin, ramax, decmin, decmax) + sources_from_wrap = query_sources(ramin, ramax, decmin, decmax) for key in sources_in_healpix: sources_in_healpix[key] += sources_from_wrap[key] else: - sources_in_healpix = query_sources(db_address, ramin, ramax, decmin, decmax) + sources_in_healpix = query_sources(ramin, ramax, decmin, decmax) return sources_in_healpix @@ -287,36 +287,3 @@ def make_unique_id(scale): unique_id = unique_id.replace('+', '1') unique_id = unique_id.replace('-', '0') return unique_id - - -def plot_source_density(): - source_density = np.zeros((1000, 1000)) - merged_catalogs = glob('gaia-fullcatalog*.fits') - for catalog in merged_catalogs: - print(catalog) - # Read in the file - data = fits.getdata(catalog, 1) - # Convert the decs to thetas: theta = 90 - dec - theta = 90 - data['dec'] - # calculate z: z = cos(theta) - z = np.cos(np.deg2rad(theta)) - # if abs(z) < 2/3: - # xs = ra (in radians) - xs = np.deg2rad(data['ra'].copy()) - # ys = (3 pi / 8) z - ys = 3.0 * np.pi / 8.0 * z - # else: - # sigma = 2 - (3 (1 - abs(z)))^1/2 - sigma = 2 - (3 * (1.0 - np.abs(z)))**0.5 - # sigma(-z) = -sigma(z) - sigma[z < 0] = -sigma[z < 0] - # phi_l = ra % pi/2 - phi_l = np.deg2rad(data['ra']) % (np.pi / 2.0) - # xs = ra - (abs(sigma(z)) - 1)(phi_l - pi / 4) - polar_cap = np.abs(z) > (2.0 / 3.0) - xs[polar_cap] = np.deg2rad(data['ra'][polar_cap]) - (np.abs(sigma[polar_cap]) - 1.0)*(phi_l[polar_cap] - (np.pi / 4.0)) - # ys = pi / 4 * sigma(z) - ys[polar_cap] = np.pi / 4.0 * sigma[polar_cap] - source_density += np.histogram2d(xs, ys, bins=(1000, 1000), - range=[[0.0, 2.0 * np.pi], [-np.pi / 2.0, np.pi / 2.0]])[0] - return source_density diff --git a/gaia_astrometry_index_files/plots.py b/gaia_astrometry_index_files/plots.py new file mode 100644 index 0000000..168a857 --- /dev/null +++ b/gaia_astrometry_index_files/plots.py @@ -0,0 +1,37 @@ +from glob import glob + +import numpy as np +from astropy.io import fits + + +def plot_source_density(): + source_density = np.zeros((1000, 1000)) + merged_catalogs = glob('gaia-fullcatalog*.fits') + for catalog in merged_catalogs: + print(catalog) + # Read in the file + data = fits.getdata(catalog, 1) + # Convert the decs to thetas: theta = 90 - dec + theta = 90 - data['dec'] + # calculate z: z = cos(theta) + z = np.cos(np.deg2rad(theta)) + # if abs(z) < 2/3: + # xs = ra (in radians) + xs = np.deg2rad(data['ra'].copy()) + # ys = (3 pi / 8) z + ys = 3.0 * np.pi / 8.0 * z + # else: + # sigma = 2 - (3 (1 - abs(z)))^1/2 + sigma = 2 - (3 * (1.0 - np.abs(z)))**0.5 + # sigma(-z) = -sigma(z) + sigma[z < 0] = -sigma[z < 0] + # phi_l = ra % pi/2 + phi_l = np.deg2rad(data['ra']) % (np.pi / 2.0) + # xs = ra - (abs(sigma(z)) - 1)(phi_l - pi / 4) + polar_cap = np.abs(z) > (2.0 / 3.0) + xs[polar_cap] = np.deg2rad(data['ra'][polar_cap]) - (np.abs(sigma[polar_cap]) - 1.0)*(phi_l[polar_cap] - (np.pi / 4.0)) + # ys = pi / 4 * sigma(z) + ys[polar_cap] = np.pi / 4.0 * sigma[polar_cap] + source_density += np.histogram2d(xs, ys, bins=(1000, 1000), + range=[[0.0, 2.0 * np.pi], [-np.pi / 2.0, np.pi / 2.0]])[0] + return source_density diff --git a/gaia_astrometry_index_files/utils.py b/gaia_astrometry_index_files/utils.py index 8e5b7c5..0b532fd 100644 --- a/gaia_astrometry_index_files/utils.py +++ b/gaia_astrometry_index_files/utils.py @@ -28,4 +28,4 @@ def great_circle_distance(ra1, dec1, ra2, dec2): ra1_rad, dec1_rad = np.deg2rad([ra1, dec1]) ra2_rad, dec2_rad = np.deg2rad([ra2, dec2]) distance_rad = np.arccos(np.sin(dec1_rad) * np.sin(dec2_rad) + np.cos(dec1_rad) * np.cos(dec2_rad) * np.cos(ra2_rad - ra1_rad)) - return np.rad2deg(distance_rad) \ No newline at end of file + return np.rad2deg(distance_rad) diff --git a/setup.py b/setup.py index bbfbcc9..62ae75a 100644 --- a/setup.py +++ b/setup.py @@ -20,5 +20,5 @@ package_dir={'gaia_astrometry_index_files': 'gaia_astrometry_index_files'}, include_package_data=True, install_requires=['numpy', 'sqlalchemy', 'astrometry', 'astropy', 'lcogt_logging', 'sep', 'pyyaml'], - entry_points={'console_scripts': ['create_gaia_index_files=gaia_astrometry_index_files.main:create_index_files', + entry_points={'console_scripts': ['create_gaia_index_files=gaia_astrometry_index_files.index_file_creator:create_index_files', 'solve-astrometry-gaia=gaia_astrometry_index_files.solve:solve_frame']})