wcssph2xy.pro

basic median filter simulation

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

; Convert from standard coordinate system to "native" coordinate system
; if the CRVAL keyword is set.  Otherwise, assume the latitude and longitude 
; given are in "native" coordinates already (this is  essentially what is done
; in the procedure AITOFF).

 PV2_1 = N_elements(pv2) GT 0 ? pv2[0] : 0
 PV2_2 = N_elements(pv2) GT 1 ? pv2[1] : 0

 if N_elements(crval) GE 2 then begin
           if N_elements(map_type) EQ 0 then begin
           	wmt      = where(projection_type EQ map_types)
		map_type = wmt[0]
	   endif
           conic = (map_type GE 13) and (map_type LE 16) 
           zenithal = ((map_type GE 1) and (map_type LE 8)) or (map_type EQ 26) 

; Rotate from standard celestial coordinates into the native system.
        if conic then theta0 = PV2_1 else if zenithal then theta0 = 90 $
                 else theta0 = 0                 
        wcs_rotate, lng, lat, phi, theta, crval, $
                latpole = latpole, longpole=longpole, theta0 = theta0
        phi = phi/radeg
        theta = theta/radeg
 endif else begin
        phi = lng/radeg
        theta = lat/radeg
 endelse

; BRANCH BY MAP PROJECTION TYPE
case strupcase(projection_type) of
  'AZP':begin
     if (PV2_1 lt 0) then message,$
      'AZP map projection requires the keyword PV2_1 >= 0'
    gamma = PV2_2/radeg
    mu = PV2_1

    r_theta = radeg*cos(theta)*(mu + 1.d0)/ $
             ( (mu + sin(theta)) + cos(theta)*cos(phi)*tan(gamma))
    x = r_theta*sin(phi)
    y = -r_theta*cos(phi)/cos(gamma)
  end  
  'SZP': begin
     mu = N_elements(PV2) GT 0 ? PV2[0] : 0
     phi_c = N_elements(PV2) GT 1 ? PV2[1] : 0
     theta_c = N_elements(PV2) GT 1 ? PV2[2] : 90
     phi_c = phi_c/radeg & theta_c = theta_c/radeg
     xp = -mu*cos(theta_c)*sin(phi_c)
     yp =  mu*cos(theta_c)*cos(phi_c)
     zp =  mu*sin(theta_c) + 1.
     denom = zp - (1-sin(theta))
     x = radeg*( zp*cos(theta)*sin(phi) - xp*(1-sin(theta)) )/ denom
     y = -radeg*( zp*cos(theta)*cos(phi) + yp*(1-sin(theta)) )/ denom

     end
  'TAN':begin
    sz_theta = size(theta,/dimen)
    if sz_theta[0] EQ 0 then x = !Values.D_NAN else $
          x = make_array(value = !values.D_NAN, dimen=sz_theta)
    y = x
    g = where(theta GT 0, Ng)
    if Ng GT 0 then begin
        r_theta = radeg/tan(theta[g])
        x[g] = r_theta*sin(phi[g])
        y[g] = -r_theta*cos(phi[g])
    endif
  end

  'SIN':begin
    if N_elements(PV2_1) EQ 0 then PV2_1 = 0
    if N_elements(PV2_2) EQ 0 then PV2_2 = 0
    if (PV2_1 EQ 0) and (PV2_2 EQ 0) then begin
        r_theta = radeg*cos(theta)
        x = r_theta*sin(phi)
        y = -r_theta*cos(phi)
    endif else begin                   ;NCP projection
        x =  radeg*(cos(theta)*sin(phi) + PV2_1*(1-sin(theta)) )
        y = -radeg*(cos(theta)*cos(phi) - PV2_2*(1-sin(theta)) )
    endelse
  end

  'STG':begin
    r_theta = 2.d0*radeg*tan((pi2-theta)/2.d0)
    x = r_theta*sin(phi)
    y = -r_theta*cos(phi)
  end

  'ARC':begin
    r_theta = radeg*( pi2 - theta )
    x = r_theta*sin(phi)
    y = -r_theta*cos(phi)
  end

  'ZPN':begin
    z = pi2 - theta
    g = where(pv2 NE 0, Ng)
    if Ng GT 0 then np = max(g) else np =0
    r_theta = radeg*poly(z, pv2[0:np])
    x = r_theta*sin(phi)
    y = -r_theta*cos(phi)
    end

  'ZEA':begin
    r_theta = 2.d0*radeg*sin((pi2 - theta)/2.d0)
    x = r_theta*sin(phi)
    y = -r_theta*cos(phi)
  end

  'AIR':begin
    if not(keyword_set(PV2_1)) then begin
      message,/informational,$
          'PV2_1 not set, using default of PV2_1 = 90 for AIR map projection'
      PV2_1 = 9.d1
    endif
    theta_b = PV2_1/radeg

    xi = (pi2 - theta)/2.d0

; When theta_b (aka PV2_1 in radians) is equal to pi/2 the normal equations
; for the AIR projection produce infinities.  To avoid the problem, values
; of theta_b equal to pi/2 cause a different set of equations to be used.
    if (theta_b eq pi2) then begin

; AIR produces the same radii for different latitudes, causing some overlap.  To
; avoid this problem, if latitudes which are far enough south to be a problem 
; are included in the data, the routine will stop.

      if (min(theta) lt -36/radeg) then begin
        message,'AIR produces overlap of native latitudes south of ',/continue
        print,'-36 with the PV2_1 = 90'
        return
      endif

; points with xi too small are labelled as bad to prevent poor behavior of the
; equation for r_theta
      good = where(abs(xi) ge 1.d-10, Ngood)
      r_theta = lng*0
      if (Ngood GT 0) then $
        r_theta[good] = -2*radeg*(alog(cos(xi[good]))/tan(xi[good]) - $
	                 0.5*tan(xi[good]))

    endif else begin
      xi_b = (pi2 - theta_b)/2.d0
      a = alog(cos(xi_b))/tan(xi_b)/tan(xi_b)

; AIR produces the same radii for different latitudes, causing some overlap.  To
; avoid this problem, if latitudes which are far enough south to be a problem 
; are included in the data, the routine will stop.

      xi_temp = (findgen(90) + 1)/radeg
      radius=-radeg*(alog(cos(xi_temp))/tan(xi_temp)+alog(cos(xi_b))/$
                                                      tan(xi_b)*tan(xi_temp))
      i = 0
      repeat i = i + 1 $
      until ((radius[i + 1] le radius[i]) or (i eq n_elements(radius) - 2))
      if (i lt (n_elements(radius)- 2)) then min_lat = 90 - 2*radeg*xi_temp[i] $
      else min_lat = -90
      if (min(theta) lt min_lat[0]/radeg) then begin
        message,'AIR produces overlap of native latitudes south of ',/continue
        print,format='(i3,a21,i3)',min_lat[0],' with the PV2_1 = ',PV2_1
        return
      endif

; points with xi too small are labelled as bad to prevent poor behavior of the
; equation for r_theta
 
      good = where(abs(xi) ge 1.d-10, Ngood)
      r_theta = lng*0
      if (Ngood GT 0) then r_theta[good] = -2*radeg*(alog(cos(xi[good]))/$
        tan(xi[good]) + a*tan(xi[good]))
    endelse
    x = r_theta*sin(phi)
    y = -r_theta*cos(phi)
  end

  'CYP':begin
    if (n_elements(PV2_1) eq 0) then begin
      message,/informational,$
           'PV2_1 not set, using default of PV2_1 = 0 for CYP map projection'
      PV2_1 = 0.d0
    endif
    if (n_elements(PV2_2) eq 0) then begin
      message,/informational,$
           'PV2_2 not set, using default of PV2_2 = 1 for CYP map projection'
      PV2_2 = 1.d0
    endif
    if (PV2_1 eq -PV2_2) then message,$
      'PV2_1 = -PV2_2 is not allowed for CYP map projection.'

    x = PV2_2*radeg*phi
    y = radeg*(PV2_1 + PV2_2)*sin(theta)/(PV2_1 + cos(theta))
  end
  
  'CAR':begin
    x = radeg*phi
    y = radeg*theta
  end
  
  'MER':begin
    x = radeg*phi
    y = radeg*alog(tan((pi2 + theta)/2.d0))
  end
  
  'CEA':begin
    if N_elements(PV2_1) EQ 0  then message,$
      'CEA map projection requires that PV2_1 keyword be set.'
    if ((PV2_1 le 0) or (PV2_1 gt 1)) then message,$
      'CEA map projection requires 0 < PV2_1 <= 1'
    x = radeg*phi
    y = radeg*sin(theta)/PV2_1
  end
  
  'COP':begin
    if not(keyword_set(PV2_1)) then message,$
      'COP map projection requires that PV2_1 keyword be set.'
    if not(keyword_set(PV2_2)) then begin 
      message,/informational,$
      'PV2_2 not set, using default of PV2_2 = 0 for COP map projection'
      PV2_2= 0
    endif
    if ((PV2_1 lt -90) or (PV2_2 gt 90) or (PV2_1 gt 90)) then message,$
 'PV2_1 and PV2_2 must satisfy -90<=PV2_1<=90,0<=PV2_2<=90 for COP projection'
    if (PV2_1 eq -PV2_2) then message,$
 'COP projection with PV2_1=-PV2_2 is better done as a cylindrical projection'
    theta_a = PV2_1/radeg
    alpha = PV2_2/radeg
    bad = where((theta ge theta_a + pi2) or (theta le theta_a - pi2))
    if (bad[0] ne -1) then begin 
      message,/continue,$
  'COP map projection diverges for native latitude = PV2_1 +- 90.'
      message,'Remove these points and try again.'
    endif

    r_theta = radeg*cos(alpha)*(1.d0/tan(theta_a)-tan(theta-theta_a))
    a_phi = phi*sin(theta_a)
    y_0 = radeg*cos(alpha)/tan(theta_a)
    x = r_theta*sin(a_phi)
    y = y_0 - r_theta*cos(a_phi)

  end
  
  'COD':begin
    if not(keyword_set(PV2_1)) then message,$
      'COD map projection requires that PV2_1 keyword be set.'
    if not(keyword_set(PV2_2)) then begin
      message,/informational,$
     'PV2_2 not set, using default of PV2_2 = 0 for COD map projection'
      PV2_2 = 0
    end
    if ((PV2_1 lt -90) or (PV2_2 gt 90) or (PV2_1 gt 90)) then message,$
 'PV2_1 and PV2_2 must satisfy -90<=PV2_1<=90,PV2_2<=90 for COD projection'
    if (PV2_1 eq -PV2_2) then message,$
    'COD gives divergent equations for PV2_1 = -PV2_2'
    theta_a = PV2_1/radeg

; when PV2_1 not = PV2_2 use regular equations
  if (PV2_2 NE 0) then begin
      alpha = PV2_2/radeg
      r_theta = theta_a - theta + alpha/(tan(alpha)*tan(theta_a))
      a_phi = sin(theta_a)*sin(alpha)*phi/alpha
      y_0 = radeg*alpha/(tan(alpha)*tan(theta_a))
; if the two parameters PV2_1 and PV2_2 are equal use the simpler set of
; equations
    endif else begin 
      r_theta = theta_a - theta + 1.d0/tan(theta_a)
      a_phi = phi*sin(theta_a)
      y_0 = radeg/tan(theta_a)

    endelse
    x = radeg*r_theta*sin(a_phi)
    y = y_0 - radeg*r_theta*cos(a_phi)
    
  end
  
  'COE':begin
    if N_elements(PV2_1) EQ 0 then message,$
      'COE map projection requires that PV2_1 keyword be set.'
    if N_elements(PV2_2) EQ 0 then begin
      message,/informational,$
      'PV2_2 not set, using default of PV2_2 = 0 for COE map projection'
      PV2_2 = 0
    end
    if ((PV2_1 lt -90) or (PV2_2 gt 90) or (PV2_1 gt PV2_2)) then message,$
 'PV2_1 and PV2_2 must satisfy -90<=PV2_1<=PV2_2<=90 for COE map projection'
    if (PV2_1 eq -PV2_2) then message,$
    'COE gives divergent equations for PV2_1 = -PV2_2'
   
    theta_1 = (PV2_1 - PV2_2)/radeg
    theta_2 = (PV2_1 + PV2_2)/radeg
    s_1 = sin(theta_1)
    s_2 = sin(theta_2)
    stheta_a = sin(PV2_1/radeg)
    gamma = s_1 + s_2
    r_theta=radeg*2.d0*sqrt(1.d0+ s_1*s_2-gamma*sin(theta))/gamma

     a_phi = phi*gamma/2.d0
    y_0 = radeg*2.d0*sqrt(1.d0+ s_1*s_2-gamma*stheta_a)/gamma
    x = r_theta*sin(a_phi)
    y = y_0 - r_theta*cos(a_phi)
  end
  
  'COO':begin
    if not(keyword_set(PV2_1)) then message,$
      'COO map projection requires that PV2_1 keyword be set.'
    if not(keyword_set(PV2_2)) then begin
      message,/informational,$
      'PV2_2 not set, using default of PV2_2 = 0 for COO map projection'
      PV2_2 = 0
    end
    if ((PV2_1 lt -90) or (PV2_2 gt 90) or (PV2_1 gt 90)) then message,$
 'PV2_1 and PV2_2 must satisfy -90<=PV2_1<=90,PV2_2<=90 for COO projection'
    if (PV2_1 eq -PV2_2) then message,$
    'COO gives divergent equations for PV2_1 = -PV2_2'
    theta_1 = (PV2_1 - PV2_2)/radeg
    theta_2 = (PV2_1 + PV2_2)/radeg 
    theta_a = PV2_1/radeg


; for cases where PV2_1 = 0, use a simpler formula to calculate c,
; otherwise use the regular formula
    if (PV2_2 eq 0) then c = sin(theta_1) else $
    c = alog(cos(theta_2)/cos(theta_1))/alog(tan((pi2-theta_2)/2.d0)/$
    tan((pi2-theta_1)/2.d0))

    alpha = radeg*cos(theta_1)/(c*(tan((pi2-theta_1)/2.d0))^c)
    r_theta = alpha*(tan((pi2-theta)/2.d0))^c
    y_0 = alpha*tan((pi2-theta_a)/2.)^c 
    a_phi = c*phi
    x = r_theta*sin(a_phi)
    y = y_0 - r_theta*cos(a_phi)

  end
 
  'BON':begin
    if (N_elements(PV2) LT 1) then message,$
      'BON map projection requires that PV2_1 keyword be set.'
    pv2_1 = pv2[0]
    if ((PV2_1 lt -90) or (PV2_1 gt 90)) then message,$
      'PV2_1 must satisfy -90 <= PV2_1 <= 90 for BON map projection'