Commit f187e0dc authored by karpov-sv's avatar karpov-sv
Browse files

Merge branch 'master' of https://github.com/karpov-sv/stdpipe into docker

parents a97d2195 dbc13dfa
......@@ -40,6 +40,7 @@ extensions = [
"sphinx.ext.intersphinx",
"sphinx.ext.autosectionlabel",
"sphinx.ext.viewcode",
"sphinx-prompt",
]
mathjax_config = {
......
......@@ -22,6 +22,10 @@ Both routines return the results as a standard Astropy Table, ordered by the obj
For :func:`~stdpipe.photometry.get_objects_sextractor`, the output may also contain columns related to PSF photometry (`x_psf`, `y_psf`, `flux_psf`, `fluxerr_psf`, `mag_psf`, `magerr_psf`, `chi2_psf`, `spread_model`, `spreaderr_model`), as well as any additional measuremet parameter requested through `extra_params` argument.
.. attention::
*STDPipe* (as well as SEP library) uses pixel coordinate convention with `(0, 0)` as the origin - that differs from *SExtractor* that uses `(1, 1)` as origin of coordinates! So the routine transparently converts `x` and `y` (as well as `x_psf` and `y_psf`) in the output from *SExtractor* to proper origin, so that the coordinates of objects detected by both routines are in the same system. However, if you manually add some pixel coordinate parameter to the output through `extra_params` argument, they will not be appropriately adjusted!
Below are some examples of object detection.
.. code-block:: python
......@@ -72,9 +76,9 @@ Finally, using SExtractor star/galaxy separators - `CLASS_STAR` and `SPREAD_MODE
print('SPREAD_MODEL = %.3f +/- %.3f, CLASS_STAR = %.2f' %
(cand['spread_model'], cand['spreaderr_model'], cand['CLASS_STAR']))
.. attention::
The most important problem with object detection using these routines is handling of blended objects, as the codes we are using can't properly deblend close groups, except for simplest cases.
.. note::
The most important problem with object detection using these routines is handling of blended objects, as the codes we are using can't properly deblend close groups, except for simplest cases.
.. autofunction:: stdpipe.photometry.get_objects_sextractor
:noindex:
......
=============================================
STDPipe - Simple Transient Detection Pipeline
---------------------------------------------
=============================================
*AKA: random codes noone else will ever use*
......@@ -26,10 +27,11 @@ User Guide
.. toctree::
:maxdepth: 3
About <self>
installation
usage
todo
contributing
todo
API documentation <api/modules>
Contributing
......
# Installation
Installation
============
*STDpipe* is available at [https://github.com/karpov-sv/stdpipe](https://github.com/karpov-sv/stdpipe) and is mirrored at [https://gitlab.in2p3.fr/icare/stdpipe](https://gitlab.in2p3.fr/icare/stdpipe)
*STDpipe* is available at https://github.com/karpov-sv/stdpipe and is mirrored at https://gitlab.in2p3.fr/icare/stdpipe.
In order to use it, you will need a working Python (>=3.6) installation, as well as a number of additional Python libraries and external packages. Below you will find some basic instructions how to set it up (using Anaconda environment as an example, but you may use any other as well) if you do not have it already.
## Preparing the environment (optional)
Preparing the environment (optional)
------------------------------------
The steps highlighted below are primarily for Linux and MacOS systems.
Windows users are advised to use WSL (preferrably Ubuntu 20.04) for smooth installation.
......@@ -14,127 +14,145 @@ Ubuntu 20.04 is available for free download on Microsoft Store.
You may safely skip these steps if you already have a working Python environment where you would like to install *STDPipe*.
Installing Anaconda
^^^^^^^^^^^^^^^^^^^
On your Linux/MacOS/WSL terminal, run the following commands to install `Anaconda <https://www.anaconda.com>`_ (replace 5.3.1 by the latest version, and adjust operating system name):
**Installing Anaconda**
.. prompt:: bash
On your Linux/MacOS/WSL terminal, run the following commands to install [Anaconda](https://www.anaconda.com) (replace 5.3.1 by the latest version, and adjust operating system name):
wget https://repo.anaconda.com/archive/Anaconda3-5.3.1-Linux-x86_64.sh
* $ wget https://repo.anaconda.com/archive/Anaconda3-5.3.1-Linux-x86_64.sh
.. prompt:: bash
* $ bash Anaconda3-5.3.1-Linux-x86_64.sh
bash Anaconda3-5.3.1-Linux-x86_64.sh
(For 32-bit installation, skip the ‘_64’ in both commands)
NOTE: If you already have Anaconda3 installed, please make sure that it is updated to the latest version (conda update --all). Also check that you do not have multiple
NOTE: If you already have Anaconda3 installed, please make sure that it is updated to the latest version (:code:`conda update --all`). Also check that you do not have multiple
versions of python installed in usr/lib/ directory as it can cause version conflicts while installing dependencies.
Now do:
* $ conda update --all
**Creating separate environment**
.. prompt:: bash
Create a new environment using this command (environment name is `stdpipe` in this case):
conda update --all
* $ conda create --name stdpipe
* $ conda activate stdpipe
NOTE: If this gives an error like:
CommandNotFoundError:
Your shell has not been properly configured to use 'conda activate', then run:
Installing basic dependencies
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
* $ source ~/anaconda3/etc/profile.d/conda.sh
.. prompt:: bash
then proceed with conda activate nmma_env.
conda install numpy scipy astropy matplotlib pandas
Check python and pip version like this:
Install conda-forge dependencies
* $ python --version
* $ pip --version
.. prompt:: bash
Python 3.7 and above and Pip 21.2 and above is ideal for this installation. It is recommended to update these for your installation.
conda install -c conda-forge astroscrappy
Install conda astropy dependencies
**Installing basic dependencies**
.. prompt:: bash
* $ conda install numpy scipy astropy matplotlib pandas
conda install -c astropy astroquery
Install conda-forge dependencies
* $ conda install -c conda-forge esutil astroscrappy
STDPipe installation
--------------------
Install conda astropy dependencies
Clone the STDPipe repository from GitHub at https://github.com/karpov-sv/stdpipe
* $ conda install -c astropy astroquery
.. prompt:: bash
git clone https://github.com/karpov-sv/stdpipe.git
## STDPipe installation
Change directory to the stdpipe folder:
.. prompt:: bash
Clone the STDPipe repository from GitHub at [https://github.com/karpov-sv/stdpipe](https://github.com/karpov-sv/stdpipe)
cd stdpipe
* $ git clone https://github.com/karpov-sv/stdpipe.git
Use the command below to install the rest of dependencies and the package itself in an *editable* manner so that it will be updated automatically when you update the code:
Change directory to the stdpipe folder:
.. prompt:: bash
* $ cd stdpipe
python setup.py develop
Use the commands below to install the rest of dependencies and the package itself in an *editable* manner so that it will be updated automatically when you update the code:
.. note::
* $ python setup.py develop
Alternative installation command (try it if the one above fails - they use slightly different strategies of installing the dependencies, so results may really vary!) would be
Alternative installation command (try it if the one above fails - they use slightly different strategies of installing the dependencies, so results may really vary!) would be
.. prompt:: bash
* $ pip install -e .
pip install -e .
**Keeping up to date**
Keeping up to date
^^^^^^^^^^^^^^^^^^
The command above installs the package to your Python environment in an *editable* way - it means that all changes you may make to the source tree (where you cloned the code) will immediately be reflected in the installed package, you do not need to repeat the installation.
As the code base in the repository evolves fast -- new features are being added, bugs fixed, etc -- it is a good idea to update your cloned code from the upstream often. The following command from inside stdpipe folder will do it:
* $ git pull
.. prompt:: bash
git pull
**Quick testing the installation**
Quick testing the installation
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Run the following commands:
Run the following commands inside your python (e.g. after typing :code:`ipython`):
* $ ipython
* import stdpipe
* import stdpipe.photometry
* import stdpipe.cutouts
.. prompt:: python
import stdpipe
import stdpipe.photometry
import stdpipe.cutouts
NOTE (Okay, last one!): if everything is ok, it's the end of the installation. But in case it shows that such-and-such modules are absent, feel free to install those modules by visiting their anaconda documentation and install
those with their given commands. In case modules like photutils and statsmodels are needed, don't hesitate to do it with pip (normally it shouldn't happen), but some modules may not install correctly in case of disturbance.
This instruction file will likely cover the issues you might face during your installation. However, please open issues on GitHub if there appear to be unresolvable conflicts.
This instruction page will likely cover the issues you might face during your installation. However, please open issues on GitHub if there appear to be unresolvable conflicts.
## Installation of external packages
Installation of external packages
---------------------------------
*STDPipe* makes use of a number of (optional) external packages:
- [SExtractor](https://github.com/astromatic/sextractor)
- [SCAMP](https://github.com/astromatic/scamp)
- [PSFEx](https://github.com/astromatic/psfex)
- [SWarp](https://github.com/astromatic/swarp)
- [HOTPANTS](https://github.com/acbecker/hotpants)
- [Astrometry.Net](https://github.com/dstndstn/astrometry.net)
- `SExtractor <https://github.com/astromatic/sextractor>`__
- `SCAMP <https://github.com/astromatic/scamp>`__
- `PSFEx <https://github.com/astromatic/psfex>`__
- `SWarp <https://github.com/astromatic/swarp>`__
- `HOTPANTS <https://github.com/acbecker/hotpants>`__
- `Astrometry.Net <https://github.com/dstndstn/astrometry.net>`__
Most of them are also available in the repositories of various Linux distributions, and may be conveniently installed from there (see below).
HOTPANTS image subtraction package cannot presently (as far as I know) be installed from any package manager, and has to be compiled manually.
### Ubuntu
.. attention::
If HOTPANTS compilation fails for you on the linking stage with a number of :code:`multiple definition of` error messages - that's a `known bug <https://github.com/acbecker/hotpants/issues/5>`__ related to some recent changes in GCC compiler defaults. You may easily fix it by editing the :file:`Makefile` and adding :code:`-fcommon` switch among the others in the `COPTS` options (line `30 <https://github.com/acbecker/hotpants/blob/master/Makefile#L30>`__ at the moment of writing).
* $ sudo apt install sextractor scamp psfex swarp
Ubuntu
^^^^^^
.. prompt:: bash
sudo apt install sextractor scamp psfex swarp
Astrometry.Net may also be installed from repository, but might require additional manual configuration steps (and quite a lot of disk space for larger indices!), so install it only when you really need it, and when you really know what you are doing!
* $ sudo apt install astrometry.net
.. prompt:: bash
sudo apt install astrometry.net
Anaconda
^^^^^^^^
### Anaconda
.. prompt:: bash
* $ conda install -c conda-forge astromatic-source-extractor astromatic-scamp astromatic-psfex astromatic-swarp
conda install -c conda-forge astromatic-source-extractor astromatic-scamp astromatic-psfex astromatic-swarp
Point Spread Function (PSF) models
==================================
*STDPipe* includes basic support for point spread function (PSF) construction and analysis through the interface to `PSFEx <https://github.com/astromatic/psfex>`_ code that is able to build supersampled PSF models from the object lists created using `SExtractor <https://github.com/astromatic/sextractor>`_.
Please consider checking `its documentation <https://psfex.readthedocs.io>`__ to better understand the concepts of PSFEx operation and its possible configuration options.
The :func:`stdpipe.psf.run_psfex` routine transparently calls the `SExtractor`_ on the image and then uses its output (it needs a special set of measured object parameters and postage stamps for them) to run `PSFEx`_. The results are then parsed and returned as a structure representing the PSF model file as described `here <https://psfex.readthedocs.io/en/latest/Appendices.html#psf-file-format-description>`__. This structure may later be used to visualize PSF shape at arbitrary position inside the image, as well as to generate artificial stars and inject them into the image.
.. autofunction:: stdpipe.psf.run_psfex
:noindex:
.. attention::
While `PSFEx`_ may build PSF models depending on a variety of object parameters as described `here <https://psfex.readthedocs.io/en/latest/Working.html#managing-psf-variations>`__, *STDPipe* currently allows reading and interpreting only the ones depending on the position on the image (e.g. representing optical distortions). On the other hand, the original PSFEx model file may be stored to the file using `psffile` parameter of :func:`stdpipe.psf.run_psfex`, and then supplied directly to `SExtractor`_ through `psf` parameter of :func:`stdpipe.photometry.get_objects_sextractor` to be used for SExtractor PSF photometry.
Using PSF models
----------------
The PSF model built by `PSFEx`_ may also be loaded from PSFEx model file. This model may then be used to construct PSF stamps both in supersampled PSF model pixels, or original image pixels.
.. autofunction:: stdpipe.psf.load_psf
:noindex:
.. autofunction:: stdpipe.psf.get_supersampled_psf_stamp
:noindex:
.. autofunction:: stdpipe.psf.get_psf_stamp
:noindex:
Placing artificial stars into the image
---------------------------------------
*STDPipe* also contains a couple of convenience functions to directly place a PSF model at a given image position - that would correspond to artificial point source injection into the image. Higher-level one, :func:`stdpipe.pipeline.place_random_stars`, inserts a number of randomly generated sources with realistically (log-uniform in fluxes) distributed brightness at random positions, as well as returns the catalogue of injected stars that may be used e.g. in deriving the object detection efficiency.
.. code-block:: python
# Derive PSF model assuming that it does not change over the image
psf_model = psf.run_psfex(image, mask=mask, order=0, verbose=True)
# ..and now perform the detection efficiency analysis
sims = [] # will hold simulated stars
# We will repeatedly inject the stars (20 at a time), detect objects, and
# match them in order to see whether simulated stars are detectable
for _ in tqdm(range(100)):
# We will operate on a copy of original image
image1 = image.copy()
# We will use this value as a saturation level
saturation = 50000
# Simulate 20 random stars
sim = pipeline.place_random_stars(image1, psf_model, nstars=20, minflux=10, maxflux=1e6,
wcs=wcs, gain=gain, saturation=saturation)
# Initial metadata for injected stars
sim['mag_calib'] = sim['mag'] + zero_fn(sim['x'], sim['y']) # Apply zero point from photometric calibration
sim['detected'] = False
sim['mag_measured'] = np.nan
sim['magerr_measured'] = np.nan
sim['flags_measured'] = np.nan
# Mask corresponding to the saturation level
mask1 = image1 >= saturation
obj1 = photometry.get_objects_sextractor(image1, mask=mask|mask1, r0=1, aper=5.0,
wcs=wcs, gain=gain, minarea=3, sn=5)
# Apply zero point from photometric calibration
obj1['mag_calib'] = obj1['mag'] + zero_fn(obj1['x'], obj1['y'])
# Positional match within FWHM/2 radius
oidx,sidx,dist = astrometry.spherical_match(obj1['ra'], obj1['dec'], sim['ra'], sim['dec'],
pixscale*np.median(obj1['fwhm'])/2)
# Mark matched stars
sim['detected'][sidx] = True
# Also store measured magnitude, its error and flags
sim['mag_measured'][sidx] = obj1['mag_calib'][oidx]
sim['magerr_measured'][sidx] = obj1['magerr'][oidx]
sim['flags_measured'][sidx] = obj1['flags'][oidx]
sims.append(sim)
# Now we may stack the data from all runs into a single table
from astropy.table import vstack
sims = vstack(sims)
# ..and plot it as a histograms of magnitudes
h0,b0,_ = plt.hist(sims['mag_calib'], range=[12,22], bins=50, alpha=0.3, label='All simulated stars');
h1,b1,_ = plt.hist(sims['mag_calib'][sims['detected']], range=[12,22], bins=50, alpha=0.3, label='Detected');
h2,b2,_ = plt.hist(sims['mag_calib'][idx], range=[12,22], bins=50, alpha=0.3, label='Detected and unflagged');
plt.legend()
plt.xlabel('Injected magnitude')
plt.show()
# ..or as a detection efficiency
plt.plot(0.5*(b0[1:]+b0[:-1]), h1/h0, 'o-', label='Detected')
plt.plot(0.5*(b0[1:]+b0[:-1]), h2/h0, 'o-', label='Detected and unflagged')
plt.legend()
plt.xlabel('Injected magnitude')
plt.title('Fraction of detected artificial stars')
The complete example of detection efficiency analysis may be seen in the `corresponding notebook <https://github.com/karpov-sv/stdpipe/blob/master/notebooks/simulated_stars.ipynb>`_.
.. attention::
Note that the flux that you specify for the artificial star is a total flux - it will not necessarily correspond to the one measured inside an aperture (especially if the aperture size is quite small). The difference (*aperture correction*) is due to the fraction of the total flux that falls outside of your aperture, and it has to be somehow estimated and corrected if you wish to compare the generated fluxes of injected stars with measured values.
.. autofunction:: stdpipe.psf.place_psf_stamp
:noindex:
.. autofunction:: stdpipe.pipeline.place_random_stars
:noindex:
PSF photometry in SExtractor
----------------------------
Derived PSF models may also be used for performing PSF photometry in SExtractor.
As currently *STDPipe* lacks the routine for saving PSF models back to FITS files, in order to use PSF photometry you will need to pass `psffile` argument to :func:`stdpipe.psf.run_psfex` in order to directly store the produced PSF model into the file. Then, this file may be used in SExtractor by passing its path as a `psf` argument into :func:`stdpipe.photometry.get_objects_sextractor`. It will then produce the list of detected objects having a set of additional measured parameters including `flux_psf`, `fluxerr_psf`, `mag_psf`, `magerr_psf` - they correspond to the flux derived from fitting the objects with the PSF model.
.. code-block:: python
# We probably do not have enough stars to study PSF spatial variance, so we use order=0 here
psf_model = psf.run_psfex(image, mask=mask, order=0, gain=gain, psffile='/tmp/psf.psf', verbose=True)
obj_psf = photometry.get_objects_sextractor(image, mask=mask, aper=3.0, edge=10, gain=gain,
wcs=wcs, psf='/tmp/psf.psf', verbose=False)
# Now we may calibrate the photometry using PSF fluxes
m_psf = pipeline.calibrate_photometry(obj_psf, cat, sr=1/3600,
obj_col_mag='mag_psf', obj_col_mag_err='magerr_psf',
cat_col_mag='rmag', cat_col_mag1='gmag', cat_col_mag2='rmag',
order=0, verbose=True)
......@@ -4,6 +4,7 @@ sphinx_math_dollar
numpydoc
myst_parser
sphinxcontrib-apidoc
sphinx-prompt
extension-helpers
numpy
......@@ -15,7 +16,6 @@ astroquery>0.4.1
sep
photutils
statsmodels
esutil==0.6.0
tqdm
psycopg2
regions
......@@ -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
---------------
......@@ -48,6 +49,9 @@ Data processing
subtraction
transients
cutouts
psf
utils
Common principles
-----------------
......@@ -64,6 +68,9 @@ Design principles:
- everything operates on temporary files, nothing is kept after the run unless explicitly asked for
- temporary files are created in unique temporary directories for each run, so several instances of routines may be safely run in parallel
- All image (pixel) coordinates have origin at `(0, 0)` - the ones returned from and passed to e.g. *SExtractor*, *SCAMP* and likes (that are based on `(1, 1)` origin) are transparently converted
Common conventions for routine arguments:
- Most of functions accept `verbose` argument that controls the amount of informational outputs the function produces. You may use `verbose=True` to see the details of what exactly the function is doing. Also, you may pass any `print`-like function to this argument to receive the messages instead of printing - so e.g. they may be logged.
......@@ -92,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:
This diff is collapsed.
......@@ -16,9 +16,10 @@ requirements = [
'sep',
'photutils',
'statsmodels',
'esutil',
'tqdm',
'regions'
'regions',
'python-dateutil',
'requests',
]
setup(
......
......@@ -3,13 +3,12 @@ from __future__ import absolute_import, division, print_function, unicode_litera
import os, tempfile, shutil, shlex, re, warnings
import numpy as np
from esutil import coords, htm
from astropy.wcs import WCS
from astropy.io import fits
from astropy.wcs.utils import fit_wcs_from_points
from astropy.coordinates import SkyCoord
from astropy.coordinates import SkyCoord, search_around_sky
from astropy.table import Table
from astropy import units as u
from scipy.stats import chi2
......@@ -34,10 +33,10 @@ def get_frame_center(filename=None, header=None, wcs=None, width=None, height=No
elif shape is not None:
height,width = shape
[ra1],[dec1] = wcs.all_pix2world([0], [0], 1)
[ra0],[dec0] = wcs.all_pix2world([width/2], [height/2], 1)
ra1,dec1 = wcs.all_pix2world(0, 0, 0)
ra0,dec0 = wcs.all_pix2world(width/2, height/2, 0)
sr = coords.sphdist(ra0, dec0, ra1, dec1)[0]
sr = spherical_distance(ra0, dec0, ra1, dec1)
return ra0, dec0, sr
......@@ -71,6 +70,47 @@ def xyztoradec(xyz):
return (np.rad2deg(ra), np.rad2deg(dec))
def spherical_distance(ra1, dec1, ra2, dec2):
"""Spherical distance.
:param ra1: First point or set of points RA
:param dec1: First point or set of points Dec
:param ra2: Second point or set of points RA
:param dec2: Second point or set of points Dec
:returns: Spherical distance in degrees
"""
x = np.sin(np.deg2rad((ra1 - ra2)/2))
x *= x;
y = np.sin(np.deg2rad((dec1 - dec2)/2))
y *= y;
z = np.cos(np.deg2rad((dec1 + dec2)/2))
z *= z;
return np.rad2deg(2*np.arcsin(np.sqrt(x*(z - y) + y)))
def spherical_match(ra1, dec1, ra2, dec2, sr=1/3600):
"""Positional match on the sphere for two lists of coordinates.
Aimed to be a direct replacement for :func:`esutil.htm.HTM.match` method with :code:`maxmatch=0`.
:param ra1: First set of points RA
:param dec1: First set of points Dec
:param ra2: Second set of points RA
:param dec2: Second set of points Dec
:param sr: Maximal acceptable pair distance to be considered a match, in degrees
:returns: Two parallel sets of indices corresponding to matches from first and second lists, along with the pairwise distances in degrees
"""
idx1,idx2,dist,_ = search_around_sky(SkyCoord(ra1, dec1, unit='deg'), SkyCoord(ra2, dec2, unit='deg'), sr*u.deg)
dist = dist.deg # convert to degrees
return idx1, idx2, dist
def get_objects_center(obj, col_ra='ra', col_dec='dec'):
"""
Returns the center RA, Dec, and radius in degrees for a cloud of objects on the sky.
......@@ -79,7 +119,7 @@ def get_objects_center(obj, col_ra='ra', col_dec='dec'):
xyz0 = np.mean(xyz, axis=1)
ra0,dec0 = xyztoradec(xyz0)
sr0 = np.max(coords.sphdist(ra0, dec0, obj[col_ra], obj[col_dec]))
sr0 = np.max(spherical_distance(ra0, dec0, obj[col_ra