SourceXtractorPlusPlus  0.11
Please provide a description of the project.
model_fitting.py
Go to the documentation of this file.
1 # -*- coding: utf-8 -*-
2 
3 # Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université
4 #
5 # This library is free software; you can redistribute it and/or modify it under
6 # the terms of the GNU Lesser General Public License as published by the Free
7 # Software Foundation; either version 3.0 of the License, or (at your option)
8 # any later version.
9 #
10 # This library is distributed in the hope that it will be useful, but WITHOUT
11 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 # FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
13 # details.
14 #
15 # You should have received a copy of the GNU Lesser General Public License
16 # along with this library; if not, write to the Free Software Foundation, Inc.,
17 # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 from __future__ import division, print_function
19 
20 import sys
21 from enum import Enum
22 
23 import _SourceXtractorPy as cpp
24 from .measurement_images import MeasurementGroup
25 
26 from astropy import units as u
27 from astropy.coordinates import SkyCoord
28 from astropy.coordinates import Angle
29 
30 import math
31 
32 
33 class RangeType(Enum):
34  LINEAR = 1
35  EXPONENTIAL = 2
36 
37 
38 class Range(object):
39  r"""
40  Limit, and normalize, the range of values for a model fitting parameter.
41 
42 
43  Parameters
44  ----------
45  limits : a tuple (min, max), or a callable that receives a source, and returns a tuple (min, max)
46  type : RangeType
47 
48  Notes
49  -----
50  RangeType.LINEAR
51  Normalized to engine space using a sigmoid function
52 
53  .. math::
54 
55  engine = \ln \frac{world - min}{max-world} \\
56  world = min + \frac{max - min}{1 + e^{engine}}
57 
58  RangeType.EXPONENTIAL
59  Normalized to engine space using an exponential sigmoid function
60 
61  .. math::
62 
63  engine = \ln \left( \frac{\ln(world/min)}{\ln(max /world)} \right) \\
64  world = min * e^\frac{ \ln(max / min) }{ (1 + e^{-engine}) }
65 
66  """
67 
68  def __init__(self, limits, type):
69  """
70  Constructor.
71  """
72  self.__limits = limits
73  self.__type = type
74 
75  def get_limits(self):
76  """
77  Returns
78  -------
79  callable
80  Receives a source, and returns a tuple (min, max)
81  """
82  return self.__limits if hasattr(self.__limits, '__call__') else lambda v,o: self.__limits
83 
84  def get_type(self):
85  """
86  Returns
87  -------
88  RangeType
89  """
90  return self.__type
91 
92  def __str__(self):
93  """
94  Returns
95  -------
96  str
97  Human readable representation for the object
98  """
99  res = '['
100  if hasattr(self.__limits, '__call__'):
101  res += 'func'
102  else:
103  res += '{},{}'.format(self.__limits[0], self.__limits[1])
104  type_str = {
105  RangeType.LINEAR : 'LIN',
106  RangeType.EXPONENTIAL : 'EXP'
107  }
108  res += ',{}]'.format(type_str[self.__type])
109  return res
110 
111 
112 class Unbounded(object):
113  """
114  Unbounded, but normalize, value of a model fitting parameter
115 
116  Parameters
117  ----------
118  normalization_factor: float, or a callable that receives the initial value parameter value and a source,
119  and returns a float
120  The world value which will be normalized to 1 in engine coordinates
121  """
122 
123  def __init__(self, normalization_factor=1):
124  """
125  Constructor.
126  """
127  self.__normalization_factor = normalization_factor
128 
130  """
131  Returns
132  -------
133  callable
134  Receives the initial parameter value and a source, and returns the world value which will be
135  normalized to 1 in engine coordinates
136  """
137  return self.__normalization_factor if hasattr(self.__normalization_factor, '__call__') else lambda v,o: self.__normalization_factor
138 
139  def __str__(self):
140  """
141  Returns
142  -------
143  str
144  Human readable representation for the object
145  """
146  res = '['
147  if hasattr(self.__normalization_factor, '__call__'):
148  res += 'func'
149  else:
150  res += '{}'.format(self.__normalization_factor)
151  res += ']'
152  return res
153 
154 
155 constant_parameter_dict = {}
156 free_parameter_dict = {}
157 dependent_parameter_dict = {}
158 
159 
160 def print_parameters(file=sys.stderr):
161  """
162  Print a human-readable representation of the configured model fitting parameters.
163 
164  Parameters
165  ----------
166  file : file object
167  Where to print the representation. Defaults to sys.stderr
168  """
169  if constant_parameter_dict:
170  print('Constant parameters:', file=file)
171  for n in constant_parameter_dict:
172  print(' {}: {}'.format(n, constant_parameter_dict[n]), file=file)
173  if free_parameter_dict:
174  print('Free parameters:', file=file)
175  for n in free_parameter_dict:
176  print(' {}: {}'.format(n, free_parameter_dict[n]), file=file)
177  if dependent_parameter_dict:
178  print('Dependent parameters:', file=file)
179  for n in dependent_parameter_dict:
180  print(' {}: {}'.format(n, dependent_parameter_dict[n]), file=file)
181 
182 
183 class ParameterBase(cpp.Id):
184  """
185  Base class for all model fitting parameter types.
186  Can not be used directly.
187  """
188 
189  def __str__(self):
190  """
191  Returns
192  -------
193  str
194  Human readable representation for the object
195  """
196  return '(ID:{})'.format(self.id)
197 
198 
200  """
201  A parameter with a single value that remains constant. It will not be fitted.
202 
203  Parameters
204  ----------
205  value : float, or callable that receives a source and returns a float
206  Value for this parameter
207  """
208 
209  def __init__(self, value):
210  """
211  Constructor.
212  """
213  ParameterBase.__init__(self)
214  self.__value = value
215  constant_parameter_dict[self.id] = self
216 
217  def get_value(self):
218  """
219  Returns
220  -------
221  callable
222  Receives a source and returns a value for the parameter
223  """
224  return self.__value if hasattr(self.__value, '__call__') else lambda o: self.__value
225 
226  def __str__(self):
227  """
228  Returns
229  -------
230  str
231  Human readable representation for the object
232  """
233  res = ParameterBase.__str__(self)[:-1] + ', value:'
234  if hasattr(self.__value, '__call__'):
235  res += 'func'
236  else:
237  res += str(self.__value)
238  return res + ')'
239 
240 
242  """
243  A parameter that will be fitted by the model fitting engine.
244 
245  Parameters
246  ----------
247  init_value : float or callable that receives a source, and returns a float
248  Initial value for the parameter.
249  range : instance of Range or Unbounded
250  Defines if this parameter is unbounded or bounded, and how.
251 
252  See Also
253  --------
254  Unbounded
255  Range
256 
257  Examples
258  --------
259  >>> sersic = FreeParameter(2.0, Range((1.0, 7.0), RangeType.LINEAR))
260  """
261 
262  def __init__(self, init_value, range=Unbounded()):
263  """
264  Constructor.
265  """
266  ParameterBase.__init__(self)
267  self.__init_value = init_value
268  self.__range = range
269  free_parameter_dict[self.id] = self
270 
271  def get_init_value(self):
272  """
273  Returns
274  -------
275  callable
276  Receives a source, and returns an initial value for the parameter.
277  """
278  return self.__init_value if hasattr(self.__init_value, '__call__') else lambda o: self.__init_value
279 
280  def get_range(self):
281  """
282  Returns
283  -------
284  Unbounded or Range
285  """
286  return self.__range
287 
288  def __str__(self):
289  """
290  Returns
291  -------
292  str
293  Human readable representation for the object
294  """
295  res = ParameterBase.__str__(self)[:-1] + ', init:'
296  if hasattr(self.__init_value, '__call__'):
297  res += 'func'
298  else:
299  res += str(self.__init_value)
300  if self.__range:
301  res += ', range:' + str(self.__range)
302  return res + ')'
303 
304 
306  """
307  A DependentParameter is not fitted by itself, but its value is derived from another Parameters, whatever their type:
308  FreeParameter, ConstantParameter, or other DependentParameter
309 
310  Parameters
311  ----------
312  func : callable
313  A callable that will be called with all the parameters specified in this constructor each time a new
314  evaluation is needed.
315  params : list of ParameterBase
316  List of parameters on which this DependentParameter depends.
317 
318  Examples
319  --------
320  >>> flux = get_flux_parameter()
321  >>> mag = DependentParameter(lambda f: -2.5 * np.log10(f) + args.mag_zeropoint, flux)
322  >>> add_output_column('mf_mag_' + band, mag)
323  """
324 
325  def __init__(self, func, *params):
326  """
327  Constructor.
328  """
329  ParameterBase.__init__(self)
330  self.func = func
331  self.params = [p.id for p in params]
332  dependent_parameter_dict[self.id] = self
333 
334 
336  """
337  Convenience function for the position parameter X and Y.
338 
339  Returns
340  -------
341  x : FreeParameter
342  X coordinate, starting at the X coordinate of the centroid and linearly limited to X +/- the object radius.
343  y : FreeParameter
344  Y coordinate, starting at the Y coordinate of the centroid and linearly limited to Y +/- the object radius.
345  Notes
346  -----
347  X and Y are fitted on the detection image X and Y coordinates. Internally, these are translated to measurement
348  images using the WCS headers.
349  """
350  return (
351  FreeParameter(lambda o: o.get_centroid_x(), Range(lambda v,o: (v-o.get_radius(), v+o.get_radius()), RangeType.LINEAR)),
352  FreeParameter(lambda o: o.get_centroid_y(), Range(lambda v,o: (v-o.get_radius(), v+o.get_radius()), RangeType.LINEAR))
353  )
354 
355 
357  """
358  Possible flux types to use as initial value for the flux parameter.
359  Right now, only isophotal is supported.
360  """
361  ISO = 1
362 
363 
364 def get_flux_parameter(type=FluxParameterType.ISO, scale=1):
365  """
366  Convenience function for the flux parameter.
367 
368  Parameters
369  ----------
370  type : int
371  One of the values defined in FluxParameterType
372  scale : float
373  Scaling of the initial flux. Defaults to 1.
374 
375  Returns
376  -------
377  flux : FreeParameter
378  Flux parameter, starting at the flux defined by `type`, and limited to +/- 1e3 times the initial value.
379  """
380  func_map = {
381  FluxParameterType.ISO : 'get_iso_flux'
382  }
383  return FreeParameter(lambda o: getattr(o, func_map[type])() * scale, Range(lambda v,o: (v * 1E-3, v * 1E3), RangeType.EXPONENTIAL))
384 
385 
386 prior_dict = {}
387 
388 
389 class Prior(cpp.Id):
390  """
391  Model a Gaussian prior on a given parameter.
392 
393  Parameters
394  ----------
395  param : ParameterBase
396  Model fitting parameter
397  value : float or callable that receives a source and returns a float
398  Mean of the Gaussian
399  sigma : float or callable that receives a source and returns a float
400  Standard deviation of the Gaussian
401  """
402 
403  def __init__(self, param, value, sigma):
404  """
405  Constructor.
406  """
407  cpp.Id.__init__(self)
408  self.param = param.id
409  self.value = value if hasattr(value, '__call__') else lambda o: value
410  self.sigma = sigma if hasattr(sigma, '__call__') else lambda o: sigma
411 
412 
413 def add_prior(param, value, sigma):
414  """
415  Add a prior to the given parameter.
416 
417  Parameters
418  ----------
419  param : ParameterBase
420  value : float or callable that receives a source and returns a float
421  Mean of the Gaussian
422  sigma : float or callable that receives a source and returns a float
423  Standard deviation of the Gaussian
424  """
425  prior = Prior(param, value, sigma)
426  prior_dict[prior.id] = prior
427 
428 
429 frame_models_dict = {}
430 
431 
432 def _set_model_to_frames(group, model):
433  for x in group:
434  if isinstance(x, tuple):
435  _set_model_to_frames(x[1], model)
436  else:
437  if not x.id in frame_models_dict:
438  frame_models_dict[x.id] = []
439  frame_models_dict[x.id].append(model.id)
440 
441 
442 def add_model(group, model):
443  """
444  Add a model to be fitted to the given group.
445 
446  Parameters
447  ----------
448  group : MeasurementGroup
449  model : ModelBase
450  """
451  if not isinstance(group, MeasurementGroup):
452  raise TypeError('Models can only be added on MeasurementGroup, got {}'.format(type(group)))
453  if not hasattr(group, 'models'):
454  group.models = []
455  group.models.append(model)
456  _set_model_to_frames(group, model)
457 
458 
459 def print_model_fitting_info(group, show_params=False, prefix='', file=sys.stderr):
460  """
461  Print a human-readable representation of the configured models.
462 
463  Parameters
464  ----------
465  group : MeasurementGroup
466  Print the models for this group.
467  show_params : bool
468  If True, print also the parameters that belong to the model
469  prefix : str
470  Prefix each line with this string. Used internally for indentation.
471  file : file object
472  Where to print the representation. Defaults to sys.stderr
473  """
474  if hasattr(group, 'models') and group.models:
475  print('{}Models:'.format(prefix), file=file)
476  for m in group.models:
477  print('{} {}'.format(prefix, m.to_string(show_params)), file=file)
478  for x in group:
479  if isinstance(x, tuple):
480  print('{}{}:'.format(prefix, x[0]), file=file)
481  print_model_fitting_info(x[1], show_params, prefix + ' ', file=file)
482 
483 
484 constant_model_dict = {}
485 point_source_model_dict = {}
486 sersic_model_dict = {}
487 exponential_model_dict = {}
488 de_vaucouleurs_model_dict = {}
489 params_dict = {"max_iterations": 100, "modified_chi_squared_scale": 10, "engine": ""}
490 
491 
492 def set_max_iterations(iterations):
493  """
494  Parameters
495  ----------
496  iterations : int
497  Max number of iterations for the model fitting.
498  """
499  params_dict["max_iterations"] = iterations
500 
501 
503  """
504  Parameters
505  ----------
506  scale : float
507  Sets u0, as used by the modified chi squared residual comparator, a function that reduces the effect of large
508  deviations.
509  Refer to the SourceXtractor++ documentation for a better explanation of how residuals are computed and how
510  this value affects the model fitting.
511  """
512  params_dict["modified_chi_squared_scale"] = scale
513 
514 
515 def set_engine(engine):
516  """
517  Parameters
518  ----------
519  engine : str
520  Minimization engine for the model fitting : levmar or gsl
521  """
522  params_dict["engine"] = engine
523 
524 
525 class ModelBase(cpp.Id):
526  """
527  Base class for all models.
528  """
529  pass
530 
531 
532 class CoordinateModelBase(ModelBase):
533  """
534  Base class for positioned models with a flux. It can not be used directly.
535 
536  Parameters
537  ----------
538  x_coord : ParameterBase or float
539  X coordinate (in the detection image)
540  y_coord : ParameterBase or float
541  Y coordinate (in the detection image)
542  flux : ParameterBase or float
543  Total flux
544  """
545 
546  def __init__(self, x_coord, y_coord, flux):
547  """
548  Constructor.
549  """
550  ModelBase.__init__(self)
551  self.x_coord = x_coord if isinstance(x_coord, ParameterBase) else ConstantParameter(x_coord)
552  self.y_coord = y_coord if isinstance(y_coord, ParameterBase) else ConstantParameter(y_coord)
553  self.flux = flux if isinstance(flux, ParameterBase) else ConstantParameter(flux)
554 
555 
557  """
558  Models a source as a point, spread by the PSF.
559 
560  Parameters
561  ----------
562  x_coord : ParameterBase or float
563  X coordinate (in the detection image)
564  y_coord : ParameterBase or float
565  Y coordinate (in the detection image)
566  flux : ParameterBase or float
567  Total flux
568  """
569 
570  def __init__(self, x_coord, y_coord, flux):
571  """
572  Constructor.
573  """
574  CoordinateModelBase.__init__(self, x_coord, y_coord, flux)
575  global point_source_model_dict
576  point_source_model_dict[self.id] = self
577 
578  def to_string(self, show_params=False):
579  """
580  Return a human readable representation of the model.
581 
582  Parameters
583  ----------
584  show_params: bool
585  If True, include information about the parameters.
586 
587  Returns
588  -------
589  str
590  """
591  if show_params:
592  return 'PointSource[x_coord={}, y_coord={}, flux={}]'.format(
593  self.x_coord, self.y_coord, self.flux)
594  else:
595  return 'PointSource[x_coord={}, y_coord={}, flux={}]'.format(
596  self.x_coord.id, self.y_coord.id, self.flux.id)
597 
598 
600  """
601  A model that is constant through all the image.
602 
603  Parameters
604  ----------
605  value: ParameterBase or float
606  Value to add to the value of all pixels from the model.
607  """
608 
609  def __init__(self, value):
610  """
611  Constructor.
612  """
613  ModelBase.__init__(self)
614  self.value = value if isinstance(value, ParameterBase) else ConstantParameter(value)
615  global constant_model_dict
616  constant_model_dict[self.id] = self
617 
618  def to_string(self, show_params=False):
619  """
620  Return a human readable representation of the model.
621 
622  Parameters
623  ----------
624  show_params: bool
625  If True, include information about the parameters.
626 
627  Returns
628  -------
629  str
630  """
631  if show_params:
632  return 'ConstantModel[value={}]'.format(self.value)
633  else:
634  return 'ConstantModel[value={}]'.format(self.value.id)
635 
636 
638  """
639  Base class for the Sersic, Exponential and de Vaucouleurs models. It can not be used directly.
640 
641  Parameters
642  ----------
643  x_coord : ParameterBase or float
644  X coordinate (in the detection image)
645  y_coord : ParameterBase or float
646  Y coordinate (in the detection image)
647  flux : ParameterBase or float
648  Total flux
649  effective_radius : ParameterBase or float
650  Ellipse semi-major axis, in pixels on the detection image.
651  aspect_ratio : ParameterBase or float
652  Ellipse ratio.
653  angle : ParameterBase or float
654  Ellipse rotation, in radians
655  """
656 
657  def __init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle):
658  """
659  Constructor.
660  """
661  CoordinateModelBase.__init__(self, x_coord, y_coord, flux)
662  self.effective_radius = effective_radius if isinstance(effective_radius, ParameterBase) else ConstantParameter(effective_radius)
663  self.aspect_ratio = aspect_ratio if isinstance(aspect_ratio, ParameterBase) else ConstantParameter(aspect_ratio)
664  self.angle = angle if isinstance(angle, ParameterBase) else ConstantParameter(angle)
665 
666 
668  """
669  Model a source with a Sersic profile.
670 
671  Parameters
672  ----------
673  x_coord : ParameterBase or float
674  X coordinate (in the detection image)
675  y_coord : ParameterBase or float
676  Y coordinate (in the detection image)
677  flux : ParameterBase or float
678  Total flux
679  effective_radius : ParameterBase or float
680  Ellipse semi-major axis, in pixels on the detection image.
681  aspect_ratio : ParameterBase or float
682  Ellipse ratio.
683  angle : ParameterBase or float
684  Ellipse rotation, in radians
685  n : ParameterBase or float
686  Sersic index
687  """
688 
689  def __init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle, n):
690  """
691  Constructor.
692  """
693  SersicModelBase.__init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle)
694  self.n = n if isinstance(n, ParameterBase) else ConstantParameter(n)
695  global sersic_model_dict
696  sersic_model_dict[self.id] = self
697 
698  def to_string(self, show_params=False):
699  """
700  Return a human readable representation of the model.
701 
702  Parameters
703  ----------
704  show_params: bool
705  If True, include information about the parameters.
706 
707  Returns
708  -------
709  str
710  """
711  if show_params:
712  return 'Sersic[x_coord={}, y_coord={}, flux={}, effective_radius={}, aspect_ratio={}, angle={}, n={}]'.format(
713  self.x_coord, self.y_coord, self.flux, self.effective_radius, self.aspect_ratio, self.angle, self.n)
714  else:
715  return 'Sersic[x_coord={}, y_coord={}, flux={}, effective_radius={}, aspect_ratio={}, angle={}, n={}]'.format(
716  self.x_coord.id, self.y_coord.id, self.flux.id, self.effective_radius.id, self.aspect_ratio.id, self.angle.id, self.n.id)
717 
718 
720  """
721  Model a source with an exponential profile (Sersic model with an index of 1)
722 
723  Parameters
724  ----------
725  x_coord : ParameterBase or float
726  X coordinate (in the detection image)
727  y_coord : ParameterBase or float
728  Y coordinate (in the detection image)
729  flux : ParameterBase or float
730  Total flux
731  effective_radius : ParameterBase or float
732  Ellipse semi-major axis, in pixels on the detection image.
733  aspect_ratio : ParameterBase or float
734  Ellipse ratio.
735  angle : ParameterBase or float
736  Ellipse rotation, in radians
737  """
738 
739  def __init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle):
740  """
741  Constructor.
742  """
743  SersicModelBase.__init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle)
744  global exponential_model_dict
745  exponential_model_dict[self.id] = self
746 
747  def to_string(self, show_params=False):
748  """
749  Return a human readable representation of the model.
750 
751  Parameters
752  ----------
753  show_params: bool
754  If True, include information about the parameters.
755 
756  Returns
757  -------
758  str
759  """
760  if show_params:
761  return 'Exponential[x_coord={}, y_coord={}, flux={}, effective_radius={}, aspect_ratio={}, angle={}]'.format(
762  self.x_coord, self.y_coord, self.flux, self.effective_radius, self.aspect_ratio, self.angle)
763  else:
764  return 'Exponential[x_coord={}, y_coord={}, flux={}, effective_radius={}, aspect_ratio={}, angle={}]'.format(
765  self.x_coord.id, self.y_coord.id, self.flux.id, self.effective_radius.id, self.aspect_ratio.id, self.angle.id)
766 
767 
769  """
770  Model a source with a De Vaucouleurs profile (Sersic model with an index of 4)
771 
772  Parameters
773  ----------
774  x_coord : ParameterBase or float
775  X coordinate (in the detection image)
776  y_coord : ParameterBase or float
777  Y coordinate (in the detection image)
778  flux : ParameterBase or float
779  Total flux
780  effective_radius : ParameterBase or float
781  Ellipse semi-major axis, in pixels on the detection image.
782  aspect_ratio : ParameterBase or float
783  Ellipse ratio.
784  angle : ParameterBase or float
785  Ellipse rotation, in radians
786  """
787 
788  def __init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle):
789  """
790  Constructor.
791  """
792  SersicModelBase.__init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle)
793  global de_vaucouleurs_model_dict
794  de_vaucouleurs_model_dict[self.id] = self
795 
796  def to_string(self, show_params=False):
797  """
798  Return a human readable representation of the model.
799 
800  Parameters
801  ----------
802  show_params: bool
803  If True, include information about the parameters.
804 
805  Returns
806  -------
807  str
808  """
809  if show_params:
810  return 'DeVaucouleurs[x_coord={}, y_coord={}, flux={}, effective_radius={}, aspect_ratio={}, angle={}]'.format(
811  self.x_coord, self.y_coord, self.flux, self.effective_radius, self.aspect_ratio, self.angle)
812  else:
813  return 'DeVaucouleurs[x_coord={}, y_coord={}, flux={}, effective_radius={}, aspect_ratio={}, angle={}]'.format(
814  self.x_coord.id, self.y_coord.id, self.flux.id, self.effective_radius.id, self.aspect_ratio.id, self.angle.id)
815 
816 
818  """
819  Coordinates in right ascension and declination
820 
821  Parameters
822  ----------
823  ra : float
824  Right ascension
825  dec : float
826  Declination
827  """
828 
829  def __init__(self, ra, dec):
830  """
831  Constructor.
832  """
833  self.ra = ra
834  self.dec = dec
835 
836 
838  """
839  Transform an (X, Y) in pixel coordinates on the detection image to (RA, DEC) in world coordinates.
840  Parameters
841  ----------
842  x : float
843  y : float
844 
845  Returns
846  -------
847  WorldCoordinate
848  """
849  global coordinate_system
850  wc = coordinate_system.image_to_world(cpp.ImageCoordinate(x-1, y-1)) # -1 as we go FITS -> internal
851  return WorldCoordinate(wc.alpha, wc.delta)
852 
853 
854 def get_sky_coord(x, y):
855  """
856  Transform an (X, Y) in pixel coordinates on the detection image to astropy SkyCoord.
857 
858  Parameters
859  ----------
860  x : float
861  y : float
862 
863  Returns
864  -------
865  SkyCoord
866  """
867  coord = pixel_to_world_coordinate(x, y)
868  sky_coord = SkyCoord(ra=coord.ra*u.degree, dec=coord.dec*u.degree)
869  return sky_coord
870 
871 
872 def radius_to_wc_angle(x, y, rad):
873  """
874  Transform a radius in pixels on the detection image to a radius in sky coordinates.
875 
876  Parameters
877  ----------
878  x : float
879  y : float
880  rad : float
881 
882  Returns
883  -------
884  Radius in degrees
885  """
886  return (get_separation_angle(x, y, x+rad, y) + get_separation_angle(x, y, x, y+rad)) / 2.0
887 
888 
889 def get_separation_angle(x1, y1, x2, y2):
890  """
891  Get the separation angle in sky coordinates for two points defined in pixels on the detection image.
892 
893  Parameters
894  ----------
895  x1 : float
896  y1 : float
897  x2 : float
898  y2 : float
899 
900  Returns
901  -------
902  Separation in degrees
903  """
904  c1 = get_sky_coord(x1, y1)
905  c2 = get_sky_coord(x2, y2)
906  return c1.separation(c2).degree
907 
908 
909 def get_position_angle(x1, y1, x2, y2):
910  """
911  Get the position angle in sky coordinates for two points defined in pixels on the detection image.
912 
913  Parameters
914  ----------
915  x1
916  y1
917  x2
918  y2
919 
920  Returns
921  -------
922  Position angle in degrees, normalized to -/+ 90
923  """
924  c1 = get_sky_coord(x1, y1)
925  c2 = get_sky_coord(x2, y2)
926  angle = c1.position_angle(c2).degree
927 
928  # return angle normalized to range: -90 <= angle < 90
929  return (angle + 90.0) % 180.0 - 90.0
930 
931 
933  """
934  Set the global coordinate system. This function is used internally by SourceXtractor++.
935  """
936  global coordinate_system
937  coordinate_system = cs
938 
939 
941  """
942  Convenience function for generating two dependent parameter with world (alpha, delta) coordinates
943  from image (X, Y) coordinates.
944 
945  Parameters
946  ----------
947  x : ParameterBase
948  y : ParameterBase
949 
950  Returns
951  -------
952  ra : DependentParameter
953  dec : DependentParameter
954 
955  See Also
956  --------
957  get_pos_parameters
958 
959  Examples
960  --------
961  >>> x, y = get_pos_parameters()
962  >>> ra, dec = get_world_position_parameters(x, y)
963  >>> add_output_column('mf_ra', ra)
964  >>> add_output_column('mf_dec', dec)
965  """
966  ra = DependentParameter(lambda x,y: pixel_to_world_coordinate(x, y).ra, x, y)
967  dec = DependentParameter(lambda x,y: pixel_to_world_coordinate(x, y).dec, x, y)
968  return (ra, dec)
969 
970 
971 def get_world_parameters(x, y, radius, angle, ratio):
972  """
973  Convenience function for generating five dependent parameters, in world coordinates, for the position
974  and shape of a model.
975 
976  Parameters
977  ----------
978  x : ParameterBase
979  y : ParameterBase
980  radius : ParameterBase
981  angle : ParameterBase
982  ratio : ParameterBase
983 
984  Returns
985  -------
986  ra : DependentParameter
987  Right ascension
988  dec : DependentParameter
989  Declination
990  rad : DependentParameter
991  Radius as degrees
992  angle : DependentParameter
993  Angle in degrees
994  ratio : DependentParameter
995  Aspect ratio. It has to be recomputed as the axis of the ellipse may have different ratios
996  in image coordinates than in world coordinates
997 
998  Examples
999  --------
1000  >>> flux = get_flux_parameter()
1001  >>> x, y = get_pos_parameters()
1002  >>> radius = FreeParameter(lambda o: o.get_radius(), Range(lambda v, o: (.01 * v, 100 * v), RangeType.EXPONENTIAL))
1003  >>> angle = FreeParameter(lambda o: o.get_angle(), Range((-np.pi, np.pi), RangeType.LINEAR))
1004  >>> ratio = FreeParameter(1, Range((0, 10), RangeType.LINEAR))
1005  >>> add_model(group, ExponentialModel(x, y, flux, radius, ratio, angle))
1006  >>> ra, dec, wc_rad, wc_angle, wc_ratio = get_world_parameters(x, y, radius, angle, ratio)
1007  >>> add_output_column('mf_world_angle', wc_angle)
1008  """
1009  ra = DependentParameter(lambda x,y: pixel_to_world_coordinate(x, y).ra, x, y)
1010  dec = DependentParameter(lambda x,y: pixel_to_world_coordinate(x, y).dec, x, y)
1011 
1012  def get_major_axis(x, y, radius, angle, ratio):
1013  if ratio <= 1:
1014  x2 = x + math.cos(angle) * radius
1015  y2 = y + math.sin(angle) * radius
1016  else:
1017  x2 = x + math.sin(angle) * radius * ratio
1018  y2 = y + math.cos(angle) * radius * ratio
1019 
1020  return x2, y2
1021 
1022  def get_minor_axis(x, y, radius, angle, ratio):
1023  if ratio <= 1:
1024  x2 = x + math.sin(angle) * radius * ratio
1025  y2 = y + math.cos(angle) * radius * ratio
1026  else:
1027  x2 = x + math.cos(angle) * radius
1028  y2 = y + math.sin(angle) * radius
1029 
1030  return x2, y2
1031 
1032  def wc_rad_func(x, y, radius, angle, ratio):
1033  x2, y2 = get_major_axis(x, y, radius, angle, ratio)
1034  return get_separation_angle(x, y, x2, y2)
1035 
1036  def wc_angle_func(x, y, radius, angle, ratio):
1037  x2, y2 = get_major_axis(x, y, radius, angle, ratio)
1038  return get_position_angle(x, y, x2, y2)
1039 
1040  def wc_ratio_func(x, y, radius, angle, ratio):
1041  minor_x, minor_y = get_minor_axis(x, y, radius, angle, ratio)
1042  minor_angle = get_separation_angle(x,y, minor_x, minor_y)
1043 
1044  major_x, major_y = get_major_axis(x, y, radius, angle, ratio)
1045  major_angle = get_separation_angle(x,y, major_x, major_y)
1046 
1047  return minor_angle / major_angle
1048 
1049  wc_rad = DependentParameter(wc_rad_func, x, y, radius, angle, ratio)
1050  wc_angle = DependentParameter(wc_angle_func, x, y, radius, angle, ratio)
1051  wc_ratio = DependentParameter(wc_ratio_func, x, y, radius, angle, ratio)
1052 
1053  return (ra, dec, wc_rad, wc_angle, wc_ratio)
1054 
def __init__(self, x_coord, y_coord, flux)
def __init__(self, normalization_factor=1)
def __init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle)
def get_flux_parameter(type=FluxParameterType.ISO, scale=1)
def __init__(self, init_value, range=Unbounded())
def __init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle, n)
def get_position_angle(x1, y1, x2, y2)
def print_parameters(file=sys.stderr)
def __init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle)
def __init__(self, param, value, sigma)
def add_prior(param, value, sigma)
def get_separation_angle(x1, y1, x2, y2)
def print_model_fitting_info(group, show_params=False, prefix='', file=sys.stderr)
def __init__(self, x_coord, y_coord, flux, effective_radius, aspect_ratio, angle)
def get_world_parameters(x, y, radius, angle, ratio)