wcs_demo.pro

basic median filter simulation

  printf,file_unit,"x_len = 2*45.0"
  printf,file_unit,"y_len = 2*45.0"

  printf,file_unit,""
  printf,file_unit,$
  "; Correctly adjust the x and y values for display purposes (they all start "
  printf,file_unit,$
  "; out on the same face)."
  printf,file_unit,"if (face_0[0] ne -1) then begin"
  printf,file_unit,"  x0 = -x(face_0) + 2.d0*x_len"
  printf,file_unit,"  y0 = y(face_0) + y_len"
  printf,file_unit,"  xx(face_0) = x0"
  printf,file_unit,"  yy(face_0) = y0"
  printf,file_unit,"endif"
  printf,file_unit,"if (face_1[0] ne -1) then begin"
  printf,file_unit,"  x1 = -x(face_1) + 2.d0*x_len"
  printf,file_unit,"  y1 = y(face_1)"
  printf,file_unit,"  xx(face_1) = x1"
  printf,file_unit,"  yy(face_1) = y1"
  printf,file_unit,"endif"
  printf,file_unit,"if (face_2[0] ne -1) then begin"
  printf,file_unit,"  x2 = -x(face_2) + x_len"
  printf,file_unit,"  y2 = y(face_2)"
  printf,file_unit,"  xx(face_2) = x2"
  printf,file_unit,"  yy(face_2) = y2"
  printf,file_unit,"endif"
  printf,file_unit,"if (face_3[0] ne -1) then begin"
  printf,file_unit,"  x3 = -x(face_3)"
  printf,file_unit,"  y3 = y(face_3)"
  printf,file_unit,"  xx(face_3) = x3"
  printf,file_unit,"  yy(face_3) = y3"
  printf,file_unit,"endif"
  printf,file_unit,"if (face_4[0] ne -1) then begin"
  printf,file_unit,"  x4 = -x(face_4) - x_len"
  printf,file_unit,"  y4 = y(face_4)"
  printf,file_unit,"  xx(face_4) = x4"
  printf,file_unit,"  yy(face_4) = y4"
  printf,file_unit,"endif"
  printf,file_unit,"if (face_5[0] ne -1) then begin"
  printf,file_unit,"  x5 = -x(face_5) + 2.d0*x_len"
  printf,file_unit,"  y5 = y(face_5) - y_len"
  printf,file_unit,"  xx(face_5) = x5"
  printf,file_unit,"  yy(face_5) = y5"
  printf,file_unit,"endif"

  printf,file_unit,""
  printf,file_unit,$
  "; Define plot ranges by finding which faces are actually used."
  printf,file_unit,"if (face_4[0] ne -1) then x_low = -3*x_len/2 $"
  printf,file_unit,"else if (face_3[0] ne -1) then x_low = -x_len/2 $"
  printf,file_unit,"else if (face_2[0] ne -1) then x_low = x_len/2 $"
  printf,file_unit,$
  "else if ((face_1[0] ne -1) or (face_0[0] ne -1) or (face_5[0] ne -1)) $"
  printf,file_unit,"then x_low = 3*x_len/2"
  printf,file_unit,$
  "if ((face_1[0] ne -1) or (face_5[0] ne -1) or (face_0[0] ne -1)) $"
  printf,file_unit,"  then x_high = 5*x_len/2 $"
  printf,file_unit,"else if (face_2[0] ne -1) then x_high = 3*x_len/2 $"
  printf,file_unit,"else if (face_3[0] ne -1) then x_high = x_len/2 $"
  printf,file_unit,"else if (face_4[0] ne -1) then x_high = -x_len/2"
  printf,file_unit,"if (face_5[0] ne -1) then y_low = -3*y_len/2 $"
  printf,file_unit,$
  "else if ((face_1[0] ne -1) or (face_3[0] ne -1) or (face_2[0] ne -1) or $"
  printf,file_unit,"   (face_4[0] ne -1)) then y_low = -y_len/2 $"
  printf,file_unit,"else if (face_0[0] ne -1) then y_low = y_len/2"
  printf,file_unit,"if (face_0[0] ne -1) then y_high = 3*y_len/2 $"
  printf,file_unit,$
  "else if ((face_1[0] ne -1) or (face_3[0] ne -1) or (face_2[0] ne -1) or $"
  printf,file_unit,"   (face_4[0] ne -1)) then y_high = y_len/2 $"
  printf,file_unit,"else if (face_5[0] ne -1) then y_high = -y_len/2"

  printf,file_unit,""
  printf,file_unit,"; Plot the points calculated by wcssph2xy."
  printf,file_unit,$
  "plot,xx,yy,psym=3,xrange=[x_low,x_high],yrange=[y_low,y_high],xstyle=4,$"
  printf,file_unit,"     ystyle=4"

  printf,file_unit,""
  printf,file_unit,$
  "; Set-up an array with the correct ordering of indices to connect the"
  printf,file_unit,"; latitude lines correctly on faces 1-4."
  printf,file_unit,"face_ind = intarr(1)"
  printf,file_unit,"if (face_4[0] ne -1) then face_ind = [face_ind,face_4]"
  printf,file_unit,"if (face_3[0] ne -1) then face_ind = [face_ind,face_3]"
  printf,file_unit,"if (face_2[0] ne -1) then face_ind = [face_ind,face_2]"
  printf,file_unit,"if (face_1[0] ne -1) then face_ind = [face_ind,face_1]"
  printf,file_unit,"; Draw the latitude lines on faces 1-4"
  printf,file_unit,"if (n_elements(face_ind) gt 1) then begin"
  printf,file_unit,"  face_ind = face_ind(1:*)"
  printf,file_unit,"  xxx = xx(face_ind)"
  printf,file_unit,"  yyy = yy(face_ind)"
  printf,file_unit,"  for i = 0,num_lat - 1 do begin"
  printf,file_unit,"    index = where(latitude(face_ind) eq latitude(0,i))"
  printf,file_unit,"    if (index[0] ne -1) then begin"
  printf,file_unit,"      tempx = xxx(index)"
  printf,file_unit,"      tempy = yyy(index)"
  printf,file_unit,"      index = sort(tempx)"
  printf,file_unit,$
  "         if (((360 - abs(longitude(0,0) - longitude(num_lon - 1,0))) le $"
  printf,file_unit,$
  "          lon_spacing) or (max(longitude(index)) le 135) or $"
  printf,file_unit,$
"       (min(longitude(index)) gt 135)) then oplot,tempx(index),tempy(index) $"
  printf,file_unit,"        else begin"
  printf,file_unit,"           lon_ind = 0"
  printf,file_unit,$
 "           repeat lon_ind=lon_ind+1 until (longitude(index(lon_ind)) gt 135)"
  printf,file_unit,"          index_1 = index(0:lon_ind - 1)"
  printf,file_unit,"          index_2 = index(lon_ind:*)
  printf,file_unit,"          oplot,tempx(index_1),tempy(index_1)"
  printf,file_unit,"          oplot,tempx(index_2),tempy(index_2)"
  printf,file_unit,"        endelse"
  printf,file_unit,"      endif"
  printf,file_unit,"    endfor"
  printf,file_unit,"  endif"
  printf,file_unit,""
  printf,file_unit,"; Draw latitude lines on faces 0 and 5"
  printf,file_unit,"  for i = 0,num_lat - 1 do begin"
  printf,file_unit,"    if (face_0[0] ne -1) then begin"
  printf,file_unit,"      index = where(latitude(face_0) eq latitude(0,i))"
  printf,file_unit,"      if (index[0] ne -1) then begin"
  printf,file_unit,$
  "        if ((360 - abs(longitude(0,0) - longitude(n_elements(x) - 1))) $"
  printf,file_unit,"                                    le lon_spacing) then $"
  printf,file_unit,$
  "          oplot,[x0(index),x0(index[0])],[y0(index),y0(index[0])] $"
  printf,file_unit,"        else oplot,x0(index),y0(index)"
  printf,file_unit,"      endif"
  printf,file_unit,"    endif"
  printf,file_unit,"    if (face_5[0] ne -1) then begin"
  printf,file_unit,"      index = where(latitude(face_5) eq latitude(0,i))"
  printf,file_unit,"      if (index[0] ne -1) then begin"
  printf,file_unit,$
  "        if ((360 - abs(longitude(0,0) - longitude(n_elements(x) - 1))) $"
  printf,file_unit,"                                    le lon_spacing) then $"
  printf,file_unit,$
  "          oplot,[x5(index),x5(index[0])],[y5(index),y5(index[0])] $"
  printf,file_unit,"        else oplot,x5(index),y5(index)" 
  printf,file_unit,"      endif"
  printf,file_unit,"    endif"
  printf,file_unit,"  endfor"
  printf,file_unit,""
  printf,file_unit,"; Draw boxes around each face and draw longitude lines"
  printf,file_unit,"  for i = 0,num_lon - 1 do begin"
  printf,file_unit,"    if (face_4[0] ne -1) then begin"
  printf,file_unit,"      index = where(longitude(face_4) eq longitude(i,0))"
  printf,file_unit,"      if (index[0] ne -1) then oplot,x4(index),y4(index)"
  printf,file_unit,"      plots,[-3*x_len/2,-x_len/2],[-y_len/2,-y_len/2]"
  printf,file_unit,"      plots,[-3*x_len/2,-x_len/2],[y_len/2,y_len/2]"
  printf,file_unit,"      plots,[-x_len/2,-x_len/2],[-y_len/2,y_len/2]"
  printf,file_unit,"      plots,[-3*x_len/2,-3*x_len/2],[-y_len/2,y_len/2]"
  printf,file_unit,"    endif"
  printf,file_unit,"    if (face_2[0] ne -1) then begin"
  printf,file_unit,"      index = where(longitude(face_2) eq longitude(i,0))"
  printf,file_unit,"      if (index[0] ne -1) then oplot,x2(index),y2(index)"
  printf,file_unit,"      plots,[x_len/2,3*x_len/2],[-y_len/2,-y_len/2]"
  printf,file_unit,"      plots,[x_len/2,3*x_len/2],[y_len/2,y_len/2]"
  printf,file_unit,"      plots,[x_len/2,x_len/2],[-y_len/2,y_len/2]"
  printf,file_unit,"      plots,[3*x_len/2,3*x_len/2],[-y_len/2,y_len/2]"
  printf,file_unit,"    endif"
  printf,file_unit,"    if (face_3[0] ne -1) then begin"
  printf,file_unit,"      index = where(longitude(face_3) eq longitude(i,0))"
  printf,file_unit,"      if (index[0] ne -1) then oplot,x3(index),y3(index)"
  printf,file_unit,"      plots,[-x_len/2,x_len/2],[-y_len/2,-y_len/2]"
  printf,file_unit,"      plots,[-x_len/2,x_len/2],[y_len/2,y_len/2]"
  printf,file_unit,"      plots,[-x_len/2,-x_len/2],[-y_len/2,y_len/2]"
  printf,file_unit,"      plots,[x_len/2,x_len/2],[-y_len/2,y_len/2]"
  printf,file_unit,"    endif"
  printf,file_unit,"    if (face_1[0] ne -1) then begin"
  printf,file_unit,"      index = where(longitude(face_1) eq longitude(i,0))"
  printf,file_unit,"      if (index[0] ne -1) then oplot,x1(index),y1(index)"
  printf,file_unit,"      plots,[3*x_len/2,5*x_len/2],[-y_len/2,-y_len/2]"
  printf,file_unit,"      plots,[3*x_len/2,5*x_len/2],[y_len/2,y_len/2]"
  printf,file_unit,"      plots,[3*x_len/2,3*x_len/2],[-y_len/2,y_len/2]"
  printf,file_unit,"      plots,[5*x_len/2,5*x_len/2],[-y_len/2,y_len/2]"
  printf,file_unit,"    endif"
  printf,file_unit,"    if (face_0[0] ne -1) then begin"
  printf,file_unit,"      index = where(longitude(face_0) eq longitude(i,0))"
  printf,file_unit,"      if (index[0] ne -1) then oplot,x0(index),y0(index)"
  printf,file_unit,"      plots,[3*x_len/2,5*x_len/2],[y_len/2,y_len/2]"
  printf,file_unit,"      plots,[3*x_len/2,5*x_len/2],[3*y_len/2,3*y_len/2]"
  printf,file_unit,"      plots,[3*x_len/2,3*x_len/2],[y_len/2,3*y_len/2]"
  printf,file_unit,"      plots,[5*x_len/2,5*x_len/2],[y_len/2,3*y_len/2]"
  printf,file_unit,"    endif"
  printf,file_unit,"    if (face_5[0] ne -1) then begin"
  printf,file_unit,"      index = where(longitude(face_5) eq longitude(i,0))"
  printf,file_unit,"      if (index[0] ne -1) then oplot,x5(index),y5(index)"
  printf,file_unit,"      plots,[3*x_len/2,5*x_len/2],[-3*y_len/2,-3*y_len/2]"
  printf,file_unit,"      plots,[3*x_len/2,5*x_len/2],[-y_len/2,-y_len/2]"
  printf,file_unit,"      plots,[3*x_len/2,3*x_len/2],[-3*y_len/2,-y_len/2]"
  printf,file_unit,"      plots,[5*x_len/2,5*x_len/2],[-3*y_len/2,-y_len/2]"
  printf,file_unit,"    endif"
  printf,file_unit,"  endfor"
  printf,file_unit,""
  printf,file_unit,";LABELS"
  printf,file_unit,"  if (lat_index[0] ne -1) then $"
  printf,file_unit,"  xyouts,xx(0,lat_index),yy(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]),yy(lon_index,index[0]),orientation=90,$"
  printf,file_unit,$
"       strcompress(string(long(longitude(lon_index,index[0])))),alignment=0.5"
endelse
end

; PROCEDURE FOR OPTION 2
pro inversion_error,file_unit,map,param1,param2
printf,file_unit,";CONVERSION"
printf,file_unit,$
"; Convert the x-y coordinates into spherical coordinates by using wcsxy2sph."
if (map lt 23) then begin
  if (n_elements(param1) eq 0) then begin
    printf,file_unit,"wcsxy2sph,x,y,longitude_inv,latitude_inv,map"
  endif else if (n_elements(param2) eq 0) then begin
  printf,file_unit,"wcsxy2sph,x,y,longitude_inv,latitude_inv,map,pv2=param1"
  endif else begin
    printf,file_unit,$
     "wcsxy2sph,x,y,longitude_inv,latitude_inv,map,pv2= [param1, param2] "
  endelse
endif else begin
  printf,file_unit,$
 "; The variable face must be declared with the same structure as latitude and"
  printf,file_unit,"; longitude before calling wcsxy2sph." 
  printf,file_unit,"wcsxy2sph,x,y,longitude_inv,latitude_inv,map,face=face"
endelse

printf,file_unit,""
printf,file_unit,";PLOTTING"
printf,file_unit,"; Plot the resulting map."
printf,file_unit,"lon_delta = (max(longitude_inv) - min(longitude_inv))/20"
printf,file_unit,"lat_delta = (max(latitude_inv) - min(latitude_inv))/20"
printf,file_unit,$
    "plot,longitude_inv,latitude_inv,psym = 3,xrange = [min(longitude_inv) - $"
printf,file_unit,$
"    lon_delta,max(longitude_inv) + lon_delta],yrange = [min(latitude_inv) - $"
printf,file_unit,$
"    lat_delta,max(latitude_inv) + lat_delta],xstyle = 4,ystyle = 4"
printf,file_unit,"; Draw lines connecting equal longitudes"
printf,file_unit,$
       "for i = 0,num_lon - 1 do oplot,longitude_inv[i,*],latitude_inv[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,longitude_inv[*,i],latitude_inv[*,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,[longitude_inv(index,i),$"
printf,file_unit,$
       "      longitude_inv(0:index[0]-1,i)],[latitude_inv(index,i),$"
printf,file_unit,"      latitude_inv(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,longitude_inv(0:index[0] - 1,i),latitude_inv(0:index[0] - 1,i)"
printf,file_unit,"      oplot,longitude_inv(index,i),latitude_inv(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,longitude_inv(lon_index,0),latitude_inv(lon_index,0) - lat_delta,$"
printf,file_unit,$
    "       orientation=90,strcompress(string(long(longitude(lon_index,0)))),$"
printf,file_unit,"       alignment=0.5"
printf,file_unit,"lat1_index = where(longitude[0,*] eq max(longitude[0,*]))"
printf,file_unit,"if (lat_index[0] ne -1) then $"
printf,file_unit,$
"xyouts,max(longitude_inv) + lon_delta,latitude_inv(lat1_index[0],lat_index),$"
printf,file_unit,$
"       alignment=0.5,strcompress(string(long(latitude(0,lat_index))))"

printf,file_unit,"read,'Press return to continue',key"
print,"  In order to make the scripts wcssph2xy.pro and wcsxy2sph.pro"
print,"invertible and minimize the error in the process, it was necessary to"
print,"offset the latitude of all points at the poles by a small amount."
print,"When viewing the difference between the original longitude and"
print,"latitude and the longitude and latitude after points are run through"
print,"wcssph2xy.pro and wcsxy2sph.pro, the offset at the poles will show up"
print,"as vertical lines.  This overshadows any numerical error elsewhere"
print,"by orders of magnitude.  The default is to ignore these errors, but"
print,"to include them, enter n at the prompt"
print,""
key = ""
repeat $
  read,"Ignore offset at poles when plotting vector field (y or n):",key $
until ((key eq "y") or (key eq "n")) 

if (key eq "y") then begin
  printf,file_unit,"poles = where(abs(abs(latitude_inv) - 9.d1) le 573.d-4)"