Refactored some code to make it easier to find things. Added more det…

…ail to the README.

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/)

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)

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]

gaia_astrometry_index_files/main.py → 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

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

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)
+    return np.rad2deg(distance_rad)

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']})