Hello folks,
I have a script that doesn’t output all the images. Normally it should look at the FITS header, work out RA/Dec, center it, take snaps of the galaxies in the background, annotate them and combine the two pictures into one.
After running the code it opens Figure 1 and Figure 2, as well as saving img1.png. It detected the galaxies as well but it doesn’t use the snippets to create img2 and since img2 wasn’t created, it can’t rescale the two images into one.
Maybe some of you can help me because I am at a loss
# -*- coding: utf-8 -*-
"""M81M82_visualization.ipynb
Automatically generated by Colaboratory.
Original file is located at
https://colab.research.google.com/drive/1BmtQxdByqfsm6iVLxr8s-5eEXNfVHOgw
# Load fits image and header
## Attention: in this case only the header of the fits file is used. the image originates from a png!!!
"""
from astropy.io import fits
from astropy.coordinates import SkyCoord
from astropy import coordinates as coord
from astropy.wcs.utils import skycoord_to_pixel
from astropy.table import Table
import astropy.units as u
from astropy.wcs import WCS
from astroquery.simbad import Simbad
import matplotlib.pyplot as plt
from astropy.visualization import astropy_mpl_style
import numpy as np
plt.style.use('dark_background')
# Load the FITS image and extract header information
image_file = "C:/Users/xDDra/Desktop/Annotation/M81.fit"
raw_image_file = None#'M81.png'
hdul = fits.open(image_file)
# use better quality png in this case
if not raw_image_file:
data = hdul[0].data
img = np.swapaxes(np.swapaxes(data,0,2),0,1).astype(float)
img = img/max(img.flatten())
else:
img = plt.imread(raw_image_file)[::-1]
# get fits header
header = hdul[0].header
# get world coordinates system and drop color channel
wcs = WCS(header).dropaxis(2)
# Extract the right ascension and declination of the center of the image
center_ra = header['CRVAL1']
center_dec = header['CRVAL2']
fig = plt.figure(figsize=(10,6))
ax = plt.subplot(projection=wcs, label='overlays')
ax.imshow(img, origin='lower', cmap='gray')
ax.coords.grid(True, color='white', ls='-', alpha=.5)
ax.coords[0].set_axislabel('Right Ascension (J2000)')
ax.coords[1].set_axislabel('Declination (J2000)')
ax.scatter(center_ra, center_dec, transform=ax.get_transform('icrs'),s=10, edgecolor='red', facecolor='red')
img.shape
"""# Query Objects
### 1. define the center of the circle
### 2. define the search header using header-information
### 3. query objects
### 4. filter objects outside the FOV
### 5. fitler objects according to catalog (see definition in column TYPE)
"""
Simbad.TIMEOUT = 120
# Define the specific coordinates
center_ra = header['CRVAL1']
center_dec = header['CRVAL2']
target_coord = SkyCoord(ra=center_ra, dec=center_dec, unit=(u.deg, u.deg), frame='icrs')
# define the search radius
radius_deg = max([abs(header['CDELT1']),abs(header['CDELT2'])])*max([header['NAXIS1'],header['NAXIS2']])/1.5
# Set up the query criteria
custom_query = f"region(CIRCLE, {target_coord.ra.deg} + {target_coord.dec.deg}, {radius_deg}d)"
# Query Simbad with the criteria
result_table = Simbad.query_criteria(custom_query, otype='galaxy')
df = result_table.to_pandas()
# filter by position
df['Pixel_Position'] = [skycoord_to_pixel(SkyCoord(v[0],v[1], unit=(u.hourangle, u.deg)), wcs) for v in df[['RA','DEC']].values]
df['px'] = df['Pixel_Position'].apply(lambda x: int(x[0]))
df['py'] = df['Pixel_Position'].apply(lambda x: int(x[1]))
df = df[(df.px>0)&(df.py>0)&(df.px<img.shape[1])&(df.py<img.shape[0])]
# filter by catalog
# get catalog/TYPE of object by simple string manipulation
df['TYPE'] = df['MAIN_ID'].apply(lambda x: x.split(' ')[0].split('+')[0])#.value_counts()
filter_types = ['M', 'NGC', 'IC', 'MCG', 'UGC', '2MASX']
filtered_result = Table.from_pandas(df[df['TYPE'].isin(filter_types)])
filtered_result
df[~df['TYPE'].isin(filter_types)].TYPE.value_counts()
"""# Plot first image
## Attention: here I defined a dictionary for sorting the results accordingly
"""
img.shape
import matplotlib.patches as patches
sub=1
dpi=150
dfi = filtered_result.to_pandas()
rep_dic ={
'M':1,'NGC':2,'IC':3, 'UGC':4#, 'MCG', '2MASS', '2MASX', 'LEDA'
}
dfi['sorting'] = dfi.TYPE.replace(rep_dic)
dfi.MAIN_ID = dfi.MAIN_ID.apply(lambda x: x.replace(' ',''))
dfi = dfi.sort_values(['sorting', 'MAIN_ID'], ascending=True).reset_index()
scale=20
fig = plt.figure(figsize=(scale,(img.shape[0]/img.shape[1])*scale))
ax1 = plt.subplot(projection=wcs, label='overlays')
ax1.imshow(img[::sub,::sub])
ax1.coords.grid(True, color='white', ls='-', alpha=.5)
ax1.coords[0].set_axislabel('Right Ascension (J2000)')
ax1.coords[1].set_axislabel('Declination (J2000)',minpad=-1)
all_patches = []
cmap = plt.get_cmap("tab10")
filter_idxs = []
j=0
for i,row in dfi.iterrows():
if row.MAIN_ID == 'M81':
patch_size=1024
color = cmap(0)
lw=4
alpha=1
elif row.TYPE == 'M':
patch_size=512
color = cmap(0)
lw=4
alpha=1
elif row.TYPE == 'NGC':
patch_size=256
color = cmap(2)
lw=4
alpha=1
elif row.TYPE == 'IC':
patch_size=128
color = cmap(1)
lw=4
alpha=1
elif row.TYPE == 'UGC':
patch_size=128
color = cmap(1)
lw=4
alpha=1
else:
patch_size=32
color = 'white'
lw=2
alpha=.5
lw=1
alpha=.5
color = 'white'
if (row.px-patch_size//2 > 0)&(row.px+patch_size//2 < img.shape[1])&(row.py-patch_size//2 > 0)&(row.py+patch_size//2 < img.shape[0]):
rect = patches.Rectangle((row.px-patch_size//2, row.py-patch_size//2), patch_size, patch_size, alpha=alpha, linewidth=lw, edgecolor=color, facecolor='none')
ax1.add_patch(rect)
ax1.text(row.px-patch_size//2,row.py-patch_size//2-10,str(j+1),#+' ['+row.MAIN_ID+']',
ha='left',va='top',color=color,fontsize=14)
j+=1
patch = img[row.py-patch_size//2:row.py+patch_size//2,row.px-patch_size//2:row.px+patch_size//2]
all_patches.append(patch)
filter_idxs.append(i)
#ax1.set_title("M81 Bodes Galaxy & M81 Cigar Galaxy", fontsize=18)
#points = dfi[dfi.MAIN_ID.isin(['M84','M86','NGC4477','NGC4473','NGC4461','NGC4458','NGC4438','NGC4435'])][['px','py']].values
#ax1.plot(points[:, 0], points[:, 1], 'w',lw=3, ls='--',alpha=.5)
plt.tight_layout()
plt.savefig('C:/Users/xDDra/Desktop/Annotation/img1.png', bbox_inches='tight', pad_inches=.1, dpi=dpi)
plt.show()
img1 = plt.imread('C:/Users/xDDra/Desktop/Annotation/img1.png')
all_patches = np.array(all_patches)
all_patches.shape#, img1.shape
"""# darkmatters Logo"""
logo = plt.imread('C:/Users/xDDra/Desktop/Annotation/darkmatters_logo.png')
plt.imshow(logo)
"""# Second Image"""
m=6#int(np.ceil(np.sqrt(len(all_patches))))
n=10#10#int(np.ceil(len(all_patches)/m))
print(len(all_patches), m, n, m*n)
dft = dfi.iloc[filter_idxs].reset_index()
scale=2
fig, axarr = plt.subplots(m,n,figsize=(n*scale,m*scale))
for i, row in dft.iterrows():
ax = axarr[i//n,i%n]
ax.imshow(all_patches[i][::-1])
ax.set_title(row.MAIN_ID, fontsize=10)
ax.set_xticks([])
ax.set_yticks([])
ax.text(2,4,str(i+1),ha='left',va='top',fontsize=12)
for i in np.arange(len(all_patches),m*n):
ax = axarr[i//n,i%n]
ax.axis('off')
print()
#axarr[m-1,n-1].imshow(logo1, cmap='Greys')
#axarr[m-1,n-2].imshow(logo2, cmap='gray')
axarr[m-1,n-1].imshow(logo3)
plt.tight_layout()
plt.savefig('C:/Users/xDDra/Desktop/Annotation/img2.png', bbox_inches='tight', pad_inches=.1, dpi=dpi)
plt.show()
img2 = plt.imread('C:/Users/xDDra/Desktop/Annotation/img2.png')
img2.shape
"""# Image rescaling s.t. vstack works
## ugly, I know, but it does the job.
"""
from skimage.transform import resize
from PIL import Image
img22 = (resize(img2, (img2.shape[0]*(img1.shape[1]/img2.shape[1]), img1.shape[1]))*255).astype(np.uint8)
im = Image.fromarray(np.vstack([(img1*255).astype(np.uint8),img22])[:,:,:3])
im.save("C:/Users/xDDra/Desktop/Annotation/Annotated", quality=95, subsampling=0)