wcs_demo.pro

basic median filter simulation

;+
; NAME:
;       WCS_DEMO
;
; PURPOSE:
;       Demonstrate the basic capabilities of procedures WCSSPH2XY & WCSXY2SPH
;
; CATEGORY:
;       Mapping and Auxilary FITS Demo Routine
;
; CALLING SEQUENCE:
;
;       .run wcs_demo: compiles wcs_demo and the supporting demo routines
;       wcs_demo: run the demo
;
; INPUT PARAMETERS:
;
;       none
;
; OUTPUT PARAMETERS:
;       none
;
; PROCEDURE:
;
;       This is a demo program which is meant to call the routines 
;       wcssph2xy.pro and wcsxy2sph.pro.  Since the purpose of this
;       routine is both to show what the routines can do and what the
;       user has to do, a file is created with all of the commands 
;       needed to complete the desired operation.  Wcs_demo actually 
;       executes this command file, so the user can exactly duplicate
;       the results by simply re-executing this file.  Also, this 
;       allows a user to edit an already existing file which calls 
;       wcssph2xy.pro and wcsxy2sph.pro properly and extend the file's
;       usefulness.  This demo program allows several possible tests.
;       The first option is to simply draw a grid of evenly spaced
;       latitude and longitude lines in a particular map transformation.
;       Another possibility is to do a full loop, creating a Cartesian
;       grid of latitude and longitude lines and calling wcssph2xy.pro
;       to convert them to a particular map.  Then, wcsxy2sph.pro is
;       called to invert the process and the difference between the
;       original and final latitudes and longitudes can be plotted.
;       This allows one to assess the level of the numerical errors
;       introduced by the mapping routines.  A third possible option is to
;       look at some of the map transformations and include rotations of
;       the reference points so that a different perspective is given.
;
; COMMON BLOCKS:
;       none
;
; PROCEDURES CALLED:
;       SPHDIST(), WCSXY2SPH, WCSSPH2XY
; COPYRIGHT NOTICE:
;
;       Copyright 1991, 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
;       1.2     3/19/96 - J. Bloch - LANL
;                        - Made compatible with wcslib-2.2 by Calabretta.
;       Converted to IDL V5.0   W. Landsman   September 1997
;       Updated for conical projections W. Landsman  July 2003
;-

; PROCEDURE FOR OPTION 1
pro wcssph2xy_plot,file_unit,map,param1,param2
printf,file_unit,";PLOTTING"
printf,file_unit,"; Plot the resulting map."
if ((map ge 0) and (map le 22)) then begin
  ; For all but the spherical cube projections, simply plot the results from
  ; wcssph2xy.pro as is.
  printf,file_unit,"xdelta = (max(xx) - min(xx))/20"
  printf,file_unit,"ydelta = (max(y) - min(y))/20"
  printf,file_unit,$
         "plot,xx,y,psym = 3,xrange = [min(xx) - xdelta,max(xx) + xdelta],$"
  printf,file_unit,$
         "yrange = [min(y) - ydelta,max(y) + ydelta],xstyle = 4,ystyle = 4"

  ; ZENITHAL PROJECTIONS.
  if ((map ge 1) and (map le 8)) then begin

    printf,file_unit,""
    printf,file_unit,$
    "; Only connect latitude lines in a full circle if the longitude"
    printf,file_unit,"; values cover the full circle."
    printf,file_unit,$
      "if (360 - abs(longitude(0,0) - longitude(n_elements(xx[*,0])-1)) $"
    printf,file_unit,"                             le lon_spacing) $"
    printf,file_unit,$
      "then for i = 0,num_lat - 1 do oplot,[xx[*,i],xx(0,i)],[y[*,i],y(0,i)] $"
    printf,file_unit,"else for i = 0,num_lat - 1 do oplot,xx[*,i],y[*,i]"

    printf,file_unit,""
    printf,file_unit,$
    "; Connect the longitude lines from the poles outward."
    printf,file_unit,"for i = 0,num_lon - 1 do oplot,xx[i,*],y[i,*]"

    printf,file_unit,""
    printf,file_unit,";LABELS"
    printf,file_unit,$
    "; Label the latitude and longitude lines and correctly orient the labels."
    printf,file_unit,"j = 0"
    printf,file_unit,"repeat begin"
    printf,file_unit,"  i = lon_index(j)"
    printf,file_unit,"  xyouts,xx(i,0)-xdelta*sin(longitude(i,0)/!radeg),$"
    printf,file_unit,"         y(i,0)-ydelta*cos(longitude(i,0)/!radeg),$"
    printf,file_unit,$
          "         strcompress(string(long(longitude(i,0)))),alignment=0.5,$"
    printf,file_unit,"         orientation=360-longitude(i,0)"
    printf,file_unit,"  j = j + 1"
    printf,file_unit,"endrep until (j eq n_elements(lon_index))"
    printf,file_unit,"if (lat_index[0] ne -1) then $"
    printf,file_unit,"  xyouts,xx(0,lat_index),y(0,lat_index),$"
    printf,file_unit,"       strcompress(string(long(latitude(0,lat_index))))"

  ; CYLINDRICAL PROJECTIONS
  endif else if (((map ge 9) and (map le 12)) or (map eq 0)) then begin
    printf,file_unit,""
    printf,file_unit,"; Draw lines connecting equal longitudes"
    printf,file_unit,"for i = 0,num_lon - 1 do oplot,xx[i,*],y[i,*]"
    printf,file_unit,"; Draw lines connecting equal latitudes"
    printf,file_unit,$
    "if ((min(longitude[*,0]) ge 180) or (max(longitude[*,0]) lt 180)) then $"
    printf,file_unit,"  for i = 0,num_lat - 1 do oplot,xx[*,i],y[*,i] $"
    printf,file_unit,"else begin"
    printf,file_unit,"  index = where(longitude[*,0] ge 180)"
    printf,file_unit,$
    "  if ((360 - max(longitude[*,0]) + min(longitude[*,0])) le lon_spacing) $"
    printf,file_unit,"    then begin"
    printf,file_unit,$
    "    for i = 0, num_lat - 1 do oplot,[xx(index,i),xx(0:index[0]-1,i)],$"
    printf,file_unit,$
    "                                    [y(index,i),y(0:index[0]-1,i)]"
    printf,file_unit,"  endif else begin"
    printf,file_unit,"    for i = 0,num_lat - 1 do begin"
    printf,file_unit,"      oplot,xx(0:index[0] - 1,i),y(0:index[0] - 1,i)"
    printf,file_unit,"      oplot,xx(index,i),y(index,i)"
    printf,file_unit,"    endfor"
    printf,file_unit,"  endelse"
    printf,file_unit,"endelse"

    printf,file_unit,""
    printf,file_unit,";LABELS"
    printf,file_unit,$
    "; Label the latitude and longitude lines and correctly orient the labels."
    printf,file_unit,$
    "xyouts,xx(lon_index,0),y(lon_index,0) - ydelta,orientation=90,$"
    printf,file_unit,$
    "       strcompress(string(long(longitude(lon_index,0)))),alignment=0.5"
    printf,file_unit,"y_index = where(longitude[0,*] eq max(longitude[0,*]))"
    printf,file_unit,"if (lat_index[0] ne -1) then $"
    printf,file_unit,$
    "xyouts,max(xx) + xdelta,y(y_index[0],lat_index),alignment=0.5,$"
    printf,file_unit,"       strcompress(string(long(latitude(0,lat_index))))"

  ; CONICAL PROJECTIONS
  endif else if ((map ge 13) and (map le 16)) then begin
    printf,file_unit,""
    printf,file_unit,"; Draw lines of longitude out from the poles."
    printf,file_unit,"for i = 0,num_lon - 1 do oplot,xx[i,*],y[i,*]"

    printf,file_unit,$
    "; Draw lines of latitude, making sure to break the line at 180 degrees."
    printf,file_unit,"index = where(longitude[*,0] ge 180)"
    printf,file_unit,"if (index[0] ne -1) then $"
    printf,file_unit,$
          "  for i = 0,num_lat - 1 do oplot,[xx(index,i),xx(0:index[0]-1,i)],$"
    printf,file_unit,"                        [y(index,i),y(0:index[0]-1,i)] $"
    printf,file_unit,"else begin"
    printf,file_unit,"  index = where(longitude[*,0] eq max(longitude[*,0]))"
    printf,file_unit,$
    "  for i = 0,num_lat - 1 do oplot,xx(0:index[0],i),y(0:index[0],i)"
    printf,file_unit,"endelse"

    printf,file_unit,""
    printf,file_unit,";LABELS"
    printf,file_unit,$
    "; Label latitude and longitude and correctly orient the labels."
    printf,file_unit,"j = 0"
    printf,file_unit,"if (min(longitude) lt 180) then begin"
    printf,file_unit,$
    "  lon_ind_1 = lon_index(where(longitude(lon_index,0) lt 180))"
    printf,file_unit,$
    "  lon_ind_1 = lon_ind_1(reverse(sort(longitude(lon_ind_1,0))))"
    printf,file_unit,"endif"
    printf,file_unit,"if (max(longitude) ge 180) then begin"
    printf,file_unit,$
    "  lon_ind_2 = lon_index(where(longitude(lon_index,0) ge 180))"
    printf,file_unit,$
    "  lon_ind_2 = lon_ind_2(reverse(sort(longitude(lon_ind_2,0))))"
    printf,file_unit,"endif"
    printf,file_unit,$
    "if ((n_elements(lon_ind_1) ne 0) and (n_elements(lon_ind_2) ne 0)) then $"
    printf,file_unit,"  lon_index = [lon_ind_1,lon_ind_2] $"
    printf,file_unit,"else if (n_elements(lon_ind_1) ne 0) then $"
    printf,file_unit,"  lon_index = lon_ind_1 $"
    printf,file_unit,"else if (n_elements(lon_ind_2) ne 0) then $"
    printf,file_unit,"  lon_index = lon_ind_2"
    if (param2 gt -param1) then begin
    printf,file_unit,"repeat begin"
    printf,file_unit,"  i = lon_index(j)"
    printf,file_unit,"  i1 = lon_index(j + 1)"
    printf,file_unit,"  angle = atan(y(i1,0) - y(i,0),xx(i1,0) - xx(i,0))"
    printf,file_unit,$
    "  xyouts,xx(i,0) + xdelta*sin(angle),y(i,0) - ydelta*cos(angle),$"
    printf,file_unit,$
    "         strcompress(string(long(longitude(i,0)))),alignment = 0.5,$"
    printf,file_unit,"         orientation = !radeg*angle"
    printf,file_unit,"  j = j + 1"
    printf,file_unit,"endrep until (j eq (n_elements(lon_index) - 1))"
    endif else begin
    printf,file_unit,"end_index = n_elements(xx[i,*]) - 1"
    printf,file_unit,"repeat begin"
    printf,file_unit,"  i = lon_index(j)"
    printf,file_unit,"  i1 = lon_index(j + 1)"
    printf,file_unit,"  angle = atan(y(i1,end_index) - y(i,end_index),$"
    printf,file_unit,"               xx(i1,end_index) - xx(i,end_index))"
    printf,file_unit,$
    "  xyouts,xx(i,end_index) - xdelta*sin(angle),y(i,end_index) + $"
    printf,file_unit,$
    "         ydelta*cos(angle),strcompress(string(long(longitude($"
    printf,file_unit,"i,end_index)))),alignment=0.5,orientation=!radeg*angle"
    printf,file_unit,"  j = j + 1"
    printf,file_unit,"endrep until (j eq n_elements(lon_index) - 1)"
    endelse
    printf,file_unit,$
    "if (lat_index[0] ne -1) then xyouts,xx(0,lat_index),y(0,lat_index),$"
    printf,file_unit,$
    "                        strcompress(string(long(latitude(0,lat_index))))"

  ; CONVENTIONAL PROJECTIONS
  endif else if ((map ge 17) and (map le 22)) then begin
    printf,file_unit,""
    printf,file_unit,"; Draw lines of longitude"
    printf,file_unit,"for i = 0,num_lon - 1 do oplot,xx[i,*],y[i,*]"

    printf,file_unit,$
    "; Draw lines of latitude, breaking the line at 180 degrees."
    printf,file_unit,$
    "if ((min(longitude[*,0]) ge 180) or (max(longitude[*,0]) lt 180)) then $"
    printf,file_unit,"  for i = 0,num_lat - 1 do oplot,xx[*,i],y[*,i] $"
    printf,file_unit,"else begin"
    printf,file_unit,"  index = where(longitude[*,0] ge 180)"
    printf,file_unit,$
    "  if ((360 - max(longitude[*,0]) + min(longitude[*,0])) le lon_spacing) $"
    printf,file_unit,"    then begin"
    printf,file_unit,$
    "    for i = 0, num_lat - 1 do oplot,[xx(index,i),xx(0:index[0]-1,i)],$"
    printf,file_unit,$
    "                                    [y(index,i),y(0:index[0]-1,i)]"
    printf,file_unit,"  endif else begin"
    printf,file_unit,"    for i = 0,num_lat - 1 do begin"
    printf,file_unit,"      oplot,xx(0:index[0] - 1,i),y(0:index[0] - 1,i)"
    printf,file_unit,"      oplot,xx(index,i),y(index,i)"
    printf,file_unit,"    endfor"
    printf,file_unit,"  endelse"
    printf,file_unit,"endelse"

    printf,file_unit,""
    printf,file_unit,";LABELS"
    printf,file_unit,$
    "; Label latitude and longitude lines and orient the labels correctly."
    printf,file_unit,"if (lat_index[0] ne -1) then $"
    printf,file_unit,"xyouts,xx(0,lat_index),y(0,lat_index),$"
    printf,file_unit,"       strcompress(string(long(latitude(0,lat_index))))"
    printf,file_unit,$
    "index = where(abs(latitude[0,*]) eq min(abs(latitude[0,*])))"
    printf,file_unit,$
    "xyouts,xx(lon_index,index[0]),y(lon_index,index[0]),orientation=90,$"
    printf,file_unit,$
"       strcompress(string(long(longitude(lon_index,index[0])))),alignment=0.5"
  endif

; SPHERICAL CUBE PROJECTIONS
endif else begin
  printf,file_unit,"xx = -x"
  printf,file_unit,"yy = y"

  printf,file_unit,""
  printf,file_unit,"; Make arrays with the locations of all points."
  printf,file_unit,"face_0 = where(face eq 0)"
  printf,file_unit,"face_1 = where(face eq 1)"
  printf,file_unit,"face_2 = where(face eq 2)"
  printf,file_unit,"face_3 = where(face eq 3)"
  printf,file_unit,"face_4 = where(face eq 4)"
  printf,file_unit,"face_5 = where(face eq 5)"

  printf,file_unit,""
  printf,file_unit,"; Define the size of the box around each face."