Commit 334826bc authored by karpov-sv's avatar karpov-sv
Browse files

Some documentation on smaller but still convenient utility routines added

parent a0076109
......@@ -34,6 +34,7 @@ Next, you will need to import the modules from *STDPipe* itself:
Now you have everything imported, and may start actual data analysis!
Data processing
---------------
......@@ -49,6 +50,8 @@ Data processing
transients
cutouts
psf
utils
Common principles
-----------------
......@@ -96,3 +99,12 @@ Common conventions for routine arguments:
- Functions that operate on temporary files in temporary folders may be supplied with `_workdir` argument - then they will store all temporary files related to their work in this specific folder, and will not remove them afterwards. So e.g. you will be able to directly see e.g what configuration files were created, manually re-run the failed command to experiment with its options (with `verbose=True` function call typically prints the complete command line of all calls to external programs, so you may just directly copy and paste it to terminal to repeat its invocation), etc.
- Functions that run external programs (e.g. SExtractor, HOTPANTS, or Astrometry.Net) usually accept `_exe` argument to directly specify the path to corresponding executable. If not specified, the code will try to automatically find it for you, so normally you do not need to worry about it.
Examples
--------
For some examples of actual usage of *STDPipe* routines we encourage you to check the `example notebooks <https://github.com/karpov-sv/stdpipe/blob/master/notebooks>`__, especially the `STDPipe tutorial <https://github.com/karpov-sv/stdpipe/blob/master/notebooks/stdpipe_tutorial.ipynb>`__ that demonstrates all basic steps of a typical image processing.
We also have the `tutorial on injection of simulated stars <https://github.com/karpov-sv/stdpipe/blob/master/notebooks/simulated_stars.ipynb>`__ into the image and simple analysis of object detection efficiency.
`One more tutorial <https://github.com/karpov-sv/stdpipe/blob/master/notebooks/image_stacking.ipynb>`__ shows how to do a simple image co-addition (stacking) with astrometric alignment.
Convenience utilities
=====================
*STDPipe* also contains a number of convenience utilities that help solving some of maybe minor but still quite important problems you may encounter during data analysis and visualization.
Splitting the image into sub-images
-----------------------------------
Often you have quite large image, with significant variation of PSF, background etc over it. Such frames are not very well suited for image subtraction, as convolution kernel will become unstable and / or unable to properly match the resolutions of image and template.
So, for such images, it is better to split them into smaller pieces and do image subtraction on the latters. When splitting, however, you have to also split the mask, update WCS solution and FITS header, maybe select a subset of catalogue stars, etc etc. And we have a dedicated routine for doing just that!
.. code-block:: python
# Crop the image (mask, header, WCS, object list and catalogue too) into 4 (2x2) pieces,
# adding 10 pixels wide overlap between them.
for i, x0, y0, image1, mask1, header1, wcs1, obj1, cat1 in pipeline.split_image(image, nx=2,
mask=mask, header=header, wcs=wcs, obj=obj, cat=cat, overlap=10,
get_index=True, get_origin=True, verbose=True):
# We got a sub-image
print('Subimage', i, 'has origin at', x0, y0, 'shape', image1.shape, 'and contains', len(obj1),
'objects from original image')
# Do something useful on the sub-images here!
pass
.. autofunction:: stdpipe.pipeline.split_image
:noindex:
Extracting the time from FITS header
------------------------------------
We have a routine that helps extracting the information on time of observations from FITS headers.
.. autofunction:: stdpipe.utils.get_obs_time
:noindex:
Displaying the images
---------------------
We have a convenience image plotting function :func:`stdpipe.plots.imshow` what may be used as a drop-in replacement for :func:`matplotlib.pyplot.imshow`, but also supports percentile-based intensity min/max levels, various intensity stretching options (linear, asinh, log, etc), built-in colorbar that plays nicely with sub-plots, and an option to hide the axes around the image.
.. code-block:: python
# Plot the image with asinh intensity scaling between [0.5, 99.5] quantiles, and a color bar
plots.imshow(image, qq=[0.5, 99.5], stretch='asinh', show_colorbar=True)
.. autofunction:: stdpipe.plots.imshow
:noindex:
Displaying 2d histograms of randomly distributed points
-------------------------------------------------------
We also have a convenience function for plotting the images of 2d histograms of irregularly spaced data points. It may shows various statistics of the point values - means, medians, standard deviations etc, as well as display the color bar for the values, do percentile-based intensity scaling, and overlay the positions of the data points onto the histogram.
.. code-block:: python
# Show the map of FWHM of detected objects using 8x8 bins and overlaying the positions
# of the objects. Also, force the image to have correct aspect ratio
plots.binned_map(obj['x'], obj['y'], obj['fwhm'], cmap='hot', bins=8, aspect='equal',
show_colorbar=True, show_dots=True, color='blue')
.. autofunction:: stdpipe.plots.binned_map
:noindex:
......@@ -17,7 +17,9 @@ requirements = [
'photutils',
'statsmodels',
'tqdm',
'regions'
'regions',
'dateutil',
'requests',
]
setup(
......
......@@ -436,12 +436,40 @@ def place_random_stars(image, psf_model, nstars=100, minflux=1, maxflux=100000,
return cat
def split_image(image, nx=1, ny=None, mask=None, bg=None, err=None, header=None, wcs=None, obj=None, cat=None, overlap=0, get_origin=False, verbose=False):
def split_image(image, nx=1, ny=None, mask=None, bg=None, err=None, header=None, wcs=None, obj=None, cat=None, overlap=0, get_index=False, get_origin=False, verbose=False):
"""
Generator to split the image into several (nx x ny) blocks, while also optionally providing the mask, header, wcs and object list for the sub-blocks.
The blocks may optionally be extended by 'overlap' pixels in all directions.
Generator function to split the image into several (`nx` x `ny`) blocks, while also optionally providing the mask, header, wcs, object list etc for the sub-blocks.
The blocks may optionally be extended by 'overlap' pixels in all directions, so that at least in some sub-images every part of original image is far from the edge. This parameter may be used e.g. in conjunction with `edge` parameter of :func:`stdpipe.photometry.get_objects_sextractor` to avoid detecting the same object twice.
:param image: Image to split
:param nx: Number of sub-images in `x` direction
:param ny: Number of sub-images in `y` direction
:param mask: Mask image to split, optional
:param bg: Background map to split, optional
:param err: Noise model image to split, optional
:param header: Image header, optional. If set, the header corresponding to splitted sub-image will be returned, with correctly adjusted WCS information
:param wcs: WCS solution for the image, optional. If set, the solution for sub-image will be returned
:param obj: Object list, optional. If provided, the list of objects contained in the sub-image, with accordingly adjusted pixel coordinates, will be returned
:param cat: Reference catalogue, optional. If provided, the catalogue for the stars on the sub-image will be returned
:param overlap: If set, defines how much sub-images will overlap, in pixels.
:param get_index: If set, also returns the number of current sub-image, starting from zero
:param get_origin: If set, also return the sub-image origin pixel coordinates
:param verbose: Whether to show verbose messages during the run of the function or not. May be either boolean, or a `print`-like function
:returns: Every concecutive call to the generator will return the list of cropped objects corresponding to the next sub-image, as well as some sub-image metadata.
The returned list is constructed from the following elements:
- Index of current sub-image, if :code:`get_index=True`
- `x` and `y` coordinates of current sub-image origin inside the original image
- Current sub-image
- Cropped mask corresponding to current sub-image, if `mask` is provided
- Cropped background map, if `bg` is provided
- Cropped noise model, if `err` is provided
- FITS header corresponding to the sub-image with correct astrometric solution for it, if `header` is provided
- WCS astrometric solution for the sub-image, if `wcs` is provided
- Object list containing only objects that are inside the sub-image with their pixel coordinates adjusted correspondingly, if `obj` is set
- Reference catalogue containing only stars overlaying the sub-image, if `cat` is provided and `wcs` is set
Returns the list consisting of origin x,y coordinates (if get_origin is True), the cropped image, and cropped mask, header, wcs, object list, and catalogie, if they are provided.
"""
# Simple wrapper around print for logging in verbose mode only
......@@ -467,7 +495,14 @@ def split_image(image, nx=1, ny=None, mask=None, bg=None, err=None, header=None,
else:
image1,header1 = _,None
result = [x1, y1] if get_origin else []
result = []
if get_index:
result += [i]
if get_origin:
result += [x1, y1]
result += [image1]
if mask is not None:
......
......@@ -30,7 +30,7 @@ def colorbar(obj=None, ax=None, size="5%", pad=0.1):
ax.get_figure().sca(ax)
def imshow(image, qq=None, show_colorbar=True, show_axis=True, stretch='linear', ax=None, **kwargs):
"""Simple wrapper around pyplot.imshow with histogram-based intensity scaling, optional colorbar, etc.
"""Simple wrapper around pyplot.imshow with percentile-based intensity scaling, optional colorbar, etc.
:param image: Numpy 2d array to display
:param qq: two-element tuple (or list) with quantiles that define lower and upper limits for image intensity normalization. Default is `[0.5, 99.5]`. Will be superseded by manually provided `vmin` and `vmax` arguments.
......@@ -82,7 +82,7 @@ def imshow(image, qq=None, show_colorbar=True, show_axis=True, stretch='linear',
return img
def binned_map(x, y, value, bins=16, statistic='mean', qq=[0.5, 97.5], show_colorbar=True, show_axis=True, show_dots=False, ax=None, **kwargs):
def binned_map(x, y, value, bins=16, statistic='mean', qq=[0.5, 97.5], color=None, show_colorbar=True, show_axis=True, show_dots=False, ax=None, **kwargs):
"""Plots various statistical estimators binned onto regular grid from the set of irregular data points (`x`, `y`, `value`).
:param x: Abscissae of the data points
......@@ -91,6 +91,7 @@ def binned_map(x, y, value, bins=16, statistic='mean', qq=[0.5, 97.5], show_colo
:param bins: Number of bins per axis
:param statistic: Statistical estimator to plot, may be `mean`, `median`, or a function
:param qq: two-element tuple (or list) with quantiles that define lower and upper limits for image intensity normalization. Default is `[0.5, 97.5]`. Will be superseded by manually provided `vmin` and `vmax` arguments.
:param color: Color to use for plotting the positions of data points, optional
:param show_colorbar: Whether to show a colorbar alongside the image
:param show_axis: Whether to show the axes around the image
:param show_dots: Whether to overlay the positions of data points onto the plot
......@@ -124,7 +125,7 @@ def binned_map(x, y, value, bins=16, statistic='mean', qq=[0.5, 97.5], show_colo
if show_dots:
ax.set_autoscale_on(False)
ax.plot(x, y, 'b.', alpha=0.3)
ax.plot(x, y, '.', color=color, alpha=0.3)
def plot_cutout(cutout, planes=['image', 'template', 'diff', 'mask'], fig=None, mark_x=None, mark_y=None, mark_r=5.0, title=None, additional_title=None, **kwargs):
"""Routine for displaying various image planes from the cutout structure returned by :func:`stdpipe.cutouts.get_cutout`.
......
......@@ -89,8 +89,17 @@ def download(url, filename=None, overwrite=False, verbose=False):
def get_obs_time(header=None, filename=None, string=None, get_datetime=False, verbose=False):
"""
Extract date and time of observations from FITS headers of common formats.
Returns astropy Time object.
Extract date and time of observations from FITS headers of common formats, or from a string.
Will try various FITS keywords that may contain the time information - `DATE_OBS`, `DATE`, `TIME_OBS`, `UT`, 'MJD', 'JD'.
:param header: FITS header containing the information on time of observations
:param filename: If `header` is not set, the FITS header will be loaded from the file with this name
:param string: If provided, the time will be parsed from the string instead of FITS header
:param get_datetime: Whether to return the time as a standard Python :class:`datetime.datetime` object instead of Astropy Time
:param verbose: Whether to show verbose messages during the run of the function or not. May be either boolean, or a `print`-like function.
:returns: :class:`astropy.time.Time` object corresponding to the time of observations, or a :class:`datetime.datetime` object if :code:`get_datetime=True`
"""
# Simple wrapper around print for logging in verbose mode only
......@@ -98,7 +107,22 @@ def get_obs_time(header=None, filename=None, string=None, get_datetime=False, ve
# Simple wrapper to display parsed value and convert it as necessary
def convert_time(time):
time = Time(time)
if isinstance(time, float):
# Try to parse floating-point value as MJD or JD, depending on the value
if time > 0 and time < 100000:
log('Assuming it is MJD')
time = Time(time, format='mjd')
elif time > 2400000 and time < 2500000:
log('Assuming it is JD')
time = Time(time, format='jd')
else:
# Then it is probably an Unix time?..
log('Assuming it is Unix time')
time = Time(time, format='unix')
else:
time = Time(time)
log('Time parsed as:', time)
if get_datetime:
return time.datetime
......@@ -117,7 +141,7 @@ def get_obs_time(header=None, filename=None, string=None, get_datetime=False, ve
log('Loading FITS header from', filename)
header = fits.getheader(filename)
for dkey in ['DATE-OBS']:
for dkey in ['DATE-OBS', 'DATE', 'TIME-OBS', 'UT', 'MJD', 'JD']:
if dkey in header:
log('Found ' + dkey + ':', header[dkey])
# First try to parse standard ISO time
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment