为什么我的代码在每个 运行 之后从相同图像中的相同列表中找到更多对象?
Why is my code finding more objects from the same list within the same images after each run?
我写了一个代码,它查看一个预先存在的图像文件夹,并使用一个恒定的对象名称列表和它们在天空中相应的 ra、dec 位置来定位它们在每个原始图像中并制作一个 10x10 arcsec剪切(如果对象在图像中)。它正在工作,我得到了一些漂亮的剪纸,但出于某种奇怪的原因,每次我 运行 它都会保存新的剪纸图像!我真的不知道为什么会这样,因为对象列表及其 ra、dec 位置始终完全相同,而且我总是使用图像和对象的确切名称保存剪切图像,这些名称不应更改。所有原始图像也保持不变。
我进行了很多测试,但仍然感到困惑——通过我的测试,我确认对象列表 (objs
) 在每个 运行 中保持相同,并且我达到了同样的结果原始图像列表 (all_images
) 和 ra、dec 位置(列表 ras_hms
和 decs_deg
)的结论。
原始图片和对象的数量都比较长,所以我运行我的代码在较小的子集上进行测试,并且在每个运行期间出现新的剪切图像的问题仍然存在.我 运行 下面这些图像的代码:'calexp-HSC-I-18012-1,7.fits', 'calexp-HSC-I-18114-0,0.fits', 'calexp-HSC-I-18114-1,1.fits'
保存在目录 /Users/myuser/Desktop/Original_Images/
中。我 运行 将我的代码放在另一个目录中,最后还保存了剪切图。第一次 运行,我生成了这些切口:'cutout-IMG-HSC-I-18114-0,0-OBJ-NEP175719.1+652743.2.fits',
'cutout-IMG-HSC-I-18114-1,1-OBJ-NEP175509.2+653523.9.fits'
。几分钟后,当我 运行 完全相同的代码而不做任何更改时,我还生成了这两个新图像:'cutout-IMG-HSC-I-18114-0,0-OBJ-NEP175654.7+652930.2.fits',
'cutout-IMG-HSC-I-18114-1,1-OBJ-NEP175458.4+653419.1.fits'
,类似地用于未来的 运行s.
如您所见,它一定没有在我的一张图像中找到任何对象(这很好),但是每个 运行 它都以某种方式在其他每个图像中找到了一个新对象(确实,每次我 运行 这个小子集上的这个代码,我都会看到两个新的剪切图像被保存,具有不同的对象名称)。我很难过,因为就像我说的那样,它搜索的对象和坐标在每个图像中都是相同的 运行。任何见解或猜测将不胜感激!
import astropy
from astropy.nddata.utils import Cutout2D, NoOverlapError
import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u
from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord, Angle
import re
import glob
def make_cutouts(img_file, box_len):
# Image data
hdulist = fits.open(img_file)
img_data = fits.getdata(img_file)
img_name = re.search(r'calexp\-(.+)\.fits', img_file)[1]
# Make cutouts (by locating coords in image with WCS)
wcs = WCS(hdulist[1].header)
for i in range(len(objs)):
# Skip if cutout already exists
if 'cutout-IMG-'+img_name+'-OBJ-'+objs[i]+'.fits' in glob.glob('cutout*.fits'):
print('Cutout of object already exists for this image, skipping...')
continue
# Convert ra, dec to HMS for specific object
ra_h = re.search(r'h=(\d+.?\d*)', str(ras_hms[i]))[1]
ra_m = re.search(r'm=(\d+.?\d*)', str(ras_hms[i]))[1]
ra_s = re.search(r's=(\d+.?\d*)', str(ras_hms[i]))[1]
ra_angle = Angle((float(ra_h), float(ra_m), float(ra_s)), unit='hourangle')
dec_angle = decs_deg[i]
# Coordinate transformation to pixels
center = SkyCoord(ra_angle, dec_angle, frame='fk5')
xp, yp = astropy.wcs.utils.skycoord_to_pixel(center, wcs=wcs, origin=1)
# Make cutout, skip when object is not in image
size = u.Quantity((box_len,box_len),u.arcsec)
try:
co = Cutout2D(img_data,(xp, yp),size,wcs=wcs)
except NoOverlapError:
continue
hdu = fits.PrimaryHDU(data=co.data,header=co.wcs.to_header())
hdu.writeto('cutout-IMG-'+img_name+'-OBJ-'+objs[i]+'.fits', overwrite=True)
return
# Gather all original images
all_images = glob.glob('/Users/myuser/Desktop/Original_Images/calexp*.fits')
coords_file = 'good_dataset.fits'
# Coordinates
hdul = fits.open(coords_file)
coords_data = hdul[1].data
objs = coords_data['Name']
ras = np.array(coords_data['RA']) # in decimal degrees
decs = np.array(coords_data['DEC']) # in decimal degrees
# Convert coordinate systems using astropy
decs_deg = Angle(decs, unit=u.deg)
ras_deg = Angle(ras, unit=u.deg)
ras_hms = [ra.hms for ra in ras_deg]
count=0
for image in all_images:
make_cutouts(image, 10.0)
count+=1
print('Image %d out of %d completed' % (count, len(all_images)))
这里是我刚才的运行打印语句的输出示例,它再次生成了两个新的剪切图像(不同的对象,相同的两个图像)...这里是图像 2没有发现任何物体的地方。此外,有趣的是,每个图像的 "already exists, skipping" 语句的数量增加了两个 运行.
Cutout of object already exists for this image, skipping...
Cutout of object already exists for this image, skipping...
Cutout of object already exists for this image, skipping...
Cutout of object already exists for this image, skipping...
Image 1 out of 3 completed
Image 2 out of 3 completed
Cutout of object already exists for this image, skipping...
Cutout of object already exists for this image, skipping...
Cutout of object already exists for this image, skipping...
Cutout of object already exists for this image, skipping...
Image 3 out of 3 completed
这是一个简单的错误:在 for
循环的末尾有一个 return
语句,这意味着 make_cutouts
的每个 运行 被限制为最多产生一个切口。每次你 运行 它,它都会产生第一个切口,然后下一次,它看到一个存在,用 continue
语句跳过它,然后得到下一个,然后退出。删除 return
语句,代码可能会正常工作。
但是,有几件事您应该尽量避免:
(1) 您在函数中使用了全局变量。您最好将 objs
和 ras_hms
作为参数传递给函数,而不是依赖全局状态来访问它们。
(2) 当您可以遍历对象本身时,您正在遍历索引,即 for thisobj, thisra in zip(objs, ras_hms):
(3) 更小,但 if 'cutout-IMG-'+img_name+'-OBJ-'+objs[i]+'.fits' in glob.glob('cutout*.fits'):
比 if os.path.exists(if 'cutout-IMG-'+img_name+'-OBJ-'+objs[i]+'.fits' in glob.glob('cutout*.fits'):
更有效。如果您使用 'cutout-IMG-{img_name}-OBJ-{obj_id}.fits'.format(img_name=img_name, obj_id=objs[i])
作为字符串
,您可能还会发现它更具可读性
我写了一个代码,它查看一个预先存在的图像文件夹,并使用一个恒定的对象名称列表和它们在天空中相应的 ra、dec 位置来定位它们在每个原始图像中并制作一个 10x10 arcsec剪切(如果对象在图像中)。它正在工作,我得到了一些漂亮的剪纸,但出于某种奇怪的原因,每次我 运行 它都会保存新的剪纸图像!我真的不知道为什么会这样,因为对象列表及其 ra、dec 位置始终完全相同,而且我总是使用图像和对象的确切名称保存剪切图像,这些名称不应更改。所有原始图像也保持不变。
我进行了很多测试,但仍然感到困惑——通过我的测试,我确认对象列表 (objs
) 在每个 运行 中保持相同,并且我达到了同样的结果原始图像列表 (all_images
) 和 ra、dec 位置(列表 ras_hms
和 decs_deg
)的结论。
原始图片和对象的数量都比较长,所以我运行我的代码在较小的子集上进行测试,并且在每个运行期间出现新的剪切图像的问题仍然存在.我 运行 下面这些图像的代码:'calexp-HSC-I-18012-1,7.fits', 'calexp-HSC-I-18114-0,0.fits', 'calexp-HSC-I-18114-1,1.fits'
保存在目录 /Users/myuser/Desktop/Original_Images/
中。我 运行 将我的代码放在另一个目录中,最后还保存了剪切图。第一次 运行,我生成了这些切口:'cutout-IMG-HSC-I-18114-0,0-OBJ-NEP175719.1+652743.2.fits',
'cutout-IMG-HSC-I-18114-1,1-OBJ-NEP175509.2+653523.9.fits'
。几分钟后,当我 运行 完全相同的代码而不做任何更改时,我还生成了这两个新图像:'cutout-IMG-HSC-I-18114-0,0-OBJ-NEP175654.7+652930.2.fits',
'cutout-IMG-HSC-I-18114-1,1-OBJ-NEP175458.4+653419.1.fits'
,类似地用于未来的 运行s.
如您所见,它一定没有在我的一张图像中找到任何对象(这很好),但是每个 运行 它都以某种方式在其他每个图像中找到了一个新对象(确实,每次我 运行 这个小子集上的这个代码,我都会看到两个新的剪切图像被保存,具有不同的对象名称)。我很难过,因为就像我说的那样,它搜索的对象和坐标在每个图像中都是相同的 运行。任何见解或猜测将不胜感激!
import astropy
from astropy.nddata.utils import Cutout2D, NoOverlapError
import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u
from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord, Angle
import re
import glob
def make_cutouts(img_file, box_len):
# Image data
hdulist = fits.open(img_file)
img_data = fits.getdata(img_file)
img_name = re.search(r'calexp\-(.+)\.fits', img_file)[1]
# Make cutouts (by locating coords in image with WCS)
wcs = WCS(hdulist[1].header)
for i in range(len(objs)):
# Skip if cutout already exists
if 'cutout-IMG-'+img_name+'-OBJ-'+objs[i]+'.fits' in glob.glob('cutout*.fits'):
print('Cutout of object already exists for this image, skipping...')
continue
# Convert ra, dec to HMS for specific object
ra_h = re.search(r'h=(\d+.?\d*)', str(ras_hms[i]))[1]
ra_m = re.search(r'm=(\d+.?\d*)', str(ras_hms[i]))[1]
ra_s = re.search(r's=(\d+.?\d*)', str(ras_hms[i]))[1]
ra_angle = Angle((float(ra_h), float(ra_m), float(ra_s)), unit='hourangle')
dec_angle = decs_deg[i]
# Coordinate transformation to pixels
center = SkyCoord(ra_angle, dec_angle, frame='fk5')
xp, yp = astropy.wcs.utils.skycoord_to_pixel(center, wcs=wcs, origin=1)
# Make cutout, skip when object is not in image
size = u.Quantity((box_len,box_len),u.arcsec)
try:
co = Cutout2D(img_data,(xp, yp),size,wcs=wcs)
except NoOverlapError:
continue
hdu = fits.PrimaryHDU(data=co.data,header=co.wcs.to_header())
hdu.writeto('cutout-IMG-'+img_name+'-OBJ-'+objs[i]+'.fits', overwrite=True)
return
# Gather all original images
all_images = glob.glob('/Users/myuser/Desktop/Original_Images/calexp*.fits')
coords_file = 'good_dataset.fits'
# Coordinates
hdul = fits.open(coords_file)
coords_data = hdul[1].data
objs = coords_data['Name']
ras = np.array(coords_data['RA']) # in decimal degrees
decs = np.array(coords_data['DEC']) # in decimal degrees
# Convert coordinate systems using astropy
decs_deg = Angle(decs, unit=u.deg)
ras_deg = Angle(ras, unit=u.deg)
ras_hms = [ra.hms for ra in ras_deg]
count=0
for image in all_images:
make_cutouts(image, 10.0)
count+=1
print('Image %d out of %d completed' % (count, len(all_images)))
这里是我刚才的运行打印语句的输出示例,它再次生成了两个新的剪切图像(不同的对象,相同的两个图像)...这里是图像 2没有发现任何物体的地方。此外,有趣的是,每个图像的 "already exists, skipping" 语句的数量增加了两个 运行.
Cutout of object already exists for this image, skipping...
Cutout of object already exists for this image, skipping...
Cutout of object already exists for this image, skipping...
Cutout of object already exists for this image, skipping...
Image 1 out of 3 completed
Image 2 out of 3 completed
Cutout of object already exists for this image, skipping...
Cutout of object already exists for this image, skipping...
Cutout of object already exists for this image, skipping...
Cutout of object already exists for this image, skipping...
Image 3 out of 3 completed
这是一个简单的错误:在 for
循环的末尾有一个 return
语句,这意味着 make_cutouts
的每个 运行 被限制为最多产生一个切口。每次你 运行 它,它都会产生第一个切口,然后下一次,它看到一个存在,用 continue
语句跳过它,然后得到下一个,然后退出。删除 return
语句,代码可能会正常工作。
但是,有几件事您应该尽量避免:
(1) 您在函数中使用了全局变量。您最好将 objs
和 ras_hms
作为参数传递给函数,而不是依赖全局状态来访问它们。
(2) 当您可以遍历对象本身时,您正在遍历索引,即 for thisobj, thisra in zip(objs, ras_hms):
(3) 更小,但 if 'cutout-IMG-'+img_name+'-OBJ-'+objs[i]+'.fits' in glob.glob('cutout*.fits'):
比 if os.path.exists(if 'cutout-IMG-'+img_name+'-OBJ-'+objs[i]+'.fits' in glob.glob('cutout*.fits'):
更有效。如果您使用 'cutout-IMG-{img_name}-OBJ-{obj_id}.fits'.format(img_name=img_name, obj_id=objs[i])
作为字符串