desihub / desimeter

DESI coordinates and transformations
BSD 3-Clause "New" or "Revised" License
2 stars 4 forks source link

Bug in desimeter when running with desiconda/20200801-1.4.0-spec #114

Closed julienguy closed 4 years ago

julienguy commented 4 years ago

I get only the right half of the focal plane matched when running

desi_fvc_proc --zbfit -i fvc-00052009.fits.fz --extname F0002 --field-model fm-00052009-20200730.json --expected coordinates-00052009.fits -o fvc-00052009-20200817.csv

with the new desiconda install, loaded with,

module use /global/common/software/desi/${NERSC_HOST}/desiconda/20200801-1.4.0-spec/modulefiles
module load desiconda

It works fine with the old desiconda install

module use /global/common/software/desi/${NERSC_HOST}/desiconda/20190804-1.3.0-spec/modulefiles
module load desiconda

So it's either a different interpretation of the code between different python versions, or some bug with an I/O module.

julienguy commented 4 years ago

I have isolated the bug with a script that uses only astropy.table and scipy.spatial.cKDTree. It occurs when there are NaN in the astropy table coordinates.

failing:

source /global/cfs/cdirs/desi/software/desi_environment.sh 20.7
cd /global/cfs/cdirs/desi/users/jguy/match/
./testmatch.py

ok:
source /global/cfs/cdirs/desi/software/desi_environment.sh 19.12
cd /global/cfs/cdirs/desi/users/jguy/match/
./testmatch.py

Content of testmatch.py:

#!/usr/bin/env python

import numpy as np
import  matplotlib.pyplot as plt
from scipy.spatial import cKDTree as KDTree

# not working with desiconda 20200801-1.4.0 but working with desiconda 20190804-1.3.0-spec ???
# with desiconda 20200801-1.4.0 only x<0 is matched (!)
from astropy.table import Table
t1=Table.read("xy1.csv")
t2=Table.read("xy2.csv")
#t2=Table.read("xy2-no-nan.csv") # working fine when no NaN in table
x1=t1["x"]
y1=t1["y"]
x2=t2["x"]
y2=t2["y"]

xy1=np.array([x1,y1]).T
xy2=np.array([x2,y2]).T
tree2 = KDTree(xy2)
distances,indices_2 = tree2.query(xy1,k=1)
is_matched = (distances<10.)&(indices_2>=0)

plt.figure()
plt.plot(x1,y1,"o",alpha=0.5)
plt.plot(x2,y2,".",alpha=0.5,color="k")
plt.plot(x1[is_matched],y1[is_matched],"o",alpha=0.5)
plt.plot(x2[indices_2[is_matched]],y2[indices_2[is_matched]],"x",alpha=0.5,c="r")
plt.show()
julienguy commented 4 years ago

fixed in PR #115