Working with Custom Point Clouds as Index

[Leonard Freißmuth]

I'm trying to get astrometry to work with a custom index, which is comparably small with respect to the usual universe indices: only 500 points.

Granted, this is probably a setting astrometry has not been designed for, but it would be interesting to see if the algorithm translates to this setting, too.

To reproduce, I'm generating a custom fits file of a random map like the following:

import numpy as np

import fitsio

import os

import pathlib

import astrometry

# generate a random map spanning the entire sky

num_stars = 500

full_map = np.random.rand(num_stars, 2) * np.array([180, 360]) - np.array([90, 180])

# convert map into data to be read by fitsio

data = np.zeros(len(full_map), dtype=[

('id', 'i8'), ('ra', 'f8'), ('dec', 'f8'), ('mag', 'f4')

])

for i in range(num_stars):

data[i] = (i, full_map[i, 0], full_map[i, 1], 1.0)

fitsio.write('test.fits', data)

# generate quads for .fits file

os.system(f"build-astrometry-index -i test.fits -P 14 -o test.fits")

solver = astrometry.Solver([pathlib.Path('test.fits')])

I'm seeking to lookup maps in the realm of 64 arcmins, so I'm using -P 14.

Also, I want all stars to be used for matching, so I'm giving equal mag/weights.

Now I want to find a local submap (stars within circle of radius 32 around origin) of the bigger map, which has no rotation, scale, translation or noise applied:

# generate and lookup local map

local_map = full_map[np.linalg.norm(full_map, axis=1) < 32]

print(len(local_map))

solution = solver.solve(

local_map,

size_hint=None,

position_hint=None,

solution_parameters=astrometry.SolutionParameters()

)

print(solution.has_match())

Yet, I can't make astrometry find a solution, the print statement always evaluates False.

Is there something wrong with my configuration? Or is astrometry just not designed to work with maps this sparse? In this case, do you have an idea, how I could make it work, maybe adjust some parameters, or let the full map span only a small area of the sky? (Didn't work thus far)

Any help is highly appreciated!

[Dustin Lang]

Hi,

First up: you've got RA going from -90 to +90 and Dec from -180 to +180... that's backwards.

I feel like your math on the preset size doesn't work...  the sky is about 40,000 square degrees = 144 million square arcmin.  If you divide that into 64 x 64 arcmin patches, there are 35,000 patches like that.  So 500 stars is not nearly enough... you want about 10 stars per patch.

More broadly:  Astrometry.net assumes a *pinhole camera projection*...   the way you are pulling your "local_map" stars is a very different projection (plate caree or something), so that's definitely not going to work.

If you shrink your whole 2-d plane down to a much smaller region where the projection effects are smaller, that is likely to work better.

cheers,

dustin

--
You received this message because you are subscribed to the Google Groups "astrometry" group.
To unsubscribe from this group and stop receiving emails from it, send an email to astrometry+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/astrometry/c17c9574-df95-4c22-a259-4a341cc38a8dn%40googlegroups.com.

[Leonard Freißmuth]

Hi Dustin, thanks for the quick reply!

you're totally right about the RA, I fixed that. Also, I mistook the units, I meant 64 x 64 deg patches, not arcmin, sorry. This leaves around 10 patches So my number of stars should theoretically suffice then.

Not sure why the projection matters here, I just assume a random distribution of the stars across the variable range of RA and DEC. Is something wrong with this? I only use the norm to select a circular sub-section of points, the coordinates of the local_map's stars are exactly the same as in the index.

Considering your comment about the projection effects---thanks for that---I reconfigured the experiment so the random map only spans +- 0.5 deg and the local_map is a circular section of this with radius 0.2 deg (see attached picture)

As such, the query stars span 24 arcmins so I now chose -P 4 as per the documentation. And indeed, build-astrometry-index now builds almost 400 quads in total, so this seems to work nicely.

However, the solver still can't manage to find a solution (see the attached code for reproducing).

Am I missing something about the solver? In my mind, everything should work now: I have loads of quads and a query set that spans a substantial amount of the full index. Is there something else about my setup that I need to take care of?

Thanks again for the excellent support!

import numpy as np

import fitsio

import os

import pathlib

import astrometry

# generate a random map spanning the entire sky

num_stars = 500

full_map = np.random.rand(num_stars, 2) * np.array([1, 1]) - np.array([0.5, 0.5])

# convert map into data to be written by fitsio

data = np.zeros(len(full_map), dtype=[

('id', 'i8'), ('ra', 'f8'), ('dec', 'f8'), ('mag', 'f4')

])

for i in range(num_stars):

data[i] = (i, full_map[i, 1], full_map[i, 0], 1.0)

fitsio.write('test.fits', data)

# generate quads for .fits file

os.system(f"build-astrometry-index -i test.fits -P 4 -o test.fits")

solver = astrometry.Solver([pathlib.Path('test.fits')])

# generate and lookup local map

local_map = full_map[np.linalg.norm(full_map, axis=1) < 0.2]

solution = solver.solve(

local_map,

size_hint=None,

position_hint=None,

solution_parameters=astrometry.SolutionParameters()

)

print("Has solution:", solution.has_match())

[Dustin Lang]

Hi,

Remember that astrometry.net is meant for matching images to the reference catalog ... so the "local_map" coordinates are in pixel coordinates.  There are some parameters whose defaults only make sense for units of pixels -- for example, there is an assumed uncertainty on the position of each star.  +-1 pixel is reasonable, but your local_map coordinates all have |values| less than 1, so +-1 is not reasonable any more :)

Maybe try shifting and scaling your local_map values so that they fit in a 1000 x 1000 box?

One other thing ... you have made everything the same brightness.  This is going to make things difficult.  Astrometry.net searches through the stars in your local_map in brightness order.  By default, the index files use 4-dimensional features - quadrangles.  So there are, very roughly, N^4 quadrangles to choose from given N stars in your local_map.  That's a big number.

If you give it hints about the scale, then that N^4 drops down a lot.

You could build the index using triangles (build-astrometry-index -d 3).

Also, for bulid-astrometry-index on just a small region of sky, the "-E" flag should make it run faster.

cheers,

dustin

To view this discussion on the web visit https://groups.google.com/d/msgid/astrometry/51d416c1-1094-48d2-82c4-f3557935e892n%40googlegroups.com.

[Leonard Freißmuth]

Hi Dustin,

thanks again for your answer. The scaling to pixel space makes a lot of sense, and solved the problem partially but not satisfyingly: Now I can sometimes match sub_maps but only if they are covering a large portion of the index (empirically more than 30 %).

Today, I have dug a bit deeper and plotted the quads using the produced fits file. I think I now know the core problem: Somehow, I can't make the sampling of the quads remain local.

I played around with setting the flags and reading the fits file: 

• I again used a random index of 500 stars spanning [-0.5, 0.5] x [-0.5, 0.5] as described in previous posts and the following command string:
  build-astrometry-index -i test.fits -N 1823 -l 0 -u 6 -p 128 -R 32 -L 64 -n 10 -E -o test.fits
• I tried to get a reasonable size for Healpixes using the flag -N, but am not sure if what I'm targeting is correct. If I understood the docs correctly, I want 10 stars per healpix, so with 500 stars per deg² this would yield a total number of 1823 healpixes (41253 deg² / (10 stars / 500 stars/deg²)). That's why I chose -N 1823.
  
• Crucially, I want to force the quads to remain local and not span the entire index. For this, I want to set an upper limit for quad size of say a tenth of the entire index, which would be 6 arcmins. Not wanting to restrict how small quads are, I choose -l 0 -u 6 as flags.
• I read the generated fits file to read the indices of the quads. I assumed row major ordering, so I converted the quads to a numpy array using
  
  quads = np.frombuffer(quads_fits_file[1].data, dtype=np.int32).reshape(-1, 4)
  I looked up the indices in the original map.
  

With all these settings (and many variations through trial and error) I always achieved quads looking like this:

I also tried different sizes of the index also spanning e.g.10x10 deg²  but the quads always span the entire map, which might explain the still very poor lookup performance. 

My questions now are: 

1. Is this expected behaviour?
2. Do the Parameters I chose make sense?
3. Do you know of ways how I can achieve an ideal quad coverage that looks more like the one you show in your docs?
   

Many, many thanks again for your excellent support on this.

Cheers,

Leonard

P.S. putting the code for reproducing here:

import os

import numpy as np

from matplotlib import pyplot as plt

import fitsio

import astropy

import astrometry

# generate a random map spanning the entire sky

num_stars = 500

full_map = np.random.rand(num_stars, 2) - 0.5

# convert map into data to be written by fitsio

data = np.zeros(len(full_map), dtype=[

('id', 'i8'), ('ra', 'f8'), ('dec', 'f8'), ('mag', 'f4')

])

for i in range(num_stars):

data[i] = (i, full_map[i, 1], full_map[i, 0], 1.0)

fitsio.write('test.fits', data)

# generate quads for .fits file

os.system(

f"build-astrometry-index -i test.fits -N 1823 -l 0 -u 6 -p 128 -R 32 -L 64 -n 10 -T -E -o test.fits"

)

# read quads (field for quads is at index 1 of fits_file)

fits_file = astropy.io.fits.open('test.fits')

quads = np.frombuffer(fits_file[1].data, dtype=np.int32).reshape(-1, 4)

# plot index and quads

plt.scatter(full_map[:, 0], full_map[:, 1], s=1)

for quad in quads[np.random.permutation(len(quads))[:20]]:

plt.plot(

[

full_map[quad[0], 0],

full_map[quad[1], 0],

full_map[quad[2], 0],

full_map[quad[3], 0],

full_map[quad[0], 0]

],

[

full_map[quad[0], 1],

full_map[quad[1], 1],

full_map[quad[2], 1],

full_map[quad[3], 1],

full_map[quad[0], 1]

],

linewidth=0.5

)

plt.xlabel("RA (deg)")

plt.ylabel("DEC (deg)")

plt.show()

[Dustin Lang]

Hi,

That healpix number isn't right.  For healpix number N, there are 12 * N*2 healpixel on the sky.  According to the "build-astrometry-index -h" help text, the preset for scale 0 (-P 0) uses -N 1760.

I haven't looked carefully at your quad-plotting code, but note that the order of the stars gets permuted when building the index, so the indices in the quad file do NOT refer to the index in your original reference catalog.

(You could verify whether it's working correctly by limiting the allowed quad scale to a small range and checking whether the quads you read out of the file satisfy this criterion.)

cheers,

dustin

To view this discussion on the web visit https://groups.google.com/d/msgid/astrometry/40895113-1f84-407c-b99a-d3e707e8a6d1n%40googlegroups.com.