In [58]:
def shift_gfa_points(deltax, deltay):
import numpy as np
x = [-125.10482863, -129.83038525, -159.04283509, -154.31646944]
y = [-370.01790486, -384.56223777, -375.05643893, -360.51151824]
point1 = [x[2], y[2]]
point2 = [x[1], y[1]]
vector1 = [(point2[0] - point1[0]), (point2[1] - point1[1])]
vector2 = [1, 0]
# Angle between vector1 and vector 2 using dot product
angle = np.arccos((np.dot(vector1, vector2))/(np.sqrt((vector1[0]**2) + (vector1[1]**2))))
shiftmat = np.zeros(shape=(2,2))
shiftmat[0] = [np.cos(angle), -np.sin(angle)]
shiftmat[1] = [np.sin(angle), np.cos(angle)]
reverseshift= np.zeros(shape=(2,2))
reverseshift[0] = [np.cos(angle), np.sin(angle)]
reverseshift[1] = [-np.sin(angle), np.cos(angle)]
# Shifts the initial coordinates to be parallel to the vector [1, 0]
coord = np.zeros(shape=(2,1))
oldxcoord = x
oldycoord = y
for i in range(4):
coord[0] = oldxcoord[i]
coord[1] = oldycoord[i]
newcoord = np.matmul(shiftmat, coord)
oldxcoord[i] = newcoord[0]
oldycoord[i] = newcoord[1]
if(i == 0 or i == 1):
x[i] = newcoord[0] + deltax
else:
x[i] = newcoord[0] - deltax
if(i == 1 or i == 2):
y[i] = newcoord[1] - deltay
else:
y[i] = newcoord[1] + deltay
oldxcoord = x
oldycoord = y
for i in range(4):
coord[0] = oldxcoord[i]
coord[1] = oldycoord[i]
newcoord = np.matmul(reverseshift, coord)
oldxcoord[i] = newcoord[0]
oldycoord[i] = newcoord[1]
x[i] = newcoord[0]
y[i] = newcoord[1]
rotatemat = np.zeros(shape=(2,2))
rotatemat[0] = [np.cos(np.radians(36)), -np.sin(np.radians(36))]
rotatemat[1] = [np.sin(np.radians(36)), np.cos(np.radians(36))]
return find_new_gfa_coordinates(x, y, rotatemat)
def find_new_gfa_coordinates(x, y, rotatemat):
import numpy as np
x_all = np.zeros(shape=(40,1))
y_all = np.zeros(shape=(40,1))
coord = np.zeros(shape=(2,1))
gfacoord = np.zeros(shape=(4, 2))
oldxcoord = x
oldycoord = y
counter = 0
for j in range(10):
for i in range(4):
coord[0] = oldxcoord[i]
coord[1] = oldycoord[i]
newcoord = np.matmul(rotatemat, coord)
oldxcoord[i] = newcoord[0]
oldycoord[i] = newcoord[1]
gfacoord[i] = [newcoord[0], newcoord[1]]
x_all[counter] = newcoord[0]
y_all[counter] = newcoord[1]
counter += 1
return x_all, y_all
In [59]:
def degrees2xytolerance(degrees_tolerance):
# Uses the center of a given GFA from DESI-0530-v13 Excel Spreadsheet to find the tolerance
import desimodel.io
import scipy.interpolate
platescale = desimodel.io.load_platescale()
fn = scipy.interpolate.interp1d(platescale['radius'], platescale['radial_platescale'], kind = 'quadratic')
fn1 = scipy.interpolate.interp1d(platescale['radius'], platescale['az_platescale'], kind = 'quadratic')
# Center of a given GFA from DESI-0530-v13 Excel Spreadsheet
x = 333.738
y = 217.766
radius = np.sqrt(x**2 + y**2)
# Platescales are in units of microns per arcsecond
r_ps = fn(radius)
az_ps = fn(radius)
x_tolerance = degrees_tolerance / (10 **3) * 3600 * r_ps
y_tolerance = degrees_tolerance / (10**3) * 3600 * az_ps
return x_tolerance, y_tolerance
In [60]:
def retrieve_minimum_boundary(x_tolerance, y_tolerance):
import desimodel.footprint
import desimodel.focalplane
targetx = 116.279135121
targety = -372.885546514
#6.644525362152656, -9.055425745149217 GUARANTEED TO BE IN GFA (RA, DEC)
#x, y = desimodel.focalplane.radec2xy(7.11, -10.53, targetx, targety)
#print(x, y)
# If any calculated area is under the threshold area, it is mathematically impossible
THRESHOLD_AREA = 469.7
MIN_TOLERANCE = 0.001
# The area boundary's value is the area of the gfa plus some tolerance.
# x and y hold the 40 new GFA coordinates
x, y = shift_gfa_points(x_tolerance, y_tolerance)
targetx = np.asarray(targetx)
targety = np.asarray(targety)
#print(targetx)
# Method to check if point is inside the rectangle
for gfaid in range(0, 4, 4):
# a1 through a4 are edge lengths of the rectangle formed by corners of the GFAs
a1 = np.sqrt((x[gfaid] - x[gfaid + 1])**2
+ (y[gfaid] - y[gfaid + 1])**2)
a2 = np.sqrt((x[gfaid + 1] - x[gfaid + 2])**2
+ (y[gfaid + 1] - y[gfaid + 2])**2)
a3 = np.sqrt((x[gfaid + 2] - x[gfaid + 3])**2
+ (y[gfaid + 2] - y[gfaid + 3])**2)
a4 = np.sqrt((x[gfaid + 3] - x[gfaid])**2
+ (y[gfaid + 3] - y[gfaid])**2)
# b1 through b4 are the line segments from each corner to the target location
b1 = np.sqrt((x[gfaid] - targetx)**2
+ (y[gfaid] - targety)**2)
b2 = np.sqrt((x[gfaid + 1] - targetx)**2
+ (y[gfaid + 1] - targety)**2)
b3 = np.sqrt((x[gfaid + 2] - targetx)**2
+ (y[gfaid + 2] - targety)**2)
b4 = np.sqrt((x[gfaid + 3] - targetx)**2
+ (y[gfaid + 3] - targety)**2)
# Calculating areas of triangles using Heron's Formula
u1 = (a1 + b1 + b2) / 2.0
u2 = (a2 + b2 + b3) / 2.0
u3 = (a3 + b3 + b4) / 2.0
u4 = (a4 + b4 + b1) / 2.0
area1 = np.sqrt(u1 * (u1 - a1) * (u1 - b1) * (u1 - b2))
area2 = np.sqrt(u2 * (u2 - a2) * (u2 - b2) * (u2 - b3))
area3 = np.sqrt(u3 * (u3 - a3) * (u3 - b3) * (u3 - b4))
area4 = np.sqrt(u4 * (u4 - a4) * (u4 - b4) * (u4 - b1))
targetarea = area1 + area2 + area3 + area4
assert np.all(targetarea > THRESHOLD_AREA)
#print(targetarea)
return targetarea
In [61]:
# CHANGE TO WHERE IT ONLY TAKES TELRA AND TELDEC
# ALSO CHANGE TO WHERE YOURE NOT LOOPING OVER ALL TARGETS
def is_on_gfa(tileid, ra, dec, degrees_tolerance = 0.0):
"""
Checks if a target is on any of the 10 GFAs given a tileid and an array of RA and Dec pointings,
as well as a parameter for degrees of tolerance one would like to allow. When using
desimodel.footprint.find_points_in_tiles(tiles, ra, dec, radius) with this function to
check what points are on the GFAs, the default radius parameter should be set to 1.651 (degrees),
so that boundary GFA area actually encompasses points normally outside of the tile.
Parameters
----------
tileid: An integer representing the tile's unique ID which can retrieve its RA and Dec pointing
ra: An array of RA values for locations in the sky
dec: An array of declination values for locations in the sky
degrees_tolerance: A value in degrees on the sky of how much tolerance one would allow for seeing
if a target is on the gfa.
"""
import desimodel.footprint
# If any calculated area is under the threshold area, it is mathematically impossible
THRESHOLD_AREA = 469.7
MIN_TOLERANCE = 0.001
telra, teldec = desimodel.footprint.get_tile_radec(tileid)
targetx, targety = desimodel.focalplane.radec2xy(telra, teldec, ra, dec)
x_tolerance, y_tolerance = degrees2xytolerance(degrees_tolerance)
targetindices = []
gfalist = []
# x and y hold the 40 new GFA coordinates
x, y = shift_gfa_points(x_tolerance, y_tolerance)
# The area boundary's value is the area of the gfa plus some tolerance.
AREA_BOUNDARY = retrieve_minimum_boundary(x_tolerance, y_tolerance) + MIN_TOLERANCE
targetx = np.asarray(targetx)
targety = np.asarray(targety)
# Method to check if point is inside the rectangle
for gfaid in range(0, 40, 4):
# a1 through a4 are edge lengths of the rectangle formed by corners of the GFAs
a1 = np.sqrt((x[gfaid] - x[gfaid + 1])**2
+ (y[gfaid] - y[gfaid + 1])**2)
a2 = np.sqrt((x[gfaid + 1] - x[gfaid + 2])**2
+ (y[gfaid + 1] - y[gfaid + 2])**2)
a3 = np.sqrt((x[gfaid + 2] - x[gfaid + 3])**2
+ (y[gfaid + 2] - y[gfaid + 3])**2)
a4 = np.sqrt((x[gfaid + 3] - x[gfaid])**2
+ (y[gfaid + 3] - y[gfaid])**2)
# b1 through b4 are the line segments from each corner to the target location
b1 = np.sqrt((x[gfaid] - targetx)**2
+ (y[gfaid] - targety)**2)
b2 = np.sqrt((x[gfaid + 1] - targetx)**2
+ (y[gfaid + 1] - targety)**2)
b3 = np.sqrt((x[gfaid + 2] - targetx)**2
+ (y[gfaid + 2] - targety)**2)
b4 = np.sqrt((x[gfaid + 3] - targetx)**2
+ (y[gfaid + 3] - targety)**2)
# Calculating areas of triangles using Heron's Formula
#gfaarea = a1 * a2
#gfaarea2 = a2 * a3
u1 = (a1 + b1 + b2) / 2.0
u2 = (a2 + b2 + b3) / 2.0
u3 = (a3 + b3 + b4) / 2.0
u4 = (a4 + b4 + b1) / 2.0
area1 = np.sqrt((u1 * (u1 - a1) * (u1 - b1) * (u1 - b2)).clip(0))
area2 = np.sqrt((u2 * (u2 - a2) * (u2 - b2) * (u2 - b3)).clip(0))
area3 = np.sqrt((u3 * (u3 - a3) * (u3 - b3) * (u3 - b4)).clip(0))
area4 = np.sqrt((u4 * (u4 - a4) * (u4 - b4) * (u4 - b1)).clip(0))
targetarea = area1 + area2 + area3 + area4
"""TEMPORARILY REMOVING THE ASSERT BECAUSE OF THE RUNTIME WARNING"""
#assert np.all(targetarea > THRESHOLD_AREA)
#error = np.where(np.isnan(targetarea))
#print(targetarea[error])
#print(area1[error])
#print(area2[error])
#print(area3[error])
#print(u3[error])
#print(a3)
#print(b3[error])
#print(b4[error])
#print(u3[error] - b3[error])
#print(area4[error])
#print(np.where(targetarea < THRESHOLD_AREA))
if(any(targetarea < AREA_BOUNDARY) and all(targetarea > THRESHOLD_AREA)):
# Does not need to print what target the GFA is on, but should consider
#returning this information, along with its tileid
#print("The target is on GFA", int(gfaid / 4))
gfalist.extend([gfaid / 4])
newtargetindices = np.where(targetarea < AREA_BOUNDARY)
targetindices.append(newtargetindices[0])
"""
Returns:
targetindices is a list of arrays of targets with their respective indices in the
RA and Dec list passed in that fall on certain GFAs denoted by the index
in the targetindices list.
gfaid is a list of values from 0-9 representing the gfa that a target is on.
"""
#print(gfatable['X'], gfatable['Y'])
return targetindices, gfalist
In [62]:
import matplotlib.pyplot as plt
import numpy as np
#telra and teldec = 7.11 -10.53
mindec = -12.5
maxdec = -7
minra = 5.4
maxra = 8.9
ra = []
dec = []
import desimodel.focalplane
import desimodel.footprint
while mindec < maxdec:
startra = minra
while startra < maxra:
ra.append(startra)
dec.append(mindec)
startra += .02
mindec += .02
import desimodel.io
tiles = desimodel.io.load_tiles()
ra = np.asarray(ra)
dec = np.asarray(dec)
"""index = 0
for specific_tile in tiles:
if(tiles[index]['TILEID'] == 23658):
print(index)
index += 1
#INDEX IS TILES[8254]"""
#points = desimodel.footprint.find_points_in_tiles(tiles[8254], ra, dec, 1.651)
"""for lists in points:
if lists:
print(lists)"""
#print(points)
new_ra = []
new_dec = []
"""for point_indices in points:
#plt.plot(ra[point_indices], dec[point_indices], 'b.')
new_ra.append(ra[point_indices])
new_dec.append(dec[point_indices])"""
"""targetindices = is_on_gfa(23658, new_ra, new_dec, .1)
for gfaid in targetindices:
for indices in gfaid:
plt.plot(new_ra[indices], new_dec[indices], 'r.')
targetindices = is_on_gfa(23658, new_ra, new_dec, 0)
for gfaid in targetindices:
for indices in gfaid:
plt.plot(new_ra[indices], new_dec[indices], 'g.')"""
Out[62]:
In [63]:
# invalid value encountered in sqrt
# invalid value encountered in less
# invalid value encountered in greater
# Could this be a NaN problem?
# Identify the list of GFAs with the 100 brightest objects on them
def findBrightestObjects(sweepname):
from astropy.io import fits
from astropy.table import Table
sweep = Table.read(sweepname)
import desimodel.io
import desimodel.footprint
tiles = desimodel.io.load_tiles()
points = desimodel.footprint.find_points_in_tiles(tiles, sweep['RA'], sweep['DEC'])
counter = 0
for lists in points:
if lists:
tileid = tiles[counter]['TILEID']
#print(tileid)
#print(lists)
testra = []
testdec = []
telra, teldec = desimodel.footprint.get_tile_radec(tileid)
#print(telra, teldec)
for indices in lists:
testra.append(sweep[indices]['RA'])
testdec.append(sweep[indices]['DEC'])
#print(type(sweep[indices]['RA']))
if(type(sweep[indices]['RA']) is float):
print(sweep[indices]['RA'])
targetindices, gfaid = is_on_gfa(tileid, testra, testdec, 0)
#print(testra)
#print(testdec)
#print(targetindices)
counter += 1
sweepname = 'sweep-280p050-290p055.fits'
findBrightestObjects(sweepname)
#tiles
In [ ]:
# Identify the list of GFAs with the 100 brightest objects on them
def brightest_objects(tiles, ra, dec, flux_r):
import desimodel.footprint
points = desimodel.footprint.find_points_in_tiles(tiles, ra, dec)
alltargetindices = []
brightestindices = []
counter = 0
for lists in points:
if lists:
tileid = tiles[counter]['TILEID']
testra = []
testdec = []
#telra, teldec = desimodel.footprint.get_tile_radec(tileid)
for indices in lists:
testra.append(ra[indices])
testdec.append(dec[indices])
"""if(type(sweep[indices]['RA']) is not float):
print(sweep[indices]['RA'])"""
targetindices, gfaid = is_on_gfa(tileid, testra, testdec, 0)
#blankarray = np.asarray([])
#entirelist = sum(targetindices, blankarray)
alltargetindices.append(targetindices)
counter += 1
#print(alltargetindices)
for i in range(100):
fluxindex = 0
largestfluxr = 0
for listofarrays in alltargetindices:
for gfaarrays in listofarrays:
for celestialobjects in gfaarrays:
#print(flux_r[celestialobjects])
if(celestialobjects not in brightestindices and flux_r[celestialobjects] > largestfluxr):
largestfluxr = flux_r[celestialobjects]
fluxindex = celestialobjects
brightestindices.append(fluxindex)
print(brightestindices)
return brightestindices
from astropy.io import fits
from astropy.table import Table
sweepname = 'sweep-280p050-290p055.fits'
sweep = Table.read(sweepname)
import desimodel.io
tiles = desimodel.io.load_tiles()
indices = brightest_objects(tiles, sweep['RA'], sweep['DEC'], sweep['FLUX_R'])
In [ ]:
def find_valid_gfas(tiles, ra, dec, flux_r):
import numpy as np
import desimodel.footprint
points = desimodel.footprint.find_points_in_tiles(tiles, ra, dec)
allvalidgfas = []
allinvalidgfas = []
# Using equation 22.5 - 2.5*log10(flux_r) < 15
MAX_BRIGHTNESS = 15
indexcounter = 0
for lists in points:
validgfas = []
invalidgfas = []
tileid = tiles[indexcounter]['TILEID']
testra = []
testdec = []
for indices in lists:
testra.append(ra[indices])
testdec.append(dec[indices])
targetindices, gfaid = is_on_gfa(tileid, testra, testdec)
gfacounter = 0
for gfaarrays in targetindices:
for celestialobjects in gfaarrays:
if((22.5 - 2.5*np.log10(flux_r[celestialobjects])) > MAX_BRIGHTNESS):
#print("This target is too bright and might cause issues for observing!")
invalidgfas.extend([gfaid[gfacounter]])
continue
if(gfaid[gfacounter] not in invalidgfas):
validgfas.extend([gfaid[gfacounter]])
gfacounter += 1
for gfas in range(10):
if gfas not in gfaid:
invalidgfas.extend([gfas])
indexcounter += 1
allvalidgfas.append(validgfas)
allinvalidgfas.append(invalidgfas)
return allvalidgfas, allinvalidgfas
from astropy.io import fits
from astropy.table import Table
sweepname = 'sweep-290p065-300p070.fits'
sweep = Table.read(sweepname)
import desimodel.io
tiles = desimodel.io.load_tiles()
validgfas, invalidgfas = find_valid_gfas(tiles, sweep['RA'], sweep['DEC'], sweep['FLUX_R'])
print(validgfas)
In [ ]:
%pylab
from astropy.io import fits
from astropy.table import Table
import desimodel.io
tiles = desimodel.io.load_tiles()
sweepname = 'sweep-280p050-290p055.fits'
#sweepname = 'sweep-290p065-300p070.fits'
sweep = Table.read(sweepname)
len(sweep['RA'])
indesi = desimodel.footprint.is_point_in_desi(tiles, sweep['RA'], sweep['DEC'])
print(indesi)
print(np.count_nonzero(indesi))
plot(sweep['RA'], sweep['DEC'], 'b.')
plot(sweep['RA'][indesi], sweep['DEC'][indesi], 'g.')