from astropy.io import fits
import matplotlib.pyplot as plt 
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
from astropy.nddata import Cutout2D
import numpy as np

hst = fits.open('hst_10596_02_acs_wfc_f814w_drz.fits')
legacy = fits.open('legacysurvey-0422m082-image-z.fits.fz')

f = plt.figure()
position = SkyCoord(42.30362*u.deg, -8.25714*u.deg, frame='icrs')
wcs1 = WCS(hst[1].header)
size = u.Quantity((25, 25), u.arcsec)

hst_cutout = Cutout2D(hst[1].data, position, size, wcs=wcs1)
f.add_subplot(1,2,1)
plt.imshow(hst_cutout.data, Norm=LogNorm())
f.add_subplot(1,2,2)
wcs2 = WCS(legacy[1].header)
legacy_cutout = Cutout2D(legacy[1].data, position, size, wcs=wcs2.celestial)
plt.imshow(legacy_cutout.data)
plt.show()