wcsxy2sph.pro

basic median filter simulation

;+
; NAME:
;      WCSXY2SPH  
;
; PURPOSE:
;      Convert x and y (map) coordinates to spherical coordinates
; EXPLANATION:
;      To convert x and y (map) coordinates to spherical (longitude and 
;      latitude or sky) coordinates.    This procedure is the inverse of
;      WCSSPH2XY.
;
;     This is a lower level procedure -- given a FITS header, the user will 
;     usually use XYAD which will then call WCSXY2SPH with the appropriate 
;     parameters.
; CATEGORY:
;      Mapping and Auxilary FITS Routine
;
; CALLING SEQUENCE:
;
;      wcsxy2sph, x, y, longitude, latitude, [map_type], [ CTYPE = ,$
;             FACE = ,PV2 = ,CRVAL =, CRXY =, LONGPOLE=, LATPOLE=]
;
; INPUT PARAMETERS:
;
;       x - x coordinate of data, scalar or vector, in degrees, NOTE: x 
;               increases to the left, not the right
;       y - y coordinate of data, same number of elements as x, in degrees
;       map_type - optional positional parameter, scalar corresponding to a 
;               particular map projection.  This is not a FITS standard, it is 
;               simply put in to allow function similar to that of less general 
;               map projection procedures (eg AITOFF).  The following list gives
;               the map projection types and their respective numbers.
;
;  FITS  Number  Name                       Comments
;  code   code
;  ----  ------  -----------------------    -----------------------------------
;   DEF     0    Default = Cartesian
;   AZP     1    Zenithal perspective       pv2_1 required
;   TAN     2    Gnomic                     AZP w/ pv2_1 = 0
;   SIN     3    Orthographic               pv2_1, pv2_2 optional
;   STG     4    Stereographic              AZP w/ pv2_1 = 1
;   ARC     5    Zenithal Equidistant
;   ZPN     6    Zenithal polynomial        PV2_0, PV2_1....PV2_20 possible
;   ZEA     7    Zenithal equal area
;   AIR     8    Airy                       pv2_1 required
;   CYP     9    Cylindrical perspective    pv2_1 and pv2_2 required
;   CAR    10    Cartesian
;   MER    11    Mercator
;   CEA    12    Cylindrical equal area     pv2_1 required
;   xy    13    Conical perspective        pv2_1 and pv2_2 required
;   COD    14    Conical equidistant        pv2_1 and pv2_2 required
;   COE    15    Conical equal area         pv2_1 and pv2_2 required
;   COO    16    Conical orthomorphic       pv2_1 and pv2_2 required
;   BON    17    Bonne's equal area         pv2_1 required
;   PCO    18    Polyconic
;   SFL    19    Sanson-Flamsteed
;   PAR    20    Parabolic
;   AIT    21    Hammer-Aitoff
;   MOL    22    Mollweide
;   CSC    23    Cobe Quadrilateralized     inverse converges poorly
;                Spherical Cube
;   QCS    24    Quadrilateralized
;                Spherical Cube
;   TSC    25    Tangential Spherical Cube
;   SZP    26    Slant Zenithal perspective  PV2_1,PV2_2, PV2_3 optional 
;
; OPTIONAL KEYWORD PARAMETERS:
;
;       CTYPE - One, two, or three element vector containing 8 character 
;               strings corresponding to the CTYPE1, CTYPE2, and CTYPE3 
;               FITS keywords: 
;
;               CTYPE[0] - first four characters specify standard system
;               ('RA--','GLON' or 'ELON' for right ascension, galactic 
;               longitude or ecliptic longitude respectively), second four 
;               letters specify the type of map projection (eg '-AIT' for 
;               Aitoff projection)
;               CTYPE[1] - first four characters specify standard system
;               ('DEC-','GLAT' or 'ELAT' for declination, galactic latitude
;               or ecliptic latitude respectively; these must match 
;               the appropriate system of ctype1), second four letters of 
;               ctype2 must match second four letters of ctype1.
;               CTYPE[2] - if present must be the 8 character string,'CUBEFACE',
;                only used for spherical cube projections to identify an axis 
;               as containing the face on which each x and y pair of 
;               coordinates lie.
;       FACE - a input variable used for spherical cube projections to 
;               designate the face of the cube on which the x and y 
;               coordinates lie.   Must contain the same number of elements
;               as X and Y.
;       CRVAL - 2 element vector containing standard system coordinates (the 
;               longitude and latitude) of the reference point
;       CRXY - 2 element vector giving the x and y coordinates of the 
;               reference point, if this is not set the offset of the x 
;               coordinate is assumed to be 0.
;       LATPOLE -  native latitude of the standard system's North Pole
;       LONGPOLE - native longitude of standard system's North Pole, default
;               is 180 degrees, numeric scalar
;       PV2  - Vector of projection parameter associated with latitude axis
;             PV2 will have up to 21 elements for the ZPN projection, up to 3
;             for the SIN projection and no more than 2 for any other
;             projection.   The first element corresponds to PV2_1, the 
;             second to PV2_2, etc.
;
; OUTPUT PARAMETERS:
;
;       longitude - longitude of data, same number of elements as x, in degrees
;       latitude - latitude of data, same number of elements as x, in degrees
;
;       Longitude and latitude will be set to NaN, wherever elements of X,Y
;       have no corresponding longitude, latitude values. 
; NOTES:
;       The conventions followed here are described in more detail in the paper
;      "Representations of Celestial Coordinates in FITS" by Calabretta &
;       Greisen (2002, A&A, 395, 1077, also see 
;       http://www.aoc.nrao.edu/~egreisen).   The general scheme 
;       outlined in that article is to convert x and y coordinates into a 
;       "native" longitude and latitude and then rotate the system into one of 
;       three generally recognized systems (celestial, galactic or ecliptic).
;
;       This procedure necessitates two basic sections.  The first converts 
;       x and y coordinates to "native" coordinates while the second converts 
;       "native" to "standard" coordinates.  The first section contains the 
;       guts of the code in which all of the map projection is done.  The 
;       second step is performed by WCS_ROTATE and only involves rotation of 
;       coordinate systems.  WCSXY2SPH can be called in a form similar to 
;       AITOFF, EQPOLE, or QDCB by calling wcsxy2sph with a fifth parameter 
;       specifying the map projection by number and by not using any of the 
;       keywords related to the map projection type (eg ctype1 and ctyp2).
;
; PROCEDURE:
;       The first task of the procedure is to do general error-checking to 
;       make sure the procedure was called correctly and none of the 
;       parameters or keywords conflict.  This is particularly important 
;       because the procedure can be called in two ways (either using 
;       FITS-type keywords or using a number corresponding a map projection 
;       type).  All variables are converted into double precision values.
;
;       The second task of the procedure is to take x and y coordinates and 
;       convert them into "native" latitude and longitude coordinates.  
;       Map-specific error-checking is done at this time.  All of the 
;       equations were obtained from "Representations of Celestial 
;       Coordinates in FITS" and cases needing special attention are handled 
;       appropriately (see the comments with individual map projections for 
;       more information on special cases).     WCS_ROTATE is then called to 
;       convert the "native" coordinates to "standard" coordinates by rotating
;       the coordinate system.  This rotation is governed by the keywords 
;       CRVAL, and LONGPOLE.  The transformation is a straightforward 
;       application of euler angles.  Finally, longitude values are converted 
;       into the range from 0 to 360 degrees.
;
; COMMON BLOCKS:
;       none
; PROCEDURES CALLED:
;       WCS_ROTATE
;
; AUTHOR:
;
;       Rick Balsano
;
; MODIFICATIONS/REVISION LEVEL:
;
; 1.1    8/31/93
; 1.2    9/12/93   W. Landsman Vectorized CRXY, CRVAL, CTYPE
; 1.3    29/12/93  I. Freedman Eliminated LU decomposition
; 1.4    22/09/94  W. Landsman If scalar input, then scalar output
; 1.5    02/03/05  W. Landsman Change variable name BETA for V4.0 compatibility
; 1.6    06/07/05  W. Landsman Change loop index from integer to long
; 1.7    02/18/99  W. Landsman Fixed implementation of ARC algorithm
; 1.8    June 2003 W. Landsman Update conic projections, add LATPOLE keyword
; 1.81   Sep 2003 W. Landsman Avoid divide by zero 
; 1.82   Sep 2003 W. Landsman CTYPE keywords need not be 8 characters
; 1.83   Sep 2003 W. Landsman Preserve input array sizes
; 1.9    Jan 2004 W. Landsman don't modify scalars, fix PARabolic code
; 2.0    Feb 2004 W. Landsman Fix AIR and AZP projections
; 2.1    Feb 2004 W. Landsman Fix tangent projection for matrix input
; 3.0    May 2004  W. Landsman Support extended SIN (=NCP), slant zenithal
;                  (SZP), and zenithal polynomial (ZPN) projections, use
;                   PV2 keyword vector instead of PROJP1, PROJP2
; 3.1    May 2004 W. Landsman/J. Ballet Handle NaN values, flag invalid output
;                   for AITOFF projection
; 3.1.1  Dec 2005 W. Landsman/W. Thompson Fixed problem with Airy projection
;                   centered on 90 degree latitude
; 3.1.2  May 2006 W. Landsman/Y.Sato Fix problem selecting the correct root
;                    for the ZPN projection
; 3.2    Aug 2007  W. Landsman Correct treatment of PVi_j parameters
; 3.3    Oct 2007  Sergey Koposov Support HEALPIX projection
;-

PRO wcsxy2sph, x, y, longitude, latitude, map_type, ctype=ctype, $
              face=face,pv2 = pv2,$
              crval=crval,crxy = crxy,longpole=longpole, Latpole = latpole

; Define angle constants 

 compile_opt idl2
 pi = !DPI
 radeg = 57.295779513082323d0
 pi2 = pi/2.0d
 map_types=['DEF','AZP','TAN','SIN','STG','ARC','ZPN','ZEA','AIR','CYP',$
            'CAR','MER','CEA','COP','COD','COE','COO','BON','PCO','SFL',$
            'PAR','AIT','MOL','CSC','QSC','TSC','SZP']

; check to see that enough parameters (at least 4) were sent
 if ( N_params() lt 4 ) then begin
    print,'Syntax - WCSXY2SPH, x, y, longitude, latitude,[ map_type,  '
    print,'             CTYPE= , FACE=, PV2= , CRVAL= , CRXY= , '
    print,'             LATPOLE = , LONGPOLE = ]'
    return
 endif else if (n_params() eq 5) then begin
    if keyword_set(ctype) then message,$
  'Use either the MAP_TYPE positional parameter or set the projection type' + $
    'with CTYPE, but not both.'

; set projection_type string using map_type parameter (a number) 
  if (n_elements(map_type) ne 0) then begin 
         projection_type = map_types[map_type]
  endif else message,'MAP_TYPE must be >= 0 and <= 26, it was set to '+map_type

endif else if (n_params() eq 4) then begin

; check to see that CTYPE is set correctly 

  if N_elements( ctype ) GE 1 then begin
        ctype1 = strtrim(ctype[0],2)     
        if strlen(ctype1) LT 8 then message,'ERROR - ' + strupcase(ctype1) + $
               ' is not a valid spherical projection type.'  
        projection_type = strupcase(strmid(ctype1,5,3))
  endif

  if N_elements( ctype ) GE 2  then begin
        ctype2 = ctype[1]
         if (projection_type ne strupcase(strmid(ctype2,5,3))) then begin
      message,'The same map projection type must be in characters',/continue
      message,'     5-8 of both CTYPE1 and CTYPE2.',/con
      return
    endif
    if (((strupcase(strmid(ctype1,1,2)) eq 'RA') and $
         (strupcase(strmid(ctype2,1,3)) ne 'DEC')) or $
        ((strupcase(strmid(ctype1,1,4)) eq 'GLON') and $
         (strupcase(strmid(ctype2,1,4)) ne 'GLAT')) or $
        ((strupcase(strmid(ctype1,1,4)) eq 'ELON') and $
         (strupcase(strmid(ctype2,1,4)) ne 'ELAT'))) $
    then begin message,'The same standard system must be in the first 4',$
                       /continue
      print,'characters of both CTYPE1 and CTYPE2.'
      return
    endif
  endif else projection_type = 'DEF'

endif 

; GENERAL ERROR CHECKING

; find the number of elements in each of the data arrays
 n_x = n_elements(x)
 n_y = n_elements(y)
 sz_x = size(x)
 sz_y = size(y)

; check to see that the data arrays have the same size
 if (n_x ne n_y) then $
        message,'The arrays X and Y must have the same number of elements.'

; this sets the default map projection type for the cases when map_type or
; projection_type is set to 'DEF' or if projection_type is not set at this
; point.  As suggested in 'Representations of Celestial Coordinates in FITS'
; the default type is set to CAR (Cartesian) the simplest of all projections.
if ((n_elements(projection_type) eq 0) or (projection_type eq 'DEF')) then $
   projection_type='CAR'
 
; Check to make sure all the correct parameters and keywords are set for
; spherical projections.
 if ((N_elements(ctype) EQ 3) or keyword_set(face) or  $
    (projection_type eq 'CSC') or $
    (projection_type eq 'QSC') or (projection_type eq 'TSC')) then begin

  if not(keyword_set(face)) then noface = 1 else noface = 0

endif

; check to see if the x and y offsets are set properly.  If not, break out
; of program.  If so, apply offsets.  If the x and y offsets are not set,
; then assume they are zero.

 if ( (keyword_set(crxy)) and N_elements(crxy) NE 2) then $
     message,'Offset keyword CRXY must contain 2 elements'

 if keyword_set(crxy) then begin
        xx = double(x - crxy[0] )
        yy = double(y - crxy[1] )
 endif else begin
        xx = double(x)
        yy = double(y)
 endelse

 if  ( N_elements(crval) eq 1 ) then $
        message,'CRVAL keyword must be a 2 element vector'


; BRANCH BY MAP PROJECTION TYPE
case strupcase(projection_type) of
  'AZP':begin
  PV2_1 = N_elements(PV2) GE 1 ? PV2[0] : 0.0  ;PV2_1 =mu (spherical radii)
  PV2_2 = N_elements(PV2) GE 2 ? PV2[1] : 0.0   ; PV2_2 = gamma (degrees)

    if (pv2_1 lt 0) then message,$
      'AZP map projection requires the keyword pv2_1 >= 0'
    gamma = pv2_2/radeg
    mu = pv2_1
    r = sqrt(xx^2 + yy^2*cos(gamma)^2)
    rho = r/(radeg*(mu+1) + yy*sin(gamma) )
    omega = asin( rho*mu/ sqrt( rho^2 + 1.d0) )
    xsi = atan(1.d0, rho)
    phi = atan(xx, -yy*cos(gamma) )
    theta1 = xsi - omega
    theta2 = xsi + omega + !dpi
    theta = theta1*0.0
    if abs(mu) LT 1 then begin
          g = where(abs(theta1) LT pi2, Ng)
          if Ng GT 0 then theta[g] = theta1[g]
          g = where(abs(theta2) LT pi2, Ng)
          if Ng GT 0 then theta[g] = theta2[g]
    endif else begin
          diff1 = abs(pi2 - theta1)