pro satellite_image
; This procedure prompts the user to select a visible MODIS file
; Red, Green, Blue reflectances are put into an array for display
; a smaller image is produced for viewing,
; then the user is prompted to select a geographical subset.
; Adjustments to contrast is performed.
; Color mapping is demonstrated.
; finally, the whole image is mapped using parameters selected above.

;------ ABOUT THE DATASET ----------------
;
; A MODIS 1KM visible/IR set of channels stored in
; HDF-EOS format is used for this demonstration.  
; To read them, the procedures modradread.pro and latlonread.pro
; (written by Brian Vant-Hull) are used.  The data set was
; chosen because it has RGB channels in the visible, but there
; are some idiosyncracies:
; 1. The 1 km resolution visible bands have 1354 pixels in the x direction,
;    while the lat and lon data is at 5 km resolution with 271 pixels in x.
;    The data lines up at the [0,0] corner, so we will trim the last 4 pixels
;    off the 1km and the last 1 pixel off the 5 km resolution data.  To fit on
;    the computer screen, the visible data is first plotted at 5 km resolution.
; 2. The Terra satellite is in descending mode, so that latitude decreases
;    with y.  We will have to invert the order to plot it.

;------ SETUP FOR THIS DATASET ------------
path = '/applications/rsi/'
scale = 5
byte = 255
xdim = 1350
ydim = 2030

device,retain=2,/decomposed
;------------------------------------------
;        GET DATA
;------------------------------------------ 
; select image file
visfile = dialog_pickfile(path=path,filter='*.hdf')

; pick out lat, lon data, trim
modlatlonread,visfile, lat,lon, geodims
lat = lat[0:(geodims[0]-2),*]
lon = lon[0:(geodims[0]-2),*]
; decide if ascending or descending latitude
if (lat[0,0] gt lat[0,geodims[1]-1]) then descending=1 else descending=0

; pick out red, blue, green reflectance
modradread,visfile,red,'EV_250_Aggr1km_RefSB',band=1,/refl,dims
modradread,visfile,grn,'EV_500_Aggr1km_RefSB',band=4,/refl,dims
modradread,visfile,blu,'EV_500_Aggr1km_RefSB',band=3,/refl,dims

;----------------------------------------------
;       INITIAL IMAGE PROCESSING
;----------------------------------------------
; put into RGB array and trim
image = fltarr(3,dims[0],dims[1])
image[0,*,*] = red
image[1,*,*] = grn
image[2,*,*] = blu
; now that it's in the array, save some memory
red = 0 & grn = 0 & blu = 0

; MODIS data only has lat-lon for every 5 pixels
; plus some pixels left over that we want to trim off
image = image[*,0:(xdim-1),0:(ydim-1)]

; adjusting brightness to full scale = 0-1, railing off extremes
raillo = where(image lt 0)
if (raillo[0] ne -1) then image[raillo] = 0
railhi = where(image gt 1)
if (railhi[0] ne -1) then image[railhi] = 1

;--------------------------------------------
;   INITIAL DISPLAY
;-------------------------------------------
; display smaller scale of full image
window,0,xsize=xdim/scale,ysize=ydim/scale
scaleimage = byte*congrid(image,3,xdim/scale,ydim/scale)
TV,scaleimage,/true, order=descending

scaleimage=0 ; save memory

;------------------------------------------
;    LOOK AT LAT and LON for SELECTED POINTS
;------------------------------------------
; display lat/lon by choosing a point with the cursor (while loop)
latlonpt = 'y'
repeat begin
  print,'click point on screen to geolocate'
  cursor,x,y,/down,/device
  if (descending) then y = (geodims[1]-1)-y
  print,'lat=',lat[x,y]
  print,'lon=',lon[x,y]
  print,'another point? (y=yes,n=no)'
  read,latlonpt
endrep until (latlonpt eq 'n')

;-----------------------------------------------
;  NOW REFINE IMAGE DISPLAY OF SUBSETTED SECTION
;-----------------------------------------------
;select and display subsetted image (big while loop)
;each adjustment will repeat until a zero is entered

squareselect = 'y'
repeat begin
  
  print, 'click and drag mouse to select smaller area'
  zoomread=1.0 & contrastread=1.0 & gammaread=1.0
  zoom=1.0 & contrast=1.0 & gamma=1.0
  wset, 0
  cursor,x1,y1,/down,/device
  cursor,x2,y2,/up,/device

  ;draw rectangle on big image
  X = [x1,x2,x2,x1,x1]
  Y = [y1,y1,y2,y2,y1]
  plots, X,Y, /device

 ; expand points to full size scale
   x1 = scale*x1 & x2 = scale*x2
   y1 = scale*y1 & y2 = scale*y2

 ; If necessary fix the y order problem
   if (descending) then begin
      y1 = (ydim-1) - y1
      y2 = (ydim-1) - y2
   endif

  ; make initial plot
   subset_image,image,x1,y1,x2,y2,zoom,contrast,gamma,descending

  ; ZOOM ADJUSTMENT
    while (zoomread ne 0) do begin
       print,'enter zoom factor, zero to continue'
       read,zoomread
       if (zoomread ne 0) then begin
          zoom=zoomread
          subset_image,image,x1,y1,x2,y2,zoom,contrast,gamma,descending
       endif
    endwhile

  ; CONTRAST ADJUSTMENT
    while (contrastread ne 0) do begin
       print,'enter contrast, zero to continue'
       read,contrastread
       if (contrastread ne 0) then begin
          contrast = contrastread
          subset_image,image,x1,y1,x2,y2,zoom,contrast,gamma,descending
       endif
    endwhile

  ; GAMMA ADJUSTMENT
    while (gammaread ne 0) do begin
       print,'enter gamma, zero to continue'
       read,gammaread
       if (gammaread ne 0) then begin
          gamma = gammaread
          subset_image,image,x1,y1,x2,y2,zoom,contrast,gamma,descending
       endif       
    endwhile

   print,'Do you wish to select another area? (y/n)'
   read, squareselect

endrep until (squareselect eq 'n') ; end subsetted image adjustment loop

;------------------------------------------------
; ANNOTATE THE SUBSETTED IMAGE
;------------------------------------------------

; Put a marker at a selected point
wset,1
print,'click on point you want marked'
cursor,x,y,/down,/device
plots,x,y,psym=2,symsize=2,color='0000ff'XL,/device

; print text at selected point
text = ' '
print,'enter your text'
read,text
print,'click where you want the center'
cursor,x,y,/down,/device
xyouts,x,y,text,charsize=2,color='00ff00'XL,alignment=0.5,/device

;------------------------------------------------
; STUDYING A SELECTED REGION OF THE PLOT
;------------------------------------------------

roiloop = 'y'

repeat begin

 ; select a set of points to be studied
  xdimsubset = zoom*(abs(x2-x1)+1)
  ydimsubset = zoom*(abs(y2-y1)+1)
  indices = defroi(xdimsubset,ydimsubset,/nofill)

 ; get unmodified data, first in subset image, then selection
  subimage = image[*,x1:x2,y1:y2]
  subimage = congrid(subimage,3,xdimsubset,ydimsubset)
  subdata = subimage[indices]
 ;NOTE: this won't really work in this case because of the 
 ;      y-ordering issue.  We'd have to convert indices into
 ;      x and y indices, then subtract y from ydim.  Too complicated
 ;      for this little demonstration.

  print,'average value of selected data = ',mean(subdata)

  print,'repeat select region? (y/n)'
  read,roiloop

endrep until (roiloop eq 'n')

;-----------------------------------------------
;  CREATE AND USE A COLORMAP
;-----------------------------------------------
; We will use the subimage data and convert to colors
; We want different ranges of the data to plot to
; different colors.  First we design our colormap.

device,decomposed=0   ; not using RGB vectors

; make color table
;        R    G    B   #
TVlct,   0, 255,   0,  0  ; Green - ground
TVlct,   0,   0, 255,  1  ; Blue -  ambigous
TVlct, 255, 255, 255,  2  ; White - Cloud
TVlct, 180, 180, 180,  3  ; Gray - saturated

; we need to get a one-level version of subimage: pick green
subimage1 = reform(subimage[1,*,*])

; An example of equally spaced intervals for colormap (CM)
; we force levels of data to the integers 0-2,
; since reflectance goes 0-1, multiply and round

subimageCM = floor(3*subimage1)
; so 0-.33 => 0
;    .34-.66 => 1
;    .66-.99 => 2
;      1.0   => 3

; now map it

window,4,xsize=xdimsubset,ysize=ydimsubset
TV,subimageCM, order=descending  ;NOTE: true=0

; now use the colormap with variable ranges
greenrange = [0.0,0.0]
bluerange = [0.0,0.0]
whiterange = [0.0,0.0]
colormap='y'

repeat begin

  print,'enter greenrange low,high'
  read,greenrange
  print,'enter bluerange low,high'
  read,bluerange
  print,'enter whiterange low,high'
  read,whiterange

 ; now put data in ranges to integers
   
   greenindex = where((subimage1 ge greenrange[0]) and (subimage1 lt greenrange[1]))
   if (greenindex[0] ne -1) then subimageCM[greenindex] = 0
   blueindex = where((subimage1 ge bluerange[0]) and (subimage1 lt bluerange[1]))
   if (blueindex[0] ne -1) then subimageCM[blueindex] = 1
   whiteindex = where((subimage1 ge whiterange[0]) and (subimage1 lt whiterange[1]))
   if (whiteindex[0] ne -1) then subimageCM[whiteindex] = 2
   grayindex = where(subimage1 ge whiterange[1])
   if (grayindex[0] ne -1) then subimageCM[grayindex] = 3

   TV,subimageCM,order=descending
   
   print,'would you like to change the ranges? (y/n)'
   read,colormap

endrep until (colormap eq 'n')

device,decomposed=1  ; set back to normal 

;------------------------------------------------
; use final values to make mapped plot
;----------------------------------------------;

; we need to pick a central point for the map
latcenter = lat[geodims[0]/2,geodims[1]/2]
loncenter = lon[geodims[0]/2,geodims[1]/2]

; now use this to set the coordinate system
; cylindrical projection is default, could also set /mercator, etc
window,3
map_set,latcenter,loncenter,limit=[min(lat),min(lon),max(lat),max(lon)], $
    title='My Satellite Map'

; now get lat and lon dimensions to match the image
lat = congrid(lat,xdim,ydim)
lon = congrid(lon,xdim,ydim)

; adjust full image to our selected parameters
image = contrast*image^(1/gamma)
lorail = where(image lt 0)
if (lorail[0] ne -1) then image = image[lorail]
hirail = where(image gt 1)
if (hirail[0] ne -1) then image[hirail] = 1

;unfortunately we can't map RGB, so must cut out all but one color
image = byte*reform(image[1,*,*])
TV,map_patch(image,lon,lat)
map_continents
map_grid,/label

;-------------------------------------------------
;   SAVING IMAGES
;-------------------------------------------------

; save RGB image
rgbsave = ''
filename = ' '
print,'save RGB image? (y/n)'
read,rgbsave
if (rgbsave eq 'y') then begin
   print,'enter a filename'
   print,'.tif extension will be added automatically'
   read,filename
   wset,1
   rgb = tvrd(true=1)
   write_tiff,filename+'.tif',reverse(rgb,3)
   print,'file saved as ',filename+'.tif'
endif

; save map
mapsave=''
filename = ''
print,'save mapped image? (y/n)'
read,mapsave
if (mapsave eq 'y') then begin
   print,'enter a filename'
   print,'.tif extension will be added automatically'
   read,filename
   wset,1
   map = tvrd(true=1)
   write_tiff,filename+'.tif',reverse(map,3)
   print,'file saved as ',filename+'.tif'
endif

end
;***************************************
pro subset_image, image,x1,y1,x2,y2,zoom,contrast,gamma,descending
; creates small image based on cursor clicks

byte = 255

xdimsmall = zoom*(abs(x2-x1)+1)
ydimsmall = zoom*(abs(y2-y1)+1)

imageSubset = image[*,x1:x2,y1:y2]

; apply image processing
imageSubset = contrast*(imageSubset)^(1/gamma)
hirail = where(imageSubset gt 1)
if (hirail[0] ne -1) then imageSubset[hirail] = 1
lorail = where(imagesubset lt 0)
if (lorail[0] ne -1) then imageSubset[lorail] = 0
imageSubset = congrid(imageSubset,3,xdimsmall,ydimsmall)

window,1,xsize=xdimsmall,ysize=ydimsmall
TV,byte*imagesubset,/true,order=descending

end