wcssph2xy.pro

basic median filter simulation

;+
; NAME:
;     WCSSPH2XY 
; PURPOSE:
;     Convert spherical coordinates to x and y (map) angular coordinates
; EXPLANATION:
;     Convert spherical (longitude and latitude -- sky) coordinates to x 
;     and y (map) angular coordinates.  This procedure is the inverse of 
;     WCSXY2SPH.    See WCS_DEMO for example of use.
;
;     This is a lower level procedure -- given a FITS header, the user will 
;     usually use ADXY which will then call WCSSPH2XY with the appropriate 
;     parameters.
; CATEGORY:
;     Mapping and Auxiliary FITS Routine
;
; CALLING SEQUENCE:
;      wcssph2xy, longitude, latitude, x, y, [ map_type , CTYPE = ,
;               FACE =,PV2= , CRVAL = , CRXY = , LONGPOLE = ,
;               LATPOLE = , NORTH_OFFSET =, SOUTH_OFFSET =, BADINDEX =]
;
; INPUT PARAMETERS:
;     longitude - longitude of data, scalar or vector, in degrees 
;     latitude - latitude of data, same number of elements as longitude, 
;               in degrees
;     map_type - optional positional parameter, numeric scalar (0-26) 
;               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/ mu = 0
;   SIN     3    Orthographic               PV2_1,PV2_2 optional
;   STG     4    Stereographic              AZP w/ mu = 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
;   COP    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     convergence of inverse is poor
;                Spherical Cube
;   QSC    24    Quadrilateralized 
;                Spherical Cube
;   TSC    25    Tangential Spherical Cube
;   SZP    26    Slant Zenithal Projection   PV2_1,PV2_2, PV2_3 optional
;
; OPTIONAL INPUT 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.
;       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.
;       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 is [0,0]
;               This is not a FITS standard -- it is similar to CRPIX but in 
;               angular X,Y coordinates (degrees) rather than pixel coordinates
;       LATPOLE -  native latitude of the standard system's North Pole
;       LONGPOLE - native longitude of standard system's North Pole, default
;               is 180 degrees for Zenithal systems
;       NORTH_OFFSET - offset (radians) added to input points near north pole.
;       SOUTH_OFFSET - offset (radians) added to input points near south pole.
;       BADINDEX     - vector, list of transformed points too close to poles.
;
;
; OUTPUT PARAMETERS:
;
;       x - x coordinate of data, same number of elements as longitude, in 
;               degrees; if CRXY is set, then x will be returned offset by 
;               crxy(0).  NOTE: x in all map projections increases to the 
;               left, not the right.
;       y - y coordinate of data, same number of elements as longitude, in 
;               degrees; if CRXY is set, y will be returned offset by crxy[1]
;       bad - vector returning index to transformed points close to pole.
;
; OPTIONAL OUTPUT KEYWORD PARAMETERS:
;       FACE - a output variable used for spherical cube projections to 
;               designate the face of the cube on which the x and y 
;               coordinates lie.   Will contain the same number of elements as
;               X and Y.    Must contain at least 1 arbitrary element on input
;               If FACE is NOT defined on input, it is assumed that the
;               spherical cube projection is laid out over the whole sky
;               in the "sideways T" configuration.
; NOTES:
;       The conventions followed here are described in more detail in 
;       "Representations of Celestial Coordinates in FITS" by Mark Calabretta
;       and Eric Greisen (2002, A&A, 395, 1077; also  see  
;       http://www.aoc.nrao.edu/~egreisen).  The general 
;       scheme outlined in that article is to first use WCS_ROTATE to convert 
;       coordinates in one of three standard systems (celestial, galactic, 
;       or ecliptic) into a "native system" of latitude and longitude.  The 
;       latitude and longitude are then converted into x and y coordinates 
;       which depend on the map projection which is performed.   The rotation 
;       from standard to native coordinates can be skipped if one so desires.
;       This procedure necessitates two basic sections.  The first converts 
;       "standard" coordinates to "native" coordinates while the second converts
;       "native" coordinates to x and y coordinates.  The first section is 
;       simply a call to WCS_ROTATE, while the second contains the guts of 
;       the code in which all of the map projection is done.  This procedure 
;       can be called in a form similar to AITOFF, EQPOLE, or QDCB by calling 
;       wcssph2xy with a fifth parameter specifying the map projection by 
;       number and by not using any of the keywords related to the map 
;       projection type (e.g. CTYPE).
;
; 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 to a map projection
;       type).  All variables are converted into double precision values and 
;       angular measurements are converted from degrees into radians.
;       If necessary, longitude values are converted into the range -pi to pi.
;       Any latitude points close to the  of the poles are mapped to a specific
;       latitude of  from the pole so that the map transformations become
;       completely invertible.  The magnitude of this correction is given by 
;       the keywords NORTH_OFFSET and SOUTH_OFFSET and a list of affected 
;       points is optionally returned in the "badindex" output parameter.
;       The next task of the procedure is to convert the "standard" 
;       coordinates to "native" coordinates by rotating the coordinate system.
;       This rotation is performed by the procedure WCS_ROTATE and is governed
;       by the keywords CRVAL and LONGPOLE.   The final task of the WCSSPH2XY 
;       is to take "native" latitude and longitude coordinates and convert 
;       them into x and y coordinates.  Any 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). 
;
;       Note that a further transformation (using the CD matrix) is required
;       to convert the (x,y) coordinates to pixel coordinates. 
; COMMON BLOCKS:
;
;       none
;
; PROCEDURES CALLED:
;       WCS_ROTATE
;
; COPYRIGHT NOTICE:
;
;       Copyright 1993, The Regents of the University of California. This
;       software was produced under U.S. Government contract (W-7405-ENG-36)
;       by Los Alamos National Laboratory, which is operated by the
;       University of California for the U.S. Department of Energy.
;       The U.S. Government is licensed to use, reproduce, and distribute
;       this software. Neither the Government nor the University makes
;       any warranty, express or implied, or assumes any liability or
;       responsibility for the use of this software.
;
; AUTHOR:
;
;       Rick Balsano
;
; MODIFICATIONS/REVISION LEVEL:
;
;       1.1     8/31/93  
;       2.3     9/15/93  W. Landsman (HSTX) Update quad cube coords, vectorize
;                        keywords
;       2.4     12/29/93 I. Freedman (HSTX) Eliminated LU decomposition
;       2.5     1/5/93   I. Freedman (HSTX) Offset keywords / bad point index
;       2.6     Dec 94   Compute pole for transformations where the reference
;                       pixel is at the native origin    W. Landsman (HSTX)
;       2.7     May 95  Change internal variable BETA for V4.0 compatibility
;       2.8     June 95 Change loop indices from integer to long
;       2.9     3/18/96 Change FACE usage for cube projections to match WCSLIB
;                       C/FORTRAN software library.
;       2.10    02/18/99 Fixed implementation of ARC algorithm
;       2.11    June 2003 Update conic projections, add LATPOLE keyword
;	2.12	Aug 2003, N.Rich - Fix pre-V5.5 bug from previous update
;       2.13    Sep 2003, W. Landsman CTYPE keywords need not be 8 characters
;       2.14    Jan 2004, W. Landsman don't modify scalars, fix PARabolic code
;       2.15    Feb 2004, W. Landsman Fix AZP and AIR algorithms
;       3.0    May 2004  W. Landsman Support extended SIN (=NCP), slant zenithal
;                  (SZP), and zenithal polynomail (ZPN) projections, use
;                   PV2 keyword vector instead of PROJP1, PROJP2
;       3.1     Jul 2005 W.Landsman/C. Markwardt Set unprojectable points in
;                   tangent projection to NaN
;       3.1.1   Jul 2005 Fixed 3.1 mod to work for scalars
;       3.2     Dec 2005 Fixed Airy projection for latitude centered at 90 deg
;       3.3     Aug 2007 R. Munoz, W.Landsman Correct treatment of PV1_2 and 
;                        PV2_2 parameters
;       3.4    Oct 2007  Sergey Koposov Support HEALPIX projection
;-

PRO wcssph2xy,longitude,latitude,x,y,map_type, ctype=ctype,$
              face=face, pv2 = pv2,$
              crval=crval,crxy=crxy,longpole=longpole, latpole= latpole, $
              north_offset=north_offset, south_offset=south_offset, $
              badindex=badindex

; DEFINE ANGLE CONSTANTS 
 pi = !DPI
 pi2 = pi/2.d0
 radeg = 57.295779513082323d0
 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 - WCSSPH2XY, longitude, latitude, x, y, [ map_type,' 
    print,'           CTYPE= ,FACE=, PV2_1=, PV2_2=, CRVAL=, CRXY=, LATPOLE='
    print,'           LONGPOLE= ,NORTH_OFFSET=, SOUTH_OFFSET=, BADINDEX=]'
    return
 endif 
  
   
; GENERAL ERROR CHECKING
; find the number of elements in each of the data arrays

 n_long = N_elements( longitude )
 n_lat = N_elements( latitude )

; check to see that the data arrays have the same size
 if (n_long ne n_lat) then begin
     message,$
       'LONGITUDE and LATITUDE must have the same number of elements.'
 endif

 if (N_params() eq 5) then begin

  if (keyword_set(ctype1) or keyword_set(ctype2)) then message,$
'Use either the MAP_TYPE positional parameter or set the projection type with CRVAL1 and CRVAL2, but not both.'
; set projection_type string using map_type parameter (a number) 
  if ((map_type ge 0) and (map_type le 25)) then begin
      projection_type=map_types[map_type]
  endif else message,'MAP_TYPE must be >= 0 and <= 25, 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 = ctype[0]        
        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 

; 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 begin 
           projection_type='CAR'
        message,'Projection type not supplied, set to default (Cartesian)',/INF
 endif

; Check to make sure all the correct parameters and keywords are set for 
; spherical projections.
if (keyword_set(ctype3) or keyword_set(face) or (projection_type eq 'CSC') or $
    (projection_type eq 'QSC') or (projection_type eq 'TSC')) then begin

  if (n_elements(face) eq 0) 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 the x and y offsets are not set then assume they are zero.
if (((n_elements(crx) eq 1) and (n_elements(cry) eq 0)) or $
     ((n_elements(crx) eq 0) and (n_elements(cry) eq 1))) then $
    message,'If either CRX or CRY is set, the other must also be set'

if (((n_elements(crval1) eq 1) and (n_elements(crval2) eq 0)) or $
    ((n_elements(crval1) eq 0) and (n_elements(crval1) eq 1))) then $
    message,'If either CRVAL1 or CRVAL2 is set, the other must also be set'


; Convert all longitude values into the range -180 to 180 so that equations
; work properly.
  lng = double( longitude )   & lat = double( latitude )
  temp = where(lng ge 180.0, Ntemp)
  if Ntemp GT 0 then lng[temp] = lng[temp] - 360.0d0

; Make small offsets at poles to allow the transformations to be 
; completely invertible.  These generally introduce a small fractional error
; but only at the poles.  They are necessary since all maps 
; lose information at the poles when a rotation is applied, because all points 
; within NORTH_ or SOUTH_OFFSET of the poles are mapped to the same points.

  IF N_elements(north_offset) EQ 0 then north_offset = 1.d-7
  IF N_elements(south_offset) EQ 0 then south_offset = 1.d-7

  bad = where(abs(lat - 90.0) lt north_offset*radeg)
  IF (bad[0] ne -1) THEN BEGIN
      ;; commented out this message
  ;;MESSAGE,/INFORM,'Some input points are too close to the NORTH pole.'
   lat[bad] = 90.0 - north_offset*RADEG
   IF KEYWORD_SET(badindex) THEN badindex = bad
  ENDIF