FUNCTION GET_SDS_NAME, BAND

;- Get the MODIS 1 km SDS name for the requested band
case 1 of
  (band ge  1 and band le  2) : sds_name = 'EV_250_Aggr1km_RefSB'
  (band ge  3 and band le  7) : sds_name = 'EV_500_Aggr1km_RefSB'
  (band ge  8 and band le 19) : sds_name = 'EV_1KM_RefSB'
  (band eq 26)                : sds_name = 'EV_1KM_RefSB'
  (band ge 20 and band le 36) : sds_name = 'EV_1KM_Emissive'
  else : message, 'Invalid band number'
endcase
return, sds_name

END
;-------------------------------------------------------------------------------
FUNCTION GET_BAND_INDEX, HDFID, BAND

;- Get the index of the requested band in the corresponding MODIS 1 km SDS

;- Get the SDS name
sds_name = get_sds_name(band)

;- Get the list of band numbers for this SDS array
att_info = hdf_sd_attinfo(hdfid, sds_name, 'band_names')
if (att_info.name eq '') then message, 'Attribute not found: band_names'
band_list = str_sep(att_info.data, ',')

;- Convert the requested band number to a string
band_name = strcompress(long(band), /remove_all)
if (band eq 13) then band_name = '13lo'
if (band eq 14) then band_name = '14lo'
if (abs(band - 13.5) lt 0.01) then band_name = '13hi'
if (abs(band - 14.5) lt 0.01) then band_name = '14hi'

;- Get the index of the band
band_index = (where(band_name eq band_list))[0]
if (band_index lt 0) then message, 'Requested band number was not found'
return, band_index

END
;-------------------------------------------------------------------------------
PRO READ_MODIS_RADIANCE, HDFID, BAND, NCOL, NROW, DATA, RAW=RAW

;- Read a single band of MODIS 1 km radiance data

;- Select the SDS
sds_name = get_sds_name(band)
varid = hdf_sd_select(hdfid, hdf_sd_nametoindex(hdfid, sds_name))

;- Get the band index
band_index = get_band_index(hdfid, band)

;- Get the start and count vectors
start = [0L, 0L, band_index]
count = [ncol, nrow, 1L]

;- Read the scaled integers
hdf_sd_getdata, varid, data, start=start, count=count

;- Get the valid range, scales, and offsets
hdf_sd_attrinfo, varid, hdf_sd_attrfind(varid, 'valid_range'), data=range
hdf_sd_attrinfo, varid, hdf_sd_attrfind(varid, 'radiance_scales'), data=scale
hdf_sd_attrinfo, varid, hdf_sd_attrfind(varid, 'radiance_offsets'), data=offset

;- Deselect the SDS
hdf_sd_endaccess, varid

;- Convert scaled integers to radiance
if (keyword_set(raw) eq 0) then $
  data = scale[band_index] * (data - offset[band_index])

END
;-------------------------------------------------------------------------------
PRO WRITE_MODIS_RADIANCE, HDFID, BAND, NCOL, NROW, DATA, RAW=RAW

;- Write a single band of MODIS 1 km radiance data

;- Select the SDS
sds_name = get_sds_name(band)
varid = hdf_sd_select(hdfid, hdf_sd_nametoindex(hdfid, sds_name))

;- Get the band index
band_index = get_band_index(hdfid, band)

;- Get the start and count vectors
start = [0L, 0L, band_index]
count = [ncol, nrow, 1L]

;- Get the scales and offsets
hdf_sd_attrinfo, varid, hdf_sd_attrfind(varid, 'radiance_scales'), data=scale
hdf_sd_attrinfo, varid, hdf_sd_attrfind(varid, 'radiance_offsets'), data=offset

;- Convert radiances to scaled integers
if (keyword_set(raw) eq 0) then begin
  scaled = uint((data / scale[band_index]) + offset[band_index] + 0.5)
endif else begin
  scaled = data
endelse

;- Write the scaled integers
hdf_sd_adddata, varid, scaled, start=start, count=count

;- Deselect the SDS
hdf_sd_endaccess, varid

END
;-------------------------------------------------------------------------------
PRO HDF_DESTRIPE_NEW, RADFILE, SATELLITE=SATELLITE, QUIET=QUIET

;- Set version string
rcs_id = '$Id: hdf_destripe_new.pro,v 1.10 2004/05/26 21:45:37 gumley Exp $'
if (keyword_set(quiet) eq 0) then print, rcs_id

;- Check arguments
if (n_elements(radfile) eq 0) then message, 'Argument RADFILE is undefined'
if (n_elements(satellite) eq 0) then begin
  message, 'Satellite name not set: assuming TERRA', /continue
  satellite = 'TERRA'
endif

;-------------------------------------------------------------------------------
;- SET DESTRIPING CONFIGURATION
;-------------------------------------------------------------------------------

;- Terra MODIS destriping configuration
terra_cfg = [$
[20, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[22, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[23, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[24, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[25, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[26, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[27, 4, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0], $
[28, 4, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1], $
[29, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[30, 5, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0], $
[31, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[32, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[33, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[34, 3, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0], $
[35, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[36, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]

;- Aqua MODIS destriping configuration
aqua_cfg = [$
[20, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[22, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[23, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[24, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[25, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[26, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[27, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[28, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[29, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[30, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[31, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[32, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[33, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[34, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[35, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], $
[36, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]

;- Select satellite
case strupcase(satellite) of
  'TERRA' : begin
    band_list        = reform(terra_cfg[0, *])
    ref_list         = reform(terra_cfg[1, *])
    detector_replace = terra_cfg[2:*, *]
  end
  'AQUA' : begin
    band_list        = reform(aqua_cfg[0, *])
    ref_list         = reform(aqua_cfg[1, *])
    detector_replace = aqua_cfg[2:*, *]
  end
  else : message, 'Satellite name not recognized => ' + strcompress(satellite)
endcase

;-------------------------------------------------------------------------------
;- GET ANCILLARY INFORMATION
;-------------------------------------------------------------------------------

;- Get the mirror side data (in Vdata mode)
hdfid = hdf_open(radfile)
index = hdf_vd_find(hdfid, 'Level 1B Swath Metadata')
vdataid = hdf_vd_attach(hdfid, index)
nread = hdf_vd_read(vdataid, mirror_side, fields='Mirror Side')
hdf_vd_detach, vdataid
hdf_close, hdfid

;- Open the file (in SDS mode)
hdfid = hdf_sd_start(radfile, /rdwr)

;- Get the number of rows and columns
sds_name = 'EV_1KM_Emissive'
varinfo = hdf_sd_varinfo(hdfid, sds_name)
ncol = varinfo.dims[0]
nrow = varinfo.dims[1]
nscan = nrow / 10L

;- Rebin the mirror side data (1 value per scan) to the number of rows
mirror_side = rebin(mirror_side, nrow, /sample)

;-------------------------------------------------------------------------------
;- CHECK IF FILE IS ALREADY DESTRIPED
;-------------------------------------------------------------------------------

;- Get all global attributes
hdf_sd_fileinfo, hdfid, nvars, natts
att_names = strarr(natts)
for index = 0L, natts - 1L do begin
  hdf_sd_attrinfo, hdfid, index, name=name
  att_names[index] = name
endfor

;- Check for specific global attribute
loc = where(att_names eq 'UW_DESTRIPE_LWIR', count)
if (count gt 0) then begin
  hdf_sd_end, hdfid
  message, 'This file is already destriped => ' + radfile, /continue
  return
endif

;-------------------------------------------------------------------------------
;- DESTRIPE INFRARED BANDS
;-------------------------------------------------------------------------------

;- Set print flag for destripe procedure
quiet = 0

;- Start loop over bands
for band_index = 0L, n_elements(band_list) - 1L do begin

  ;- Get the band number
  band = band_list[band_index]
  
  if (keyword_set(quiet) eq 0) then $
    print, 'Processing band ', strcompress(band)
  
  ;- Get the reference detector number
  ref = ref_list[band_index]
   
  ;- Read the radiance data
  read_modis_radiance, hdfid, band, ncol, nrow, data, /raw

  ;- Destripe this band if at least one scan is good
  loc = where(data le 32767, count)
  if (count gt 13540) then begin
  
    ;- Destripe the data
    data = modis_edf_destripe(data, ref, mirror_side, version=version)

    ;- Replace bad detectors with nearest good neighbor
    for detector_index = 0L, 9L do begin

      ;- Check if this detector should be replaced
      if (detector_replace[detector_index, band_index]) then begin

        ;- Get the index of the replacement
        loc = where(detector_replace[*, band_index] eq 0)
        val = min(abs(detector_index - loc), min_index)
        replace_index = loc[min_index]
        if (keyword_set(quiet) eq 0) then $
          print, detector_index, replace_index, $
            format='(2x, "Replacing detector ", i2, " with detector ", i2)'

        ;- Replace the data
        data[*, (lindgen(nscan) * 10L) + detector_index] = $
          data[*, (lindgen(nscan) * 10L) + replace_index]

      endif

    endfor

    ;- Write the destriped radiances
    write_modis_radiance, hdfid, band, ncol, nrow, data, /raw

  endif else begin
  
    ;- Print diagnostic message
    print, '  Skipping this band because of insufficient valid data'

  endelse
  
endfor

;-------------------------------------------------------------------------------
;- CLOSE THE FILE
;-------------------------------------------------------------------------------

;- Write a new global attribute to show this granule has been destriped
hdf_sd_attrset, hdfid, 'UW_DESTRIPE_LWIR', rcs_id + ';' + version

;- Close the HDF file
hdf_sd_end, hdfid

if (keyword_set(quiet) eq 0) then $
  print, 'Successfully destriped ' + radfile

END