pro misr_envi,event

;  This routine makes data from a MISR Level 1B2 radiance file
;  in the original stacked block grid available in ENVI.
;
;  It allows the user to specify a MISR file to work with and to
;  select contiguous blocks of interest within the file.
;
;  For each of the four bands, the blocks of interest are appended
;  together, applying the horizontal block offsets in the process.
;
;  The block subset is geolocated by providing the appropriate
;  projection information calculated for the user's start block.
;
;  The bands are added to the ENVI Available Bands List.
;
;  If the AGP file corresponding to the data file's path is available,
;  the latitude, longitiude, and elevation fields are read, resized if
;  for the nadir camera, geolocated, and added to the Available Bands
;  List.
;
;  These data are created in ENVI in memory. The user can then choose
;  to save them to a file by right-clicking on the data set heading
;  in the Available Bands List and selecting "Save Selected File to
;  Disk..."
;
;  External routines:
;      dms_to_ddeg.pro - converts lat/lon from HDF-EOS DDDMMMSSS.SS format
;                        to decimal degrees
;      get_offsets.pro - reads the block offset values from a MISR
;                        file
;      get_user_blocks.pro - GUI/event handler routine pair to get user's
;                        desired block range
;      read_global_attr.pro - gets an HDF global Ifile) attribute value
;      read_grid_attr.pro - gets an HDF-EOS grid attribute value
;      som_block_corner_meters.pro - returns the corner values in meters
;                        for a block in a MISR SOM stacked-block grid file
;
;  Modification history:
;
;      April 2004 - original version
;
;  Contact information:
;      email:  larc@eos.nasa.gov
;      phone:  757-864-8656
;      mail:   User and Data Services
;              NASA Langley Atmospheric Sciences Data Center
;              MS 157D, 2 S. Wright St.
;              Hampton, VA  23681-2199
;
;
;------------------------------------------------------------------------

COMPILE_OPT STRICTARR

;  edit the values below to specify the locations of the MISR Level 1B2
;  data files and the AGP files

DATApath = "G:\data\misr\Level1B2"
AGPpath = "G:\data\misr\AGP"

;  these are the band names that will appear in the ENVI Available
;  Bands List--edit them as desired

bnames = ['Blue Radiance','Green Radiance','Red Radiance','NIR Radiance']

radfilepath = DIALOG_PICKFILE(PATH=DATApath,FILTER="*.hdf",TITLE="Select MISR L1B2 File")
if (radfilepath eq '') then return

radfile = FILE_BASENAME(radfilepath)

;  note--extracting the path label in this manner depends on the MISR
;  L1B2 files retaining their original names; if the files have been renamed
;  this code may need to be modified

fparts = strsplit(radfile,'_',/EXTRACT)
path = fparts[5]

num_grid = EOS_GD_INQGRID(radfilepath,gridlist)
grids = strsplit(gridlist,',',/EXTRACT)
fid = EOS_GD_OPEN(radfilepath,/READ)

;  get dimension information and projection parameters
;  use the blue grid so that appropriate information will be a
;    obtained for both nadir and off-nadir camera files

gid = EOS_GD_ATTACH(fid,grids[0])
status = EOS_GD_GRIDINFO(gid,xdimsize,ydimsize,upleft,lowright)
status = EOS_GD_PROJINFO(gid,projcode,zonecode,spherecode,projparam)
status = EOS_GD_DETACH(gid)

case ydimsize of
512: pixsize = 1100
2048: pixsize = 275
else :
endcase

;  get the start and end blocks for the file

status = read_global_attr(radfilepath,"Start_block",start_block)
status = read_global_attr(radfilepath,"End block",end_block)
start_block = start_block[0]
end_block = end_block[0]


;  get the block offsets to properly calculate the xsize of the
;  final image and the positioning of the blocks within the
;  aggregated image data array

offsets = get_offsets(radfilepath,grids[0])

;  get the user's desired block range
;  pass in the file start and end blocks, get back user's start
;    and end blocks

get_user_blocks,start_block,end_block
num_blocks = end_block - start_block + 1

WIDGET_CONTROL,/HOURGLASS

;  make array big enough to hold the specified blocks including offsets

pos = lonarr(num_blocks)
pos[0] = 0
for i=1,num_blocks-1 do pos[i] = pos[i-1] + long(offsets[i+start_block-2])
maxpos = max(pos,min=minpos)
offarray = abs(minpos) + pos

total_lines = num_blocks * xdimsize

banddata = fltarr(maxpos-minpos+ydimsize,total_lines,4)

;  read data for user's blocks for the four wavelengths
;  if it is a non-nadir camera, reduce red band resolution to match other bands
;  later may want to look at using sharpening to use red band at full res
;    to improve lower res bands

start = [0,0,start_block-1]

for i=0,3 do begin

  ;  read the data (there is one grid for each band, and that
  ;  grid contains one field with the radiance data)

  gid = EOS_GD_ATTACH(fid,grids[i])
  status = EOS_GD_INQFIELDS(gid,fieldlist,ranks,fieldtypes)
  fields = strsplit(fieldlist,',',/EXTRACT)
  status = EOS_GD_FIELDINFO(gid,fields[0],rank,dims,numtype,dimlist)

  edge = [dims[0],dims[1],num_blocks]
  status = EOS_GD_READFIELD(gid,fields[0],bdata,START=start,EDGE=edge)

;  convert from raw values to scaled values, removing quality indicator
;  in low order 2 bits and multiplying by the scale factor, use zeros for
;  bad data values since this is for imaging purposes

  bdata = bdata/4
  notgood = where(bdata ge 16378,count)
  if (count gt 0) then bdata [notgood] = 0

  status = EOS_GD_DETACH(gid)
  status = read_grid_attr(radfilepath,grids[i],"Scale factor",scale)

  bdata = float(bdata) * scale[0]

;  here's where the red band data for an off-nadir camera is redimensioned
;  down to the resolution of the other three bands

  if (dims[0] gt ydimsize) then bdata = CONGRID(bdata,ydimsize,xdimsize,num_blocks,/CENTER)

  dims = size(bdata,/dimensions)
  if (n_elements(dims) eq 3) then bdata = reform(bdata,dims[0],dims[1]*dims[2],/overwrite)

  ;  organize band in a temporary array using the offset values to properly position the
  ;  blocks in the array
  ;  (note: recode to eliminate temporary array to reduce memory needed)

  t = fltarr(maxpos-minpos+ydimsize,total_lines)
  stline = 0
  endline = xdimsize - 1

  t[0,0] = bdata[*,*]

  for j=0,num_blocks -1 do begin
    if (offarray[j] ne 0) then begin
      t[0,stline] =  shift(t[*,stline:endline],offarray[j],0)
    endif
    stline = endline + 1
    endline = endline + xdimsize
  endfor

  banddata[0,0,i] = t

endfor   ; loop through the four bands

status = EOS_GD_CLOSE(fid)

;  set up map info for this subregion with ENVI_MAP_INFO_CREATE

projname = "MISR SOM PATH "+ strtrim(string(path),2) + " BLOCK " + $
             strtrim(string(start_block),2)
descrip = FILE_BASENAME(radfilepath,'.hdf')+"_B"+strtrim(string(start_block),2)+"-"+$
            strtrim(string(end_block),2)

som_block_corner_meters,radfilepath,start_block,somupleft,somlowright

;  projection parameters needed for the generic SOM A projection:
;  a, b, IncAng, AscLong, x0, y0, PSRev, LRat

somprojparm = [6378137.0, 6356752.3, 98.303819, dms_to_ddeg(projparam[4]), $
               0.0, 0.0, 98.880000, 0.000000]

map_info = ENVI_MAP_INFO_CREATE(DATUM="WGS-84",PS=[pixsize,pixsize],$
           ROTATION=90.0,NAME=projname,$
           MC=[offarray[0],0,somupleft[0],somupleft[1]],PARAMS=somprojparm,TYPE=37)

;  use ENVI_ENTER_DATA to add data to the Available Bands List
;  note--change default stretch as desired below

ENVI_ENTER_DATA,banddata,BNAMES=bnames,DEF_STRETCH=ENVI_DEFAULT_STRETCH_CREATE(/GAUSSIAN),$
  MAP_INFO=map_info,DESCRIP=descrip

; use corresponding AGP file if available

agpfilter = "MISR_AM1_AGP_" + path + "*.hdf"
agpfile = DIALOG_PICKFILE(PATH=AGPpath,FILTER=agpfilter,TITLE="Select AGP File")
WIDGET_CONTROL,/HOURGLASS

if agpfile ne '' then begin

;    extract lat, lon and elevation for the user's specifed blocks
;    concatenate the selected blocks, geolocate, and add to the available
;      bands list
;    if radiance data was from the nadir camera, interpolate the AGP data
;      up to that resolution (AGP data is on a 1.1 km or 512x128 pixel grid)

  agpfid = EOS_GD_OPEN(agpfile,/READ)
  num_agpgrid = EOS_GD_INQGRID(agpfile,agpgridlist)
  agpgrids = strsplit(agpgridlist,',',/EXTRACT)
  agpgid = EOS_GD_ATTACH(agpfid,agpgrids[0])
  agpfields = ["AveSceneElev","GeoLatitude","GeoLongitude"]

  agpdata = fltarr(maxpos-minpos+ydimsize,total_lines,3)

  for f=0,n_elements(agpfields)-1 do begin
    status = EOS_GD_FIELDINFO(agpgid,agpfields[f],rank,dims,numtype,dimlist)
    start = [0,0,start_block-1]
    edge = [dims[0],dims[1],num_blocks]
    status = EOS_GD_READFIELD(gid,agpfields[f],fdata,START=start,EDGE=edge)

    ;  regrid data if associated radiance data is from a nadir camera

      if (dims[0] lt ydimsize) then fdata = CONGRID(fdata,ydimsize,xdimsize,num_blocks,/CENTER)

    ;  reform data from three dimensions to two--append blocks into second dimension
      dims = size(fdata,/dimensions)
      if (n_elements(dims) eq 3) then fdata = REFORM(fdata,dims[0],dims[1]*dims[2],/OVERWRITE)

      ;  organize data in a temporary array using the offset values to properly position the
      ;  blocks in the array
      ;  (note: recode to eliminate temporary array to reduce memory needed)

      t = fltarr(maxpos-minpos+ydimsize,total_lines)
      stline = 0
      endline = xdimsize - 1

      t[0,0] = fdata[*,*]

      for j=0,num_blocks -1 do begin
        if (offarray[j] ne 0) then begin
          t[0,stline] =  shift(t[*,stline:endline],offarray[j],0)
        endif
        stline = endline + 1
        endline = endline + xdimsize
      endfor

      agpdata[0,0,f] = t

    endfor  ; loop through the AGP fields

  status = EOS_GD_DETACH(agpgid)
  status = EOS_GD_CLOSE(agpfid)

  ;  set up map info for this subregion with ENVI_MAP_INFO_CREATE

  projname = "MISR AGP PATH "+ strtrim(string(path),2) + " BLOCK " + $
               strtrim(string(start_block),2)
  descrip = FILE_BASENAME(agpfile,'.hdf')+"_B"+strtrim(string(start_block),2)+"-"+$
              strtrim(string(end_block),2)

  som_block_corner_meters,agpfile,start_block,somupleft,somlowright

  ;  projection parameters needed for the generic SOM A projection:
  ;  a, b, IncAng, AscLong, x0, y0, PSRev, LRat

  somprojparm = [6378137.0, 6356752.3, 98.303819, dms_to_ddeg(projparam[4]), $
                 0.0, 0.0, 98.880000, 0.000000]

  map_info = ENVI_MAP_INFO_CREATE(DATUM="WGS-84",PS=[pixsize,pixsize],$
             ROTATION=90.0,NAME=projname,$
             MC=[offarray[0],0,somupleft[0],somupleft[1]],PARAMS=somprojparm,TYPE=37)

  ;  use ENVI_ENTER_DATA to add data to the Available Bands List

  ENVI_ENTER_DATA,agpdata,BNAMES=agpfields,DEF_STRETCH=ENVI_DEFAULT_STRETCH_CREATE(/NONE),$
    MAP_INFO=map_info,DESCRIP=descrip

endif ; AGP info available

end