/
desi-sky-locations.py
192 lines (164 loc) · 6.79 KB
/
desi-sky-locations.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
import os
import sys
import functools
from collections import Counter
import pylab as plt
import numpy as np
import fitsio
from astrometry.util.fits import fits_table, merge_tables
from astrometry.util.util import Tan
from astrometry.util.starutil_numpy import degrees_between
from astrometry.util.plotutils import plothist
from astrometry.libkd.spherematch import match_radec, tree_build_radec, tree_search_radec, tree_open
from astrometry.util.resample import resample_with_wcs, NoOverlapError
from astrometry.util.multiproc import multiproc
sys.path.insert(0, 'legacypipe/py')
from legacypipe.gaiacat import GaiaCatalog
from legacypipe.reference import fix_tycho, fix_gaia, merge_gaia_tycho
from legacypipe.survey import get_git_version
sys.path.insert(0, 'desiutil/py')
from desiutil.brick import Bricks
def run_one(X):
k, sb, bricks, version = X
print(k, sb.brickname)
outfn = 'skybricks/sky-%s.fits.fz' % sb.brickname
if os.path.exists(outfn):
return True
I = np.flatnonzero((bricks.ra2 > sb.ra1) * (bricks.ra1 < sb.ra2) * (bricks.dec2 > sb.dec1) * (bricks.dec1 < sb.dec2))
if len(I) == 0:
print('No bricks overlap')
return False
# 3600 + 1% margin on each side
w,h = 3672,3672
if sb.dec >= 78.:
# In order to have fully-overlapping tiles at high Decs, we
# need larger maps. DR9 reaches a max skytile of Dec=+85,
# where the required size is 3724. Add a little margin.
w,h = 3744,3744
binning = 4
# pixscale
cd = 1./3600.
fullw,fullh = w*binning, h*binning
fullcd = cd/binning
# There are really three states: no coverage, blob, no blob.
# Since blobs outside each brick's unique area do not appear in
# the files, we start skyblobs as zero, but also track the
# coverage so we can set !coverage to blob at the end.
skyblobs = np.zeros((fullh, fullw), bool)
covered = np.zeros((fullh, fullw), bool)
skywcs = Tan(sb.ra, sb.dec, (fullw+1)/2., (fullh+1)/2., -fullcd, 0., 0., fullcd, float(fullw), float(fullh))
for i in I:
brick = bricks[i]
#print('Blob', brickname)
#fn = 'cosmo/data/legacysurvey/dr9/%s/metrics/%s/blobs-%s.fits.gz' % (brick.hemi, brick.brickname[:3], brick.brickname)
#fn = 'dr9-south-blobs/blobs-%s.fits.gz' % (brick.brickname)
fn = '/global/cfs/cdirs/cosmo/data/legacysurvey/dr9/%s/metrics/%s/blobs-%s.fits.gz' % (brick.hemi, brick.brickname[:3], brick.brickname)
blobs,hdr = fitsio.read(fn, header=True)
wcs = Tan(hdr)
blobs = (blobs > -1)
try:
Yo,Xo,Yi,Xi,_ = resample_with_wcs(skywcs, wcs)
except NoOverlapError:
continue
skyblobs[Yo,Xo] |= blobs[Yi,Xi]
# coverage: nexp > 0 in any band
for band in ['g','r','z']:
fn = ('/global/cfs/cdirs/cosmo/data/legacysurvey/dr9/%s/coadd/%s/%s/legacysurvey-%s-nexp-%s.fits.fz' %
(brick.hemi, brick.brickname[:3], brick.brickname, brick.brickname, band))
if not os.path.exists(fn):
continue
nexp = fitsio.read(fn)
covered[Yo,Xo] |= (nexp[Yi,Xi] > 0)
if not np.any(covered):
print('No coverage')
return False
# No coverage = equivalent to there being a blob there (ie,
# conservative for placing sky fibers)
skyblobs[covered == False] = True
# bin down, counting how many times 'skyblobs' is set
subcount = np.zeros((h,w), np.uint8)
for i in range(binning):
for j in range(binning):
subcount += skyblobs[i::binning, j::binning]
subwcs = Tan(sb.ra, sb.dec, (w+1)/2., (h+1)/2., -cd, 0., 0., cd, float(w), float(h))
hdr = fitsio.FITSHDR()
hdr.add_record(dict(name='SB_VER', value=version, comment='desi-sky-locations git version'))
subwcs.add_to_header(hdr)
fitsio.write(outfn, subcount, header=hdr, clobber=True,
compress='GZIP', tile_dims=(256,256))
print('Wrote', outfn)
return True
# if not os.path.exists(sbfn):
# skybricks = Bricks(bricksize=1.0)
# skybricks.to_table().write(sbfn)
def main():
sbfn = 'skybricks.fits'
SB = fits_table(sbfn)
Bnorth = fits_table('/global/cfs/cdirs/cosmo/data/legacysurvey/dr9/north/survey-bricks-dr9-north.fits.gz')
Bsouth = fits_table('/global/cfs/cdirs/cosmo/data/legacysurvey/dr9/south/survey-bricks-dr9-south.fits.gz')
Bnorth.cut(Bnorth.survey_primary)
Bsouth.cut(Bsouth.survey_primary)
#Bsouth.cut(Bsouth.dec > -30)
Bnorth.hemi = np.array(['north']*len(Bnorth))
Bsouth.hemi = np.array(['south']*len(Bsouth))
B = merge_tables([Bnorth, Bsouth])
# Rough cut the skybricks to those near bricks.
I,J,d = match_radec(SB.ra, SB.dec, B.ra, B.dec, 1., nearest=True)
SB.cut(I)
import argparse
parser = argparse.ArgumentParser()
#parser.add_argument('--brick', help='Sky brick name')
parser.add_argument('--minra', type=float, help='Cut to a minimum RA range of sky bricks')
parser.add_argument('--maxra', type=float, help='Cut to a maximum RA range of sky bricks')
parser.add_argument('--threads', type=int, help='Parallelize on this many cores')
opt = parser.parse_args()
if opt.minra:
SB.cut(SB.ra >= opt.minra)
if opt.maxra:
SB.cut(SB.ra <= opt.maxra)
version = get_git_version(os.getcwd())
print('Version string:', version)
# Find bricks near skybricks, as a rough cut.
# (magic 1. degree > hypot(skybrick radius, brick radiu) ~= 0.9)
Inear = match_radec(SB.ra, SB.dec, B.ra, B.dec, 1., indexlist=True)
args = []
k = 0
Isb = []
for isb,(sb,inear) in enumerate(zip(SB, Inear)):
if inear is None:
continue
args.append((k, sb, B[np.array(inear)], version))
k += 1
Isb.append(isb)
print(len(args), 'sky bricks')
SB.cut(np.array(Isb))
if opt.threads:
mp = multiproc(opt.threads)
else:
mp = multiproc()
exist = mp.map(run_one, args)
if opt.minra is None and opt.maxra is None:
exist = np.array(exist)
SB[exist].writeto('skybricks-exist.fits')
return
if __name__ == '__main__':
#main()
#B = fits_table('bricks-near.fits')
#SB = fits_table('skybricks.fits')
B = fits_table('/global/cfs/cdirs/cosmo/data/legacysurvey/dr9/north/survey-bricks-dr9-north.fits.gz')
B.cut(B.dec > 77)
B.hemi = np.array(['north']*len(B))
SB = fits_table('/global/cfs/cdirs/desi/target/skybricks/v3/skybricks-exist.fits')
#I = np.flatnonzero(np.isin(SB.brickname, [
# '2187p340', '2175p340',
#'1500p290', '1507p300', '1519p300',
#'1505p310', '1490p320'
#]))
I = np.flatnonzero(SB.dec >= 78.)
version = get_git_version(os.getcwd())
#version = 4
mp = multiproc(32)
args = []
for i in I:
args.append((0, SB[i], B, version))
mp.map(run_one, args)