The MASH planeraty nebulae catalog: Advanced usage of HiPS and MOCs to explore complex regions of interest#

Stefania Amodea¹, Matthieu Baumann¹, Thomas Boch¹, Caroline Bot¹, Katarina A. Lutz¹, Manon Marchand¹.

  1. Université de Strasbourg, CNRS, Observatoire Astronomique de Strasbourg, UMR 7550, F-67000, Strasbourg, France

Thomas Boch and Caroline Bot wrote the original version of this tutorial, available on the EURO-VO tutorials page. It was presented at the workshop “Detecting the unexpected, Discovery in the Era of Astronomically Big Data”. The version here is an adaptation to jupyter notebooks by the Strasbourg astronomical Data Center (CDS) team.


This tutorial demonstrates an advanced usage of Hierarchical Progressive Surveys (HiPS) and Multi-Order Coverage (MOC) maps. Wis this tutorial, you will:

  1. learn how to find fits files

  2. build a map of the coverage of all the fits files

  3. select only the parts of this map that correspond to low-extinction regions

  4. retrieve objects from large catalogs inside the obtained map – that has a non-trivial shape and non-necessarily-connected regions

  5. combine catalogs to visualize a color-color diagram

# Standard Library
from pathlib import Path

# Astronomy tools
import astropy.units as u
from astropy.coordinates import Angle, SkyCoord, match_coordinates_sky
from import fits

# Access astronomical databases
import pyvo
from astroquery.vizier import Vizier

# Moc and HEALPix tools
import cdshealpix
import mocpy

# For plots
import matplotlib.pyplot as plt

# Data handling
import numpy as np

Step 1: Finding the images#

We want to find all Short-Red images in the Macquarie/AAO/Strasbourg Hα Planetary Galactic catalog (MASH) using the VizieR associated data service.


The VizieR service at CDS inventories astronomical catalogs published in the literature. Some of these catalogs contain data associated with publications and the tables therein. This data can be browsed and explored through the VizieR-associated data service, linked to the traditional VizieR table service. Here we look for images associated with the MASH catalog of planetary nebulae (Parker et. al. 2006-2008). The MASH fits files are cut-outs extracted from a larger Hα and Short Red survey to constitute a set of regions of interest around planetary nebulae.

To find VizieR-associated data, we use the Table Access Protocol (TAP) with the VizieR endpoint. Through the VizieR TAP endpoint, we can search for tables, their content, and information on associated data.

First, we search for the MASH catalog:

# give the address of the service, you can also directly visit the website
tap_vizier = pyvo.dal.TAPService("")

# a query that searches for all tables with the words MASH and Parker in their description
query = """
        SELECT  *  FROM tables 
        WHERE description LIKE '%MASH%Parker%'

mash_catalogues =
Table length=5
J_MNRASJ/MNRAS/412/223/table4tableThe nine MASH PNe detected and possibly detected in the PMN survey ( Bojicic I.S., Parker Q.A., Filipovic M.D., Frew D.J.)9
J_MNRASJ/MNRAS/412/223/table1tableMASH PNe detected in the NVSS ( Bojicic I.S., Parker Q.A., Filipovic M.D., Frew D.J.)201
V_combinedV/127A/mash2tableThe MASH-II Supplement (from paper II) ( Parker Q.A., Acker A., Frew D.J., Hartley M., Peyaud A.E.J., Phillipps S., Russeil D., Beaulieu S.F., Cohen M., Koppen J., Marcout J., Miszalski B., Morgan D.H., Morris R.A.H., Ochsenbein F., Pierce M.J.,)335
V_combinedV/127A/mash1tableThe MASH Catalog of Planetary Nebulae (paper I) ( Parker Q.A., Acker A., Frew D.J., Hartley M., Peyaud A.E.J., Phillipps S., Russeil D., Beaulieu S.F., Cohen M., Koppen J., Marcout J., Miszalski B., Morgan D.H., Morris R.A.H., Ochsenbein F., Pierce M.J.,)903
J_MNRASJ/MNRAS/412/223/mpgs2tableMASH PNe detected in the MPGS-2 (Cat. VIII/82) ( Bojicic I.S., Parker Q.A., Filipovic M.D., Frew D.J.)81

In this tutorial, we are interested in the tables belonging to the catalogs V/127A. This includes tables V/127A/mash1 and V/127A/mash2. To have a look at the content of these tables, we look at their first lines like this:

query = """
        SELECT TOP 5 * FROM \"V/127A/mash1\" 
mash1_head =
Table length=5
29LG234.7-02.2PHR0724-2021111.05458333333331-20.36361111111111234.7045-2.2774134.554.0ASA2452672.0HA18201HA842Large, very faint diffuse arcuate nebula; has [NII]~2xH-alpha, nothing in blue1029img_haimg_srfits
2PG227.3-12.0PHR0633-180898.35374999999999-18.13972222222222227.3207-12.028917.015.0EaSA2452672.0HA18191HA926Very faint, partial arcuate nebula also observed M1 060100; [NII]~0.8xH-alpha, strong [SII], only weak H-beta in blue - inconclusive1002img_haimg_srfits
16TG227.2-03.4PHR0705-1419106.41041666666665-14.318055555555553227.2852-3.402915.015.0ESA2452668.0HA18244HA1017Small, circular PN around a faint central star; also observed M1 040100; [NII]~0.7 H-alpha, [OIII]>>H-beta1016img_haimg_srfits
5LG223.6-06.8PHR0646-1235101.60583333333332-12.598888888888887223.6338-6.803540.037.0ESA2453788.0HA18194HA1016Slightly oval very faint PN candidate - has [OIII] and H-alpha1005img_haimg_srfits
10PG224.3-05.5PHR0652-1240103.08458333333331-12.67611111111111224.3504-5.5463187.0180.0ISA2452670.0HA18244HA1017Faint, extended S-shaped emission nebula, possible evolved PN, also observed M1 080100; has [NII]~0.8Ha, [OIII], strong [SII], [OIII]>H-beta1010img_haimg_srfits

As you can see, the last column of this table is called AssocData and contains the entry fits. If you look at this table on the VizieR web interface, you can download the associated fits file. Within this notebook, we query the obscore database to get the URLs to the fits file. Using the module, we can open the fits files from their URLs.

obs_tap_vizier = pyvo.dal.TAPService(
query = """
        SELECT TOP 5  *  FROM obscore 
        WHERE obs_collection='V/127A'   
mash_fits =
Table length=5
610560application/fits IV/127A1032_sr.fitsivo://CDS.VizieR/V/127A?res=1032_sr.fits864975549779410945obs.imageNotSet-13.7843898548012030.0713021682865721111.43275777455307Polygon ICRS 111.46664940645094 -13.81781914701963 111.39832407190251 -13.817324215795937 111.3988758446347 -13.750955924906409 111.46718176614709 -13.7514507063586090.948991575097523----51226.0--
740160application/fits IV/127A1034_ha.fitsivo://CDS.VizieR/V/127A?res=1034_ha.fits864975549779410951obs.imageNotSet-17.120492976226730.07123944496920925111.78451610217952Polygon ICRS 111.8191479048229 -17.153711083791936 111.74974303867411 -17.15359964722061 111.74989666803143 -17.087268981666405 111.8192767884711 -17.0873803700202520.9485038734315203----52338.0--
596160application/fits IV/127A1034_sr.fitsivo://CDS.VizieR/V/127A?res=1034_sr.fits864975549779410955obs.imageNotSet-17.1205339692782450.07124314880377468111.7845314034197Polygon ICRS 111.8191637002453 -17.153753988164176 111.74975722819107 -17.153640656320004 111.7499114759699 -17.08730806317958 111.81929320101966 -17.0874213466143270.9485285228101485----50871.0--
731520application/fits IV/127A1036_ha.fitsivo://CDS.VizieR/V/127A?res=1036_ha.fits864975549779410959obs.imageNotSet-21.8620994279805560.07130134708963763112.69528726569132Polygon ICRS 112.73146814532774 -21.89490077932626 112.65997554014551 -21.895645624610193 112.659123005937 -21.829290184523508 112.73058240266795 -21.8285457137899180.9488088411876365----51163.0--
590400application/fits IV/127A1005_sr.fitsivo://CDS.VizieR/V/127A?res=1005_sr.fits864975549779410960obs.imageNotSet-12.598781947184450.06695119600443825101.605805625331Polygon ICRS 101.64022494465036 -12.631531553985315 101.57226287333897 -12.63235954506209 101.57139509987323 -12.56602794013158 101.63933959085792 -12.5652001703247580.9485060738987898----51156.0--

As you can see, the result from this query provides information on the fits files associated with the MASH catalogs. In particular, the column access_url provides the location of the data. To get the first image we could do:

image =['access_url'][0])

and then work on the image, plot it or save it to our machine. However, downloading all the data takes quite some time. For this tutorial, we prepared a subsample of 335 of these Short Red images that will run promptly but we encourage you to try accessing the full Short Red sample on your own later. The subsample is available either at or in the Data Folder of this repository.

Step 2: Create a MOC of the MASH images#

The multi-order coverage (MOC) map of a set of images represents their sky coverage. MOCs can describe arbitrary zones in the sky which do not need to be connected. You’ll see that the union or intersection of two MOCs requires few time and computational effort. Catalogs can also be filtered by MOCs.

Here we want to use the fits files just downloaded to create a MOC map corresponding to the coverage of the MASH images.

Organising data#

# Where to find fits images downloaded from the archive above
datadir = Path("Data/MASH_Sample/")
datadir.mkdir(parents=True, exist_ok=True)
# Make new directory where to save MOCs
outdir = Path("Data/Result_moc/")
outdir.mkdir(parents=True, exist_ok=True)

In most cases, we could ignore the next cell. However, some possible deprecated keywords in the fits header would hamper the MOC creation and would cause errors in the underlying astropy.wcs.WCS module. This is why we rewrite the headers of the fits files so that they only contain the useful keywords needed to define the coordinate frame correctly before using mocpy.

Create the MOC#

Here we can create the MOC of the MASH images with the MOCPy module. This can take a few seconds.

# Get a list of the fits files and create the MOC
mash_file_list = datadir.glob(
)  # glob allows to find all paths that end with _sr_modified.fits in datadir
moc_mash = mocpy.MOC.from_fits_images(mash_file_list, max_norder=15) / "mash_moc.fits"), overwrite=True)
# this takes a bit of time because a lot of files are opened and processed to extract the moc
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]

Plot the MOC#

We can chose between two of the MOCpy methods to plot the MOC

# A one-liner for a very fast vizualisation
# With a bit more control over the output, the MOC.wcs method
fig = plt.figure(figsize=(10, 7))
wcs = moc_mash.wcs(fig)
ax = fig.add_subplot(projection=wcs)
ax.grid(True)  # noqa: FBT003
moc_mash.fill(ax, wcs, color="purple", alpha=0.5)
moc_mash.border(ax, wcs, color="purple")

You can see how the MOC has arbitrary shapes and not all regions are connected.

And for more control over the plot parameters, there is also the mocpy.WCS method (!)

Step 3: Load an archival extinction map and create the MOC of the low extinction regions#

Different works (e.g. Schlegel et al. 1998, Schlafly & Finkbeiner 2011, Green et al. 2015…) have created extinction maps of the sky that are publicly available. Some of these maps are all-sky maps, while others have higher resolutions, or come from different methods… They can be found in HEALPix format (among others) on the Legacy Archive for Microwave Background Data Analysis (LAMBDA) website or on the Analysis Center for Extended Data (CADE) website.

For this tutorial, we will download the well-known all-sky extinction map from Schlegel et al. from the LAMBDA website and define the low extinction area for which \(0 < E(B-V) < 0.5\) as a MOC. It has an information page.

The map is available here: and we save it to disc.

ext_map =
ext_map.writeto(outdir / "Schlegel_extinction_map.fits", overwrite=True)

We are only interested in regions with low extinction. So we aim to get a MOC of all regions where the extinction values from the Schlegel et al. map are between 0 and 0.5mag. The extinction map we got from the NASA webpage is in the HEALPix format. This is an efficient presentation of all-sky maps. The HEALPix tesselation is also used by MOCs. So to get the MOC from the extinction map, we do the following.

First, we check the coordinate system in the map header. We will need to convert to equatorial coordinates, change the projection of the map, and set the order (i.e. resolution) of the map.

# open the downloaded FITS file
hdul = / "Schlegel_extinction_map.fits")
hdr = hdul[0].header
SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                   32 / number of bits per data pixel                  
NAXIS   =                    0 / number of data axes                            
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
DATE    = '2003-02-05T00:00:00' /file creation date (YYYY-MM-DDThh:mm:ss UT)    
OBJECT  = 'ALL-SKY '           / Portion of sky given                           
COMMENT   This file contains an all-sky Galactic reddening map, E(B-V), based on
COMMENT   the derived reddening maps of Schlegel, Finkbeiner and Davis (1998).  
COMMENT   Software and data files downloaded from their website were used to    
COMMENT   interpolate their high resolution dust maps onto pixel centers        
COMMENT   appropriate for a HEALPix Nside=512 projection in Galactic            
COMMENT   coordinates. This file is distributed and maintained by LAMBDA.       
REFERENC= 'Legacy Archive for Microwave Background Data Analysis (LAMBDA)      '
REFERENC= '                              '
REFERENC= 'Maps of Dust Infrared Emission for Use in Estimation of Reddening an'
REFERENC= ' Cosmic Microwave Background Radiation Foregrounds',                '
REFERENC= ' Schlegel, Finkbeiner & Davis 1998 ApJ 500, 525                     '
REFERENC= 'Berkeley mirror site for SFD98 data:'
REFERENC= 'Princeton mirror site for SFD98 data:                               '
REFERENC= '                           http://astro/'
REFERENC= 'HEALPix Home Page:              '
RESOLUTN=                    9 / Resolution index                               
SKYCOORD= 'Galactic'           / Coordinate system                              
PIXTYPE = 'HEALPIX '           / Pixel algorithm                                
ORDERING= 'NESTED  '           / Ordering scheme                                
NSIDE   =                  512 / Resolution parameter                           
NPIX    =              3145728 / # of pixels                                    
FIRSTPIX=                    0 / First pixel (0 based)                          
LASTPIX =              3145727 / Last pixel (0 based)                           
FITS_rec([( 9.625492, 1.), (46.090515, 1.), ( 8.18071 , 1.), ...,
          (15.149189, 1.), (14.107367, 1.), (15.463686, 1.)],
         dtype=(numpy.record, [('TEMPERATURE', '>f4'), ('N_OBS', '>f4')]))

The data field here has a specific shape. It contains tuples for which the first value is the extinction (named ‘TEMPERATURE’) and the second one is the number of observations of the value (you can check that it is 1 everywhere).

extinction_values = hdul[1].data["TEMPERATURE"]

Let’s extract the information about the number of sides and the order of the healpix map from the header of the fits file

nside = hdul[0].header["NSIDE"]
norder = hdul[0].header["RESOLUTN"]

The header allows to see that this map is in galactic coordinates. We will need to convert this into equatorial coordinates to compare with our other maps.

# Creation of an HEALpix grid at order 9 in nested ordering
healpix_index = np.arange(12 * 4**norder, dtype=np.uint64)
    f"We can check the the NPIX value corresponds to the one in the header here: {len(healpix_index)}",
We can check the the NPIX value corresponds to the one in the header here: 3145728
# Get the coordinates of the centers of these healpix cells
center_coordinates_in_equatorial = cdshealpix.healpix_to_skycoord(
<SkyCoord (ICRS): (ra, dec) in deg
    [( 45.        ,  0.0746039 ), ( 45.08789062,  0.14920793),
     ( 44.91210938,  0.14920793), ..., (315.08789062, -0.14920793),
     (314.91210938, -0.14920793), (315.        , -0.0746039 )]>
# Convert this into galactic coordinates
center_coordinates_in_galactic = center_coordinates_in_equatorial.galactic
<SkyCoord (Galactic): (l, b) in deg
    [(176.8796283 , -48.85086427), (176.89078038, -48.7358142 ),
     (176.70525363, -48.86216423), ..., ( 48.82487228, -28.4122831 ),
     ( 48.7216889 , -28.26178141), ( 48.84578935, -28.29847774)]>
# Calculate the bilinear interpolation we must apply to each healpix cell to get the values in the other coordinate system
healpix, weights = cdshealpix.bilinear_interpolation(
# Apply the interpolation
ext_map_equatorial_nested = (extinction_values[] *
array([0.08981742, 0.0991632 , 0.08249644, ..., 0.08323811, 0.08352184,
       0.0820533 ])

Next we declare which pixel we want to use, let’s take all pixels with an extinction lower than 0.5:

low_extinction_index = np.where(ext_map_equatorial_nested < 0.5)[0]
    f"The low extinction criteria keeps {round((len(low_extinction_index)/ len(extinction_values)*100), 2)}% of the sky map",
The low extinction criteria keeps 86.74% of the sky map

And let’s create a MOC from this criteria

moc_low_extinction = mocpy.MOC.from_healpix_cells(

Step 4: Find out which regions are covered by the MASH short-red images in the low extinction regions defined above#

To find out the sky regions of the MASH sample that are at low extinction, we build the intersection of the two MOCs.

moc_intersection = moc_low_extinction.intersection(moc_mash)
# Once the intersection is bluit, we can for example print the sky fraction :
    f"The intersection of the two MOCs covers {round(moc_intersection.sky_fraction * 100, 4)}% of the sky",
The intersection of the two MOCs covers 0.4778% of the sky

Now we can visualize the coverage of the two MOCs and their intersection. The grey area is where the extinction is low. The blue one is the MASH coverage. The tiny red dots show the MASH coverage in low extinction regions.

fig = plt.figure(111, figsize=(10, 8))

with mocpy.WCS(
    fov=200 * u.deg,
    center=SkyCoord(200, -20, unit="deg", frame="icrs"),
) as wcs:
    ax = fig.add_subplot(1, 1, 1, projection=wcs)

        label="low extinction",
        label="MASH coverage",
        label="MASH in low extinction regions",
    # Sets labels
    # Sets ticks
    lon, lat = ax.coords[0], ax.coords[1]
    lon.set_ticks(spacing=2 * u.hourangle)
plt.legend(bbox_to_anchor=(1, 1))
<matplotlib.legend.Legend at 0x7f51482d8090>

Step 5: Query the 2MASS and Gaia Catalogues by MOC#

Without the usage of MOC, querying for sources in the low extinction regions covered by the MASH subsample would be tedious or even impossible. Indeed, one would need to load the whole catalog and make selections which would not be possible given the size of some catalogs. Alternatively, one would need to query the catalog field by field, which would take time and several queries. Instead, here we will use the power of MOC files to query large catalogs directly in the covered regions only. We will use coverages of the low extinction and MASH regions to query for sources from the Gaia and 2MASS surveys in these highly non-continuous and non-trivial shape areas.

First, let’s see which Gaia and 2MASS catalogs are available on VizieR. We could, as above, use the TAP endpoint of VizieR. But we show below the Vizier module in the astroquery package.

catalog_list_twomass = Vizier.find_catalogs("Cutri")
for k, v in catalog_list_twomass.items():
    print(k, ": ", v.description)
II/126 :  IRAS Serendipitous Survey Catalog (IPAC 1986)
II/241 :  2MASS Catalog Incremental Data Release (IPAC/UMass, 2000)
II/246 :  2MASS All-Sky Catalog of Point Sources (Cutri+ 2003)
II/281 :  2MASS 6X Point Source Working Database / Catalog (Cutri+ 2006)
II/307 :  WISE Preliminary Data Release (Cutri+ 2011)
II/311 :  WISE All-Sky Data Release (Cutri+ 2012)
II/328 :  AllWISE Data Release (Cutri+ 2013)
II/365 :  The CatWISE2020 catalog (updated version 28-Jan-2021) (Marocco+, 2021)
VII/233 :  The 2MASS Extended sources (IPAC/UMass, 2003-2006)
J/ApJ/560/566 :  K-band galaxy luminosity function from 2MASS (Kochanek+, 2001)
J/ApJ/564/421 :  Spectra of T dwarfs. I. (Burgasser+, 2002)
J/ApJ/569/23 :  Optical polarisation of 2MASS QSOs (Smith+, 2002)
J/ApJ/635/214 :  Chandra X-ray sources and NIR identifications (Ebisawa+, 2005)
J/ApJ/713/330 :  Spitzer observations of major-merger galaxies (Xu+, 2010)
J/ApJ/719/550 :  Deep NIR imaging of {rho} Oph cloud core (Marsh+, 2010)
J/ApJ/741/68 :  Main Belt asteroids with WISE/NEOWISE. I. (Masiero+, 2011)
J/ApJ/742/40 :  Jovian Trojans asteroids with WISE/NEOWISE (Grav+, 2011)
J/ApJ/743/156 :  NEOWISE observations of NEOs: preliminary results (Mainzer+, 2011)
J/ApJ/744/197 :  WISE/NEOWISE observations of Hilda asteroids (Grav+, 2012)
J/ApJ/759/L8 :  WISE/NEOWISE observations of main belt asteroids (Masiero+, 2012)
J/ApJ/760/L12 :  WISE/NEOWISE NEOs preliminary thermal fits (Mainzer+, 2012)
J/ApJ/780/92 :  Light curves of the RR Lyr SDSS J015450.17+001500.7 (Szabo+, 2014)
J/ApJ/783/122 :  AllWISE motion survey (Kirkpatrick+, 2014)
J/ApJ/784/110 :  NEOWISE observations of 105 near-Earth objects (Mainzer+, 2014)
J/ApJ/792/30 :  NEOWISE magnitudes for near-Earth objects (Mainzer+, 2014)
J/ApJ/805/90 :  WISE ELIRGs and comparison with QSOs (Tsai+, 2015)
J/ApJ/814/117 :  NEOWISE Reactivation mission: 1st yr data (Nugent+, 2015)
J/ApJ/874/82 :  Follow-up photometry & spectroscopy of PTF14jg (Hillenbrand+, 2019)
J/ApJS/95/1 :  Atlas of Quasar Energy Distributions (Elvis+ 1994)
J/ApJS/172/663 :  Infrared observations of the Pleiades (Stauffer+, 2007)
J/ApJS/175/191 :  Variables from 2MASS calibration fields (Plavchan+, 2008)
J/ApJS/190/100 :  NIR proper motion survey using 2MASS (Kirkpatrick+, 2010)
J/ApJS/199/26 :  The 2MASS Redshift Survey (2MRS) (Huchra+, 2012)
J/ApJS/224/36 :  The AllWISE motion survey (AllWISE2) (Kirkpatrick+, 2016)
J/ApJS/234/23 :  The WISE AGN candidates catalogs (Assef+, 2018)
J/AJ/112/62 :  Quasar absorption-line systems (Tanner+ 1996)
J/AJ/125/2521 :  2MASS6x survey of the Lockman Hole (Beichman+, 2003)
J/AJ/126/63 :  Host galaxies of 2MASS-QSOs with z<=3 (Hutchings+, 2003)
J/AJ/127/646 :  Unbiased census of AGN in 2MASS (Francis+, 2004)
J/AJ/135/2245 :  Absolute spectrum of the Sun and Vega for 0.2-30um (Rieke+, 2008)
J/AJ/144/148 :  Infrared photometry of brown dwarf and Hyper-LIRG (Griffith+, 2012)
J/AJ/152/63 :  NEOWISE reactivation mission: 2nd yr data (Nugent+, 2016)
J/AJ/154/53 :  WISE/NEOWISE observations of comets (Bauer+, 2017)
J/AJ/154/168 :  NEOWISE: thermal model fits for NEOs and MBAs (Masiero+, 2017)
J/AJ/156/60 :  Thermal model fits for short-arc NEOs with NEOWISE (Masiero+, 2018)
J/PASP/113/10 :  Sub-mJy radio sources complete sample (Masci+, 2001)
catalog_list_gaia = Vizier.find_catalogs("Gaia DR2", max_catalogs=1000)
for k, v in catalog_list_gaia.items():
    print(k, ": ", v.description)
WARNING: UnitsWarning: Unit 'x' not supported by the VOUnit standard.  [astropy.units.format.vounit]
WARNING: UnitsWarning: Unit 'Earth' not supported by the VOUnit standard.  [astropy.units.format.vounit]
WARNING: UnitsWarning: Unit 'Sun' not supported by the VOUnit standard. Did you mean uN? [astropy.units.format.vounit]
WARNING: UnitsWarning: Unit 'kau' not supported by the VOUnit standard. Did you mean kA, ka, kadu or ku? [astropy.units.format.vounit]
I/345 :  Gaia DR2 (Gaia Collaboration, 2018)
I/347 :  Distances to 1.33 billion stars in Gaia DR2 (Bailer-Jones+, 2018)
I/349 :  StarHorse, Gaia DR2 photo-astrometric distances (Anders+, 2019)
II/360 :  Gaia DR2 x AllWISE catalogue (Marton+, 2019)
IV/35 :  Gaia DR2-WISE Galactic Plane Matches (Wilson+, 2018)
VI/156 :  M-dwarf Lum-Temp-Radius relationships (Morrell+, 2019)
VII/285 :  Gaia DR2 quasar and galaxy classification (Bailer-Jones+, 2019)
J/ApJ/862/138 :  BANYAN. XIII. Nearby young assoc. with Gaia DR2 (Gagne+, 2018)
J/ApJ/863/89 :  Gaia DR2 PMs of stars in ultra-faint MW satellites (Simon, 2018)
J/ApJ/866/99 :  Revised radii of KIC stars & planets using Gaia DR2 (Berger+, 2018)
J/ApJ/867/151 :  YSOs in the Gould Belt regions with Gaia-DR2 data (Dzib+, 2018)
J/ApJ/868/70 :  Hot subdwarf stars from Gaia DR2 and LAMOST DR5 (Lei+, 2018)
J/ApJ/870/32 :  Kinematics in young star clusters & associations (Kuhn+, 2019)
J/ApJ/872/85 :  Detached eclipsing binaries with Gaia parallaxes (Graczyk+, 2019)
J/ApJ/873/132 :  ICRF3 radio-loud quasars with Gaia DR2 data (Makarov+, 2019)
J/ApJ/874/138 :  Gaia and LAMOST DR4 M giant members of Sgr stream (Li+, 2019)
J/ApJ/875/77 :  Proper motions of MW satellites with Gaia & DES (Pace+, 2019)
J/ApJ/877/29 :  Opt & NIR SMARTS/ANDICAM photometry for DQ Tau (Muzerolle+, 2019)
J/ApJ/878/111 :  Members in Serpens Molecular Cloud with Gaia DR2 (Herczeg+, 2019)
J/ApJ/881/7 :  Hot subdwarf stars from LAMOST DR5 & Gaia DR2 (Luo+, 2019)
J/ApJ/881/135 :  Hot subdwarf stars from Gaia DR2 and LAMOST DR5. II. (Lei+, 2019)
J/ApJ/883/8 :  Cand. young OB stars from GALEX & Gaia DR2 (Casetti-Dinescu+, 2019)
J/ApJ/884/6 :  Known members of Orion A with Gaia DR2 data (McBride+, 2019)
J/ApJ/886/100 :  High-mass white dwarfs in Gaia DR2 (Cheng+, 2019)
J/ApJ/887/126 :  Follow-up obs. of the PSJ0147+4630 lens system (Goicoechea+, 2019)
J/ApJ/888/73 :  Varstrometry of SDSS quasars with Gaia DR2 data (Hwang+, 2020)
J/ApJ/889/99 :  Gaia DR2 Blanco 1 member candidates (Zhang+, 2020)
J/ApJ/889/117 :  Hot subdwarf stars from Gaia DR2 & LAMOST DR6+7. I. (Lei+, 2020)
J/ApJ/889/176 :  Very low-mass binaries with Gaia DR2 data (Faherty+, 2020)
J/ApJ/898/64 :  Hot subdwarf stars from Gaia with LAMOST sp. II. RVs (Luo+, 2020)
J/ApJ/898/173 :  Solar analog rotations from Kepler & Gaia (Do Nascimento+, 2020)
J/ApJ/899/83 :  723 Gaia DR2 White dwarfs cand. in Local Galactic Halo (Kim+, 2020)
J/ApJ/910/46 :  Galactic anticenter substructure stars from LAMOST (Li+, 2021)
J/ApJ/913/120 :  4FGL blazars with Gaia DR2 data (Zeng+, 2021)
J/ApJ/914/16 :  NGC 5128 GCs from PISCes, Gaia DR2 and NSC (Hughes+, 2021)
J/ApJ/914/122 :  Stellar polarimetry with Gaia DR2 data (Doi+, 2021)
J/ApJ/914/123 :  Gaia DR2 and EDR3 stars with sp. follow-up (Ibata+, 2021)
J/ApJ/915/103 :  Oscillator strength ref. for Indus_0 & Indus_13 (Hansen+, 2021)
J/ApJ/921/42 :  Gaia GraL. VI. Quadruply imaged lensed QSOs (Stern+, 2021)
J/ApJS/245/32 :  Newly identified star clusters in Gaia DR2 (Liu+, 2019)
J/ApJS/246/4 :  Catalog of ultrawide binary stars from Gaia DR2 (Tian+, 2020)
J/ApJS/252/3 :  High-velocity stars in the Gal. halo from LAMOST & Gaia (Li+, 2021)
J/ApJS/253/58 :  Ages of field stars from white dwarf comp. in Gaia (Qiu+, 2021)
J/ApJS/254/20 :  Stellar groups in Taurus field from Gaia DR2 & LAMOST (Liu+, 2021)
J/ApJS/256/28 :  Hot subdwarf stars with Gaia DR2 and LAMOST DR7 data (Luo+, 2021)
J/ApJS/260/8 :  541 new open cluster candidates (He+, 2022)
J/A+A/616/A7 :  Gaia DR2 radial velocity standard stars catalog (Soubiran+, 2018)
J/A+A/616/A12 :  Gaia DR2 sources in GC and dSph (Gaia Collaboration+, 2018)
J/A+A/616/A37 :  Close encounters to the Sun in Gaia DR2 (Bailer-Jones+, 2018)
J/A+A/616/L2 :  Planetary Nebulae distances in Gaia DR2 (Kimeswenger+, 2018)
J/A+A/616/L15 :  Parallaxes and Proper Motions of OB stars (Xu+, 2018)
J/A+A/617/A135 :  20 years of photometric microlensing (Mustill+, 2018)
J/A+A/618/A44 :  Predicted microlensing events from Gaia DR2 (Bramich, 2018)
J/A+A/618/A56 :  Gaia GraL. II. Known multiply imaged quasars (Ducourant+, 2018)
J/A+A/618/A59 :  Gaia DR2 confirmed new nearby open clusters (Castro-Ginard+, 2018)
J/A+A/618/A93 :  Gaia DR2 open clusters in the Milky Way (Cantat-Gaudin+, 2018)
J/A+A/619/A8 :  Cepheid period-luminosity-metallicity relation (Groenewegen, 2018)
J/A+A/619/A106 :  3D shape of Orion A from Gaia DR2 (Grossschedl+, 2018)
J/A+A/619/A155 :  Open cluster kinematics with Gaia DR2 (Soubiran+, 2018)
J/A+A/619/A180 :  Gaia DR2 photometric sensitivity curves (Maiz Apellaniz+, 2018)
J/A+A/619/L8 :  Ultra-cool dwarfs candidates in Gaia DR2 (Reyle, 2018)
J/A+A/620/A91 :  New asteroid models (Durech+, 2018)
J/A+A/620/A128 :  Gaia DR2 study of Herbig Ae/Be stars (Vioque+, 2018)
J/A+A/620/A141 :  Physical properties of AM CVn stars (Ramsay+, 2018)
J/A+A/620/A155 :  Gaia proper motions of 7 UFD galaxies (Massari+, 2018)
J/A+A/620/A172 :  Solar neighbourhood young stars 3D mapping (Zari+, 2018)
J/A+A/621/A38 :  Gaia catalogue of hot subluminous stars (Geier+, 2019)
J/A+A/621/A48 :  Gaia-DR2 extended kinematical maps. I. (Lopez-Corredoira+, 2019)
J/A+A/621/L2 :  Hyades tidal tails revealed by Gaia DR2 (Roeser+, 2019)
J/A+A/621/L3 :  Hyades tidal tails with Gaia DR2 (Meingast+, 2019)
J/A+A/622/A60 :  Gaia DR2 misclassified RR Lyrae list (Clementini+, 2019)
J/A+A/622/A165 :  Gaia GraL. III. New lensed systems (Delchambre+, 2019)
J/A+A/622/L13 :  Stellar stream in Gaia DR2 discovery (Meingast+, 2019)
J/A+A/623/A22 :  IC 4996 Vilnius phot. and Gaia DR2 astrometry (Straizys+, 2019)
J/A+A/623/A25 :  Photometric and astrometric study of NGC 6530 (Damiani+, 2019)
J/A+A/623/A72 :  Binarity of Hipparcos stars from Gaia pm anomaly (Kervella+, 2019)
J/A+A/623/A108 :  Age of 269 GDR2 open clusters (Bossini+, 2019)
J/A+A/623/A110 :  Gaia DR2. Variable stars in CMD (Gaia Collaboration+, 2019)
J/A+A/623/A112 :  The Sco OB2 population from Gaia DR2 data (Damiani+, 2019)
J/A+A/623/A116 :  Galactic Cepheids and RR Lyrae multiplicity. I. (Kervella+, 2019)
J/A+A/623/A117 :  Galactic Cepheids and RR Lyrae multiplicity. II (Kervella+, 2019)
J/A+A/623/A129 :  New satellites of the LMC search (Fritz+, 2019)
J/A+A/624/A78 :  Masses and ages of 1059 HARPS-GTO stars (Delgado Mena+, 2019)
J/A+A/624/A126 :  New open clusters in Perseus direction (Cantat-Gaudin+, 2019)
J/A+A/625/A14 :  Reclassification of Cepheids in the Gaia DR2 (Ripepi+, 2019)
J/A+A/626/A80 :  Census of Rho Oph candidate members from Gaia DR2 (Canovas+, 2019)
J/A+A/627/A35 :  New open clusters in Galactic anti-centre (Castro-Ginard+, 2019)
J/A+A/627/A119 :  Extended halo of NGC 2682 (M 67) (Carrera+ 2019)
J/A+A/628/A81 :  Gaia DR2-based catalogue of 237 Ap stars (Scholz+, 2019)
J/A+A/629/A114 :  sigma Ori GTC+INT spectroscopy (Caballero+, 2019)
J/A+A/630/A119 :  Gaia DR2 distances to two clusters (Maiz Apellaniz, 2019)
J/A+A/630/A137 :  Structure and kinematics of the Taurus region (Galli+, 2019)
J/A+A/630/A150 :  Astrometric data for 211 GAPN sample (Gonzalez-Santamaria+, 2019)
J/A+A/631/A2 :  List of new asteroid models (Durech+, 2019)
J/A+A/631/A145 :  HIP stars DEC proper motions comparison (Damljanovic+, 2019)
J/A+A/633/A98 :  Gaia16aye microlensing event photometry (Wyrzykowski+, 2020)
J/A+A/633/A99 :  Gaia DR2 open clusters in the Milky Way. II (Cantat-Gaudin+, 2020)
J/A+A/633/A135 :  Solar neighbourhood carbon stars properties (Abia+, 2020)
J/A+A/634/A98 :  Corona-Australis DANCe. I. (Galli+, 2020)
J/A+A/635/A45 :  570 new open clusters in the Galactic disc (Castro-Ginard+, 2020)
J/A+A/635/L3 :  Candidate member stars of the Sagittarius stream (Antoja+, 2020)
J/A+A/636/A20 :  OGLE-III parallax events with Gaia DR2 (Wyrzykowski+, 2020)
J/A+A/637/A95 :  16 open clusters UBVI and Gaia DR2 photometry (Perren+, 2020)
J/A+A/638/A21 :  New Herbig Ae/Be and classical Be stars catalog (Vioque+, 2020)
J/A+A/638/A103 :  Central stars of planetary nebulae in Gaia DR2 (Chornay+, 2020)
J/A+A/638/A104 :  Gaia DR2 candidate RR Lyrae of Sgr stream & dwarf (Ramos+, 2020)
J/A+A/641/A79 :  FEDReD. Extinction map with 2MASS and Gaia DR2 (Hottier+, 2020)
J/A+A/642/A21 :  Astrometric parameters for 135 OB stars (Russeil+, 2020)
J/A+A/643/A148 :  Lupus DANCe. Census and 6D structure with Gaia-DR2 (Galli+, 2020)
J/A+A/644/A17 :  Astrometrically-selected QSO candidates (Heintz+, 2020)
J/A+A/645/A72 :  Radial velocities of UCAC4 721-037069 (Clavel+, 2021)
J/A+A/645/A116 :  NGC 2808, NGC 6266 and NGC 6397 Gaia DR2 sources (Kundu+, 2021)
J/A+A/646/A46 :  Chamaeleon DANCe. Stellar population with Gaia-DR2 (Galli+, 2021)
J/A+A/646/A99 :  Gaia DR2 Monoceros and ACS candidates (Ramos+, 2021)
J/A+A/646/A104 :  Improving the open cluster census. I. (Hunt+, 2021)
J/A+A/646/A183 :  Transit photometry of NGTS-14Ab (Smith+, 2021)
J/A+A/648/A44 :  Large-amplitude variables in Gaia DR2 (Mowlavi+, 2021)
J/A+A/653/A160 :  Updated radial velocities from Gaia DR2 (Seabroke+, 2021)
J/A+A/657/A131 :  Ten new systems with O stars (Maiz Apellaniz+, 2022)
J/AJ/156/94 :  APOGEE and Gaia DR2 analysis of IC 166 (Schiappacasse-Ulloa+, 2018)
J/AJ/156/139 :  Main belt asteroid shape distrib. from Gaia DR2 (Mommert+, 2018)
J/AJ/156/259 :  Robo-AO detected close binaries in Gaia DR2 (Ziegler+, 2018)
J/AJ/156/264 :  California-Kepler Survey. VII. Planet radius gap (Fulton+, 2018)
J/AJ/157/78 :  Double & multiple star systems from GaiaDR2 (Jimenez-Esteban+, 2019)
J/AJ/158/20 :  K-M stars of class I candidate RSGs in Gaia DR2 (Messineo+, 2019)
J/AJ/158/109 :  Occurrence rates of planets orbiting FGK stars (Hsu+, 2019)
J/AJ/158/155 :  SB candidates from the RAVE & Gaia DR2 surveys (Birko+, 2019)
J/AJ/159/84 :  Machine-learning regression of extinction in Gaia DR2 (Bai+, 2020)
J/AJ/159/95 :  IC 1369 Vilnius photometry and Gaia DR2 astrometry (Straizys+, 2020)
J/AJ/159/121 :  V-band light curve & RVs of the Cepheid V473Lyr (Evans+, 2020)
J/AJ/159/139 :  Multiple M dwarf stars with Robo-AO and Gaia DR2 (Lamman+, 2020)
J/AJ/159/166 :  Membership & properties of moving groups with Gaia (Ujjwal+, 2020)
J/AJ/159/259 :  Opt-NIR light curve of the quasar 3C 273 (Sobrino Figaredo+, 2020)
J/AJ/160/138 :  68 Gaia DR2 ultra-short-period planet host stars (Hamer+, 2020)
J/AJ/160/284 :  2590 binary stars in Hipparcos and Gaia DR2 (Makarov, 2020)
J/AJ/161/65 :  THYME. IV. 3 Exoplanets around TOI-451 B (Newton+, 2021)
J/AJ/161/231 :  A list of ~330000 stars Kepler missed (Wolniewicz+, 2021)
J/AJ/161/234 :  1.7 million K and M dwarfs cross-matching (Medan+, 2021)
J/AJ/161/237 :  Positions of Triton with Sheshan Station telescope (Zhang+, 2021)
J/AJ/162/181 :  RVel & Hipparcos positions of epsilon Eridani (Llop-sayson+, 2021)
J/AJ/163/219 :  Gaia DR2 stellar flybys in the Sco-Cen OB association (Ma+, 2022)
J/PASP/130/L4101 :  NGC 2112, 2477, 7789 and Col 261 members (Gao+ 2018)
J/PASP/131/D4101 :  Likely Pleiades members with Gaia DR2 (Gao, 2019)
J/MNRAS/480/467 :  Surface gravities for 15 000 Kepler stars (Pande+, 2018)
J/MNRAS/480/4505 :  A white dwarf catalogue from Gaia-DR2 (Jimenez-Esteban+, 2018)
J/MNRAS/480/5242 :  (non-)existence of five sparse open clusters (Kos+, 2018)
J/MNRAS/481/1195 :  Gaia DR2 RR Lyrae stars as standard candles (Muraveva+, 2018)
J/MNRAS/481/3887 :  Astrometry and photometry of 7 open clusters (Dias+, 2018)
J/MNRAS/482/4570 :  Gaia DR2 white dwarf candidates (Gentile Fusillo+, 2019)
J/MNRAS/482/5138 :  Galactic GC mean proper motions & velocities (Baumgardt+, 2019)
J/MNRAS/482/5222 :  DA and DB white dwarfs properties in Gaia DR2 (Tremblay+, 2019)
J/MNRAS/483/5508 :  Three new Galactic star clusters (Ferreira+, 2019)
J/MNRAS/484/1838 :  Spatial substructure of Cygnus OB2 (Berlanas+, 2019)
J/MNRAS/484/2832 :  Proper motions of Milky Way globular clusters (Vasiliev, 2019)
J/MNRAS/485/4906 :  GC number density profiles using Gaia DR2 (de Boer+, 2019)
J/MNRAS/485/5573 :  Gaia-DR2 100pc white dwarf population (Torres+, 2019)
J/MNRAS/486/2477 :  Catalogue of members of NGC 6530 (Wright+, 2019)
J/MNRAS/486/5726 :  80 young open clusters from Gaia DR2 (Dias+, 2019)
J/MNRAS/487/1462 :  The metal-rich halo tail extended in z (Fernandez-Alvar+, 2019)
J/MNRAS/487/2385 :  Distances and ages of 150 open clusters (Monteiro+, 2019)
J/MNRAS/487/2771 :  Gaia-DR2 distance to the W3 Complex (Navarete+, 2019)
J/MNRAS/487/4832 :  Period-luminosity relation of red supergiants (Chatys+, 2019)
J/MNRAS/488/2743 :  Giant Galactic dwarf satellite in Gaia DR2 (Torrealba+, 2019)
J/MNRAS/488/2892 :  Gaia DR2 extremely low-mass WD candidates (Pelisoli+, 2019)
J/MNRAS/488/3024 :  Globular clusters members with Gaia DR2 (Bustos Fierro+, 2019)
J/MNRAS/489/2079 :  Calibration of distances with deep learning (Leung+, 2019)
J/MNRAS/489/3093 :  Gaia DR2 parallax of globular clusters (Shao+, 2019)
J/MNRAS/490/157 :  Fastest stars in the Galaxy with Gaia DR2 (Marchetti+, 2019)
J/MNRAS/490/5647 :  Correlations in unbound star trajectories (Montanari+, 2019)
J/MNRAS/492/1164 :  Abundances of Gaia DR2 wide binaries (Hawkins+, 2020)
J/MNRAS/492/L40 :  Disentangling cataclysmic variables in Gaia DR2 (Abril+, 2020)
J/MNRAS/493/2339 :  Internal motions in OB-associations (Melnik+, 2020)
J/MNRAS/495/663 :  Gaia DR2 OB associations (Ward+, 2020)
J/MNRAS/496/2021 :  New Galactic open clusters (Ferreira+, 2020)
J/MNRAS/502/6080 :  New members of Cygnus OB2 from Gaia DR2 (Orellana+, 2021)
J/MNRAS/502/L90 :  Gaia DR2 Galactic bulge new star clusters (Ferreira+, 2021)
J/MNRAS/503/3660 :  Gaia Spectrophotometric Standard Stars. V. (Pancino+, 2021)
J/MNRAS/504/356 :  Updated parameters of 1743 open clusters (Dias+, 2021)
J/MNRAS/505/1135 :  Gaia/IPHAS catalogue of Ha-excess sources (Fratta+, 2021)
J/AZh/97/733 :  Binary Star Pop. with Common Proper Motion (Sapozhnikov+, 2020)
J/AZh/98/239 :  Orbits of Four Visual Binaries (Romanenko+, 2021)
J/AN/340/386 :  Catalogue of spectroscopic binary candidate stars (Jack, 2019)
J/AN/341/770 :  Corrected proper motion for HIP stars (Damljanovic, 2020)
J/AN/342/538 :  Frolov 1 and NGC 7510 CCD UBV photometry (Yontan+, 2021)
J/AcA/68/351 :  Predicted Microlensing Events by nearby VLM objects (Nielsen+ 2018)
J/other/ApSS/365.112 :  Kinematic data for high luminosity stars (Melnik+, 2020)
J/other/RAA/21.93 :  74 new open clusters found in Gaia DR2 (He+, 2021)
J/other/JDSO/15.434 :  New Binary Systems from Data of Gaia DR2 (Romanov, 2019)
WARNING: UnitsWarning: Unit 'k' not supported by the VOUnit standard. Did you mean K? [astropy.units.format.vounit]

For 2MASS we will want to use II/246 :  2MASS All-Sky Catalog of Point Sources (Cutri+ 2003) and for Gaia I/345 :  Gaia DR2 (Gaia Collaboration, 2018). Before we query the full two tables we only look at a few sources for each table to understand which columns are available. The query below will give us 50 sources each – the default for the get_catalogs method.

test_twomass = Vizier.get_catalogs("II/246")
TableList with 1 tables:
	'0:II/246/out' with 15 column(s) and 50 row(s) 
Table length=50
test_gaia = Vizier.get_catalogs("I/345")
TableList with 20 tables:
	'0:I/345/gaia2' with 32 column(s) and 50 row(s) 
	'1:I/345/rvstdcat' with 32 column(s) and 50 row(s) 
	'2:I/345/rvstdmes' with 7 column(s) and 50 row(s) 
	'3:I/345/allwise' with 2 column(s) and 50 row(s) 
	'4:I/345/iers' with 2 column(s) and 50 row(s) 
	'5:I/345/cepheid' with 25 column(s) and 50 row(s) 
	'6:I/345/rrlyrae' with 23 column(s) and 50 row(s) 
	'7:I/345/lpv' with 13 column(s) and 50 row(s) 
	'8:I/345/varres' with 9 column(s) and 50 row(s) 
	'9:I/345/shortts' with 9 column(s) and 50 row(s) 
	'10:I/345/tsstat' with 13 column(s) and 50 row(s) 
	'11:I/345/numtrans' with 4 column(s) and 50 row(s) 
	'12:I/345/transits' with 20 column(s) and 50 row(s) 
	'13:I/345/rm' with 9 column(s) and 50 row(s) 
	'14:I/345/rmseg' with 16 column(s) and 50 row(s) 
	'15:I/345/rmout' with 2 column(s) and 50 row(s) 
	'16:I/345/ssoobj' with 6 column(s) and 50 row(s) 
	'17:I/345/ssoorb' with 19 column(s) and 50 row(s) 
	'18:I/345/ssores' with 10 column(s) and 50 row(s) 
	'19:I/345/ssoobs' with 7 column(s) and 50 row(s) 
Table length=50
degmasdegmasmasmasmas / yrmas / yrmas / yrmas / yrmagmagmagmagmagmagmagkm / skm / sKmagmagsolRadsolLum

As you will see below, we only need coordinates, 2MASS photometry in the H and K band, and Gaia photometry in the Gaia G band. So we’ll query the tables II/246/out for 2MASS and I/345/gaia2 for Gaia DR2:

twomass = moc_intersection.query_vizier_table("II/246/out", max_rows=20000)
Table length=20000
gaia = moc_intersection.query_vizier_table("I/345/gaia2", max_rows=20000)
Table length=20000
degdegarcsecarcsecdegdegmasdegmasmasmasmas / yrmas / yrmas / yrmas / yre-/se-/smage-/se-/smage-/se-/smagmagkm / skm / sKmagmagRsunLsun

Step 6: Cross-match Gaia and WISE sources in all fields#

We now want to find sources in the selected region (observed in the MASH regions of interest and at low extinction) that are common to the Wide Infrared Survey Explorer (WISE) and Gaia catalogs. To do so, we will perform a cross-match of the Gaia and WISE catalogs. Alternatively, we could use the CDS XMatch service via the corresponding astroquery module.

To do so, let’s first inspect the match_coordinates_sky function from astropy.coordinates.

Help on function match_coordinates_sky in module astropy.coordinates.matching:

match_coordinates_sky(matchcoord, catalogcoord, nthneighbor=1, storekdtree='kdtree_sky')
    Finds the nearest on-sky matches of a coordinate or coordinates in
    a set of catalog coordinates.
    This finds the on-sky closest neighbor, which is only different from the
    3-dimensional match if ``distance`` is set in either ``matchcoord``
    or ``catalogcoord``.
    matchcoord : `~astropy.coordinates.BaseCoordinateFrame` or `~astropy.coordinates.SkyCoord`
        The coordinate(s) to match to the catalog.
    catalogcoord : `~astropy.coordinates.BaseCoordinateFrame` or `~astropy.coordinates.SkyCoord`
        The base catalog in which to search for matches. Typically this will
        be a coordinate object that is an array (i.e.,
        ``catalogcoord.isscalar == False``)
    nthneighbor : int, optional
        Which closest neighbor to search for.  Typically ``1`` is desired here,
        as that is correct for matching one set of coordinates to another.
        The next likely use case is ``2``, for matching a coordinate catalog
        against *itself* (``1`` is inappropriate because each point will find
        itself as the closest match).
    storekdtree : bool or str, optional
        If a string, will store the KD-Tree used for the computation
        in the ``catalogcoord`` in ``catalogcoord.cache`` with the
        provided name.  This dramatically speeds up subsequent calls with the
        same catalog. If False, the KD-Tree is discarded after use.
    idx : int array
        Indices into ``catalogcoord`` to get the matched points for each
        ``matchcoord``. Shape matches ``matchcoord``.
    sep2d : `~astropy.coordinates.Angle`
        The on-sky separation between the closest match for each
        ``matchcoord`` and the ``matchcoord``. Shape matches ``matchcoord``.
    dist3d : `~astropy.units.Quantity` ['length']
        The 3D distance between the closest match for each ``matchcoord`` and
        the ``matchcoord``. Shape matches ``matchcoord``.  If either
        ``matchcoord`` or ``catalogcoord`` don't have a distance, this is the 3D
        distance on the unit sphere, rather than a true distance.
    This function requires `SciPy <>`_ to be installed
    or it will fail.
# We generate the coordinates in the appropriate format
twomass_coord = SkyCoord(ra=twomass["RAJ2000"], dec=twomass["DEJ2000"], unit=u.deg)
gaia_coord = SkyCoord(ra=gaia["ra_epoch2000"], dec=gaia["dec_epoch2000"], unit=u.deg)

index, separation_2d, _ = match_coordinates_sky(twomass_coord, gaia_coord)
# Decide the maximum separation between objects to be considered acceptable matches
max_separation = 1.0 * u.arcsec
# Apply constraint on the two catalogs
sep_constraint = separation_2d < max_separation
twomass_matches = twomass[sep_constraint]
gaia_matches = gaia[index[sep_constraint]]
# Select only interesting columns from twomass_matches
match_catalog = twomass_matches["_2MASS", "RAJ2000", "DEJ2000", "Hmag", "Kmag"]
# Add column G magnitude from gaia
match_catalog["Gmag"] = gaia_matches["phot_g_mean_mag"]
Table length=6376

Step 7: Build a color-color diagram#

We now use the data we got from the cross-match to get a WISE/Gaia color-color diagram for all the sources in the low extinction sky regions covered by the MASH survey:

fig, ax = plt.subplots(figsize=(10, 8))
    match_catalog["Hmag"] - match_catalog["Kmag"],
    match_catalog["Gmag"] - match_catalog["Hmag"],
ax.set_xlabel("H - Ks [mag]", fontsize=16)
ax.set_ylabel("G - H [mag]", fontsize=16)