Skip to content

minor typos and pep8 fixes in fits cube notebook #233

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
adrn merged 1 commit into from
Jul 31, 2018
Merged

minor typos and pep8 fixes in fits cube notebook #233

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
177 changes: 90 additions & 87 deletions tutorials/notebooks/FITS-cubes/FITS-cubes.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,7 @@
"\n",
"import astropy.units as u\n",
"from astropy.utils.data import download_file\n",
"from astropy.io import fits # We use fits to open the actual data file\n",
"from astropy.io import fits # We use fits to open the actual data file\n",
"\n",
"import aplpy\n",
"from spectral_cube import SpectralCube\n",
Expand Down Expand Up @@ -90,7 +90,7 @@
},
"outputs": [],
"source": [
"#Downloads the HI data in a fits file format\n",
"# Downloads the HI data in a fits file format\n",
"hi_datafile = download_file(\n",
" 'http://cdsarc.u-strasbg.fr/vizier/ftp/cats/J/A+A/594/A116/CUBES/GAL/TAN/TAN_C14.fits',\n",
" cache=True, show_progress = True)"
Expand All @@ -116,9 +116,9 @@
},
"outputs": [],
"source": [
"hi_data = fits.open(hi_datafile) # Open the FITS file for reading\n",
"cube = SpectralCube.read(hi_data) # Initiate a SpectralCube\n",
"hi_data.close() # Close the FITS file - we already read it in and don't need it anymore!"
"hi_data = fits.open(hi_datafile) # Open the FITS file for reading\n",
"cube = SpectralCube.read(hi_data) # Initiate a SpectralCube\n",
"hi_data.close() # Close the FITS file - we already read it in and don't need it anymore!"
]
},
{
Expand All @@ -130,7 +130,7 @@
"\n",
"\n",
"`\n",
"cube = SpectralCube.read(path_to_data_file/TAN_C14.fits') \n",
"cube = SpectralCube.read('path_to_data_file/TAN_C14.fits') \n",
"`\n",
"</div>\n",
"\n",
Expand Down Expand Up @@ -177,7 +177,7 @@
},
"outputs": [],
"source": [
"cube[600,:,:].quicklook() # Slice the cube along the spectral axis, and display a quick image"
"cube[600, :, :].quicklook() # Slice the cube along the spectral axis, and display a quick image"
]
},
{
Expand All @@ -188,7 +188,7 @@
},
"outputs": [],
"source": [
"cube[:,150,150].quicklook() # Extract a single spectrum through the data cube"
"cube[:, 150, 150].quicklook() # Extract a single spectrum through the data cube"
]
},
{
Expand Down Expand Up @@ -242,8 +242,8 @@
},
"outputs": [],
"source": [
"_, b, _ = cube.world[0,:,0] #extract latitude world coordinates from cube\n",
"_, _, l = cube.world[0,0,:] #extract longitude world coordinates from cube"
"_, b, _ = cube.world[0, :, 0] #extract latitude world coordinates from cube\n",
"_, _, l = cube.world[0, 0, :] #extract longitude world coordinates from cube"
]
},
{
Expand All @@ -262,10 +262,13 @@
"outputs": [],
"source": [
"def find_nearest_idx(array, target_value): \n",
" '''Simple function to find the index closest to a target value'''\n",
" \"\"\"\n",
" Simple function to find the index closest to a target value\n",
" \"\"\"\n",
" idx = np.nanargmin(np.abs(array-target_value))\n",
" return idx\n",
"\n",
"\n",
"# Define desired latitude and longitude range\n",
"lat_range = [-46, -40] * u.deg \n",
"lon_range = [306, 295] * u.deg\n",
Expand All @@ -276,7 +279,7 @@
"lon_range_idx = sorted([find_nearest_idx(l, lon_range[0]), find_nearest_idx(l, lon_range[1])])\n",
"\n",
"# Create a sub_cube cut to these coordinates\n",
"sub_cube = cube[:,lat_range_idx[0]:lat_range_idx[1], lon_range_idx[0]:lon_range_idx[1]]\n",
"sub_cube = cube[:, lat_range_idx[0]:lat_range_idx[1], lon_range_idx[0]:lon_range_idx[1]]\n",
"\n",
"print(sub_cube)"
]
Expand Down Expand Up @@ -325,15 +328,15 @@
},
"outputs": [],
"source": [
"moment_0 = sub_cube_slab.with_spectral_unit(u.km/u.s).moment(order = 0) # Zero-th moment \n",
"moment_1 = sub_cube_slab.with_spectral_unit(u.km/u.s).moment(order = 1) # First moment\n",
"moment_0 = sub_cube_slab.with_spectral_unit(u.km/u.s).moment(order=0) # Zero-th moment \n",
"moment_1 = sub_cube_slab.with_spectral_unit(u.km/u.s).moment(order=1) # First moment\n",
"\n",
"# Write the moments as a FITS image\n",
"#moment_0.write('hi_moment_0.fits') \n",
"#moment_1.write('hi_moment_1.fits')\n",
"# moment_0.write('hi_moment_0.fits') \n",
"# moment_1.write('hi_moment_1.fits')\n",
"\n",
"print('Moment_0 has units of: ',moment_0.unit)\n",
"print('Moment_1 has units of: ',moment_1.unit)\n",
"print('Moment_0 has units of: ', moment_0.unit)\n",
"print('Moment_1 has units of: ', moment_1.unit)\n",
"\n",
"# Convert Moment_0 to a Column Density assuming optically thin media\n",
"hi_column_density = moment_0 * 1.82 * 10**18 / (u.cm * u.cm) * u.s / u.K / u.km"
Expand Down Expand Up @@ -372,28 +375,28 @@
"outputs": [],
"source": [
"# Initiate a figure \n",
"fig = plt.figure(figsize = (18,12))\n",
"fig = plt.figure(figsize=(18, 12))\n",
"\n",
"# Initiate a FITSFigure to set up axes\n",
"F = aplpy.FITSFigure(moment_1.hdu, figure = fig)\n",
"F = aplpy.FITSFigure(moment_1.hdu, figure=fig)\n",
"\n",
"# Extract the axis object that was created for future manipulation\n",
"ax = fig.gca()\n",
"\n",
"# display a colorscale map of moment_1\n",
"F.show_colorscale(cmap = 'RdBu_r', vmin = 0., vmax = 200.)\n",
"F.show_colorscale(cmap='RdBu_r', vmin=0., vmax=200.)\n",
"# display a colorbar\n",
"F.show_colorbar(axis_label_text = 'Velocity (km / s)')\n",
"F.show_colorbar(axis_label_text='Velocity (km / s)')\n",
"\n",
"# overplot contours of hi_column_density (essentially column density here)\n",
"F.show_contour(hi_column_density.hdu, cmap = 'Greys_r', levels = (1e20, 5e20, 1e21, 3e21, 5e21, 7e21, 1e22))\n",
"F.show_contour(hi_column_density.hdu, cmap='Greys_r', levels=(1e20, 5e20, 1e21, 3e21, 5e21, 7e21, 1e22))\n",
"\n",
"ax.yaxis.set_tick_params(labelsize = 16)\n",
"ax.xaxis.set_tick_params(labelsize = 16)\n",
"ax.yaxis.set_tick_params(labelsize=16)\n",
"ax.xaxis.set_tick_params(labelsize=16)\n",
"x_lab = ax.get_xlabel()\n",
"y_lab = ax.get_ylabel()\n",
"ax.set_xlabel(x_lab, fontsize = 16)\n",
"ax.set_ylabel(x_lab, fontsize = 16)\n"
"ax.set_xlabel(x_lab, fontsize=16)\n",
"ax.set_ylabel(x_lab, fontsize=16)\n"
]
},
{
Expand All @@ -416,7 +419,7 @@
},
"outputs": [],
"source": [
"print(moment_1.wcs) # Examine the WCS object associated with the moment map"
"print(moment_1.wcs) # Examine the WCS object associated with the moment map"
]
},
{
Expand All @@ -435,31 +438,31 @@
"outputs": [],
"source": [
"# Initiate a figure and axis object with WCS projection information\n",
"fig = plt.figure(figsize = (18,12))\n",
"ax = fig.add_subplot(111,projection = moment_1.wcs)\n",
"fig = plt.figure(figsize=(18, 12))\n",
"ax = fig.add_subplot(111, projection=moment_1.wcs)\n",
"\n",
"# Display the moment map image\n",
"im = ax.imshow(moment_1.hdu.data, cmap = 'RdBu_r', vmin = 0, vmax = 200)\n",
"ax.invert_yaxis() # Flips the Y axis \n",
"im = ax.imshow(moment_1.hdu.data, cmap='RdBu_r', vmin=0, vmax=200)\n",
"ax.invert_yaxis() # Flips the Y axis \n",
"\n",
"# Add axes labels\n",
"ax.set_xlabel(\"Galactic Longitude (degrees)\", fontsize = 16)\n",
"ax.set_ylabel(\"Galactic Latitude (degrees)\", fontsize = 16)\n",
"ax.set_xlabel(\"Galactic Longitude (degrees)\", fontsize=16)\n",
"ax.set_ylabel(\"Galactic Latitude (degrees)\", fontsize=16)\n",
"\n",
"# Add a colorbar\n",
"cbar = plt.colorbar(im, pad = .07)\n",
"cbar.set_label('Velocity (km/s)', size = 16)\n",
"cbar = plt.colorbar(im, pad=.07)\n",
"cbar.set_label('Velocity (km/s)', size=16)\n",
"\n",
"# Overlay set of RA/Dec Axes\n",
"overlay = ax.get_coords_overlay('fk5')\n",
"overlay.grid(color='white', ls='dotted', lw = 2)\n",
"overlay[0].set_axislabel('Right Ascension (J2000)', fontsize = 16)\n",
"overlay[1].set_axislabel('Declination (J2000)', fontsize = 16)\n",
"overlay.grid(color='white', ls='dotted', lw=2)\n",
"overlay[0].set_axislabel('Right Ascension (J2000)', fontsize=16)\n",
"overlay[1].set_axislabel('Declination (J2000)', fontsize=16)\n",
"\n",
"# Overplot column density contours \n",
"levels = (1e20, 5e20, 1e21, 3e21, 5e21, 7e21, 1e22) # Define contour levels to use\n",
"ax.contour(hi_column_density.hdu.data, cmap = 'Greys_r', alpha = 0.5, \n",
" lw = 3, levels = levels)"
"levels = (1e20, 5e20, 1e21, 3e21, 5e21, 7e21, 1e22) # Define contour levels to use\n",
"ax.contour(hi_column_density.hdu.data, cmap='Greys_r', alpha=0.5, \n",
" lw=3, levels=levels)"
]
},
{
Expand Down Expand Up @@ -495,25 +498,25 @@
},
"outputs": [],
"source": [
"lat_slice = 35 # Index of latitude dimension to slice along\n",
"lat_slice = 35 # Index of latitude dimension to slice along\n",
"\n",
"# Initiate a figure and axis object with WCS projection information\n",
"fig = plt.figure(figsize = (18,12))\n",
"ax = fig.add_subplot(111,projection = sub_cube_slab.wcs, slices = ('y', lat_slice, 'x'))\n",
"fig = plt.figure(figsize=(18, 12))\n",
"ax = fig.add_subplot(111, projection=sub_cube_slab.wcs, slices=('y', lat_slice, 'x'))\n",
"# Above, we have specified to plot the longitude along the y axis, pick just the lat_slice indicated, \n",
"# and plot the velocity along the x axis\n",
"\n",
"# Display the slice\n",
"im = ax.imshow(sub_cube_slab.hdu.data[:,lat_slice,:].transpose()) # Display the image slice\n",
"ax.invert_yaxis() # Flips the Y axis \n",
"im = ax.imshow(sub_cube_slab.hdu.data[:,lat_slice,:].transpose()) # Display the image slice\n",
"ax.invert_yaxis() # Flips the Y axis \n",
"\n",
"# Add axes labels\n",
"ax.set_xlabel(\"LSR Velocity (m/s)\", fontsize = 16)\n",
"ax.set_ylabel(\"Galactic Longitude (degrees)\", fontsize = 16)\n",
"ax.set_xlabel(\"LSR Velocity (m/s)\", fontsize=16)\n",
"ax.set_ylabel(\"Galactic Longitude (degrees)\", fontsize=16)\n",
"\n",
"# Add a colorbar\n",
"cbar = plt.colorbar(im, pad = .07, orientation = 'horizontal')\n",
"cbar.set_label('Temperature (K)', size = 16)\n",
"cbar = plt.colorbar(im, pad=.07, orientation='horizontal')\n",
"cbar.set_label('Temperature (K)', size=16)\n",
"\n"
]
},
Expand Down Expand Up @@ -579,7 +582,7 @@
"outputs": [],
"source": [
"# Query for Herschel data in a 1 degree radius around the SMC\n",
"result = ESASky.query_region_maps('SMC', radius = 1*u.deg, missions = 'Herschel')\n",
"result = ESASky.query_region_maps('SMC', radius=1*u.deg, missions='Herschel')\n",
"\n",
"print(result)"
]
Expand Down Expand Up @@ -635,16 +638,16 @@
},
"outputs": [],
"source": [
"filters = result['HERSCHEL']['filter'].astype(str) # Convert the list of filters from the query to a string\n",
"filters = result['HERSCHEL']['filter'].astype(str) # Convert the list of filters from the query to a string\n",
"\n",
"# Construct a boolean mask, searching for only the desired filters\n",
"mask = np.array(['250, 350, 500' == s for s in filters], dtype = 'bool')\n",
"mask = np.array(['250, 350, 500' == s for s in filters], dtype='bool')\n",
"\n",
"# Re-construct a new TableList object containing only our desired query entry\n",
"target_obs = TableList({\"HERSCHEL\":result['HERSCHEL'][mask]}) # This will be passed into ESASky.get_maps()\n",
"target_obs = TableList({\"HERSCHEL\":result['HERSCHEL'][mask]}) # This will be passed into ESASky.get_maps()\n",
"\n",
"IR_images = ESASky.get_maps(target_obs) # Download the images\n",
"IR_images['HERSCHEL'][0]['350'].info() # Display some information about the 350 micron image"
"IR_images = ESASky.get_maps(target_obs) # Download the images\n",
"IR_images['HERSCHEL'][0]['350'].info() # Display some information about the 350 micron image"
]
},
{
Expand All @@ -665,8 +668,8 @@
"outputs": [],
"source": [
"herschel_header = IR_images['HERSCHEL'][0]['350']['image'].header\n",
"herschel_wcs = WCS(IR_images['HERSCHEL'][0]['350']['image']) # Extract WCS information\n",
"herschel_imagehdu = IR_images['HERSCHEL'][0]['350']['image'] # Extract Image data\n",
"herschel_wcs = WCS(IR_images['HERSCHEL'][0]['350']['image']) # Extract WCS information\n",
"herschel_imagehdu = IR_images['HERSCHEL'][0]['350']['image'] # Extract Image data\n",
"print(herschel_wcs)"
]
},
Expand All @@ -686,28 +689,28 @@
"outputs": [],
"source": [
"# Initiate a figure and axis object with WCS projection information\n",
"fig = plt.figure(figsize = (18,12))\n",
"ax = fig.add_subplot(111,projection = herschel_wcs)\n",
"fig = plt.figure(figsize=(18, 12))\n",
"ax = fig.add_subplot(111, projection=herschel_wcs)\n",
"\n",
"# Display the moment map image\n",
"im = ax.imshow(herschel_imagehdu.data, cmap = 'viridis', \n",
" norm = LogNorm(), vmin = 2, vmax = 50)\n",
"#ax.invert_yaxis() # Flips the Y axis \n",
"im = ax.imshow(herschel_imagehdu.data, cmap='viridis', \n",
" norm=LogNorm(), vmin=2, vmax=50)\n",
"# ax.invert_yaxis() # Flips the Y axis \n",
"\n",
"# Add axes labels\n",
"ax.set_xlabel(\"Right Ascension\", fontsize = 16)\n",
"ax.set_ylabel(\"Declination\", fontsize = 16)\n",
"ax.grid(color = 'white', ls = 'dotted', lw = 2)\n",
"\n",
"# Add a colorbar\n",
"cbar = plt.colorbar(im, pad = .07)\n",
"cbar = plt.colorbar(im, pad=.07)\n",
"cbar.set_label(''.join(['Herschel 350'r'$\\mu$m ','(', herschel_header['BUNIT'], ')']), size = 16)\n",
"\n",
"# Overlay set of Galactic Coordinate Axes\n",
"overlay = ax.get_coords_overlay('galactic') \n",
"overlay.grid(color='black', ls='dotted', lw = 1)\n",
"overlay[0].set_axislabel('Galactic Longitude', fontsize = 14)\n",
"overlay[1].set_axislabel('Galactic Latitude', fontsize = 14)"
"overlay.grid(color='black', ls='dotted', lw=1)\n",
"overlay[0].set_axislabel('Galactic Longitude', fontsize=14)\n",
"overlay[1].set_axislabel('Galactic Latitude', fontsize=14)"
]
},
{
Expand All @@ -730,46 +733,46 @@
"outputs": [],
"source": [
"# Initiate a figure and axis object with WCS projection information\n",
"fig = plt.figure(figsize = (18,12))\n",
"ax = fig.add_subplot(111,projection = herschel_wcs)\n",
"fig = plt.figure(figsize=(18, 12))\n",
"ax = fig.add_subplot(111, projection=herschel_wcs)\n",
"\n",
"# Display the moment map image\n",
"im = ax.imshow(herschel_imagehdu.data, cmap = 'viridis', \n",
" norm = LogNorm(), vmin = 5, vmax = 50, alpha = .8)\n",
"#ax.invert_yaxis() # Flips the Y axis \n",
"im = ax.imshow(herschel_imagehdu.data, cmap='viridis', \n",
" norm=LogNorm(), vmin=5, vmax=50, alpha=.8)\n",
"# ax.invert_yaxis() # Flips the Y axis \n",
"\n",
"# Add axes labels\n",
"ax.set_xlabel(\"Right Ascension\", fontsize = 16)\n",
"ax.set_ylabel(\"Declination\", fontsize = 16)\n",
"ax.grid(color = 'white', ls = 'dotted', lw = 2)\n",
"ax.set_xlabel(\"Right Ascension\", fontsize=16)\n",
"ax.set_ylabel(\"Declination\", fontsize=16)\n",
"ax.grid(color = 'white', ls='dotted', lw=2)\n",
"\n",
"# Extract x and y coordinate limits\n",
"x_lim = ax.get_xlim()\n",
"y_lim = ax.get_ylim()\n",
"\n",
"# Add a colorbar\n",
"cbar = plt.colorbar(im, fraction=0.046, pad=-0.1)\n",
"cbar.set_label(''.join(['Herschel 350'r'$\\mu$m ','(', herschel_header['BUNIT'], ')']), size = 16)\n",
"cbar.set_label(''.join(['Herschel 350'r'$\\mu$m ','(', herschel_header['BUNIT'], ')']), size=16)\n",
"\n",
"# Overlay set of RA/Dec Axes\n",
"overlay = ax.get_coords_overlay('galactic')\n",
"overlay.grid(color='black', ls='dotted', lw = 1)\n",
"overlay[0].set_axislabel('Galactic Longitude', fontsize = 14)\n",
"overlay[1].set_axislabel('Galactic Latitude', fontsize = 14)\n",
"overlay.grid(color='black', ls='dotted', lw=1)\n",
"overlay[0].set_axislabel('Galactic Longitude', fontsize=14)\n",
"overlay[1].set_axislabel('Galactic Latitude', fontsize=14)\n",
"\n",
"hi_transform = ax.get_transform(hi_column_density.wcs) # extract axes Transform information for the HI data\n",
"hi_transform = ax.get_transform(hi_column_density.wcs) # extract axes Transform information for the HI data\n",
"\n",
"# Overplot column density contours \n",
"levels = (2e21, 3e21, 5e21, 7e21, 8e21, 1e22) # Define contour levels to use\n",
"ax.contour(hi_column_density.hdu.data, cmap = 'Greys_r', alpha = 0.8, lw = 5, levels = levels,\n",
" transform = hi_transform) # include the transform information with the keyword \"transform\"\n",
"levels = (2e21, 3e21, 5e21, 7e21, 8e21, 1e22) # Define contour levels to use\n",
"ax.contour(hi_column_density.hdu.data, cmap='Greys_r', alpha=0.8, lw=5, levels=levels,\n",
" transform=hi_transform) # include the transform information with the keyword \"transform\"\n",
"\n",
"# Overplot velocity image so we can also see the Gas velocities\n",
"im_hi = ax.imshow(moment_1.hdu.data, cmap = 'RdBu_r', vmin = 0, vmax = 200, alpha = 0.5, transform = hi_transform)\n",
"im_hi = ax.imshow(moment_1.hdu.data, cmap='RdBu_r', vmin=0, vmax=200, alpha=0.5, transform=hi_transform)\n",
"\n",
"# Add a second colorbar for the HI Velocity information\n",
"cbar_hi = plt.colorbar(im_hi, orientation = 'horizontal', fraction=0.046, pad=0.07)\n",
"cbar_hi.set_label('HI 'r'$21$cm Mean Velocity (km/s)', size = 16)\n",
"cbar_hi = plt.colorbar(im_hi, orientation='horizontal', fraction=0.046, pad=0.07)\n",
"cbar_hi.set_label('HI 'r'$21$cm Mean Velocity (km/s)', size=16)\n",
"\n",
"# Apply original image x and y coordinate limits\n",
"ax.set_xlim(x_lim)\n",
Expand Down