PRO MODIS_EXTRACT, GEOFILE, RADFILE, CLAT, CLON, $
  REFLECTANCE=REFLECTANCE, TEMPERATURE=TEMPERATURE, WIDTH=WIDTH, $
  OUTFILE=OUTFILE, UTC=UTC, SUCCESS=SUCCESS

;+
; NAME:
;    MODIS_EXTRACT
;
; PURPOSE:
;    Given a MODIS geolocation file (MOD03) and a MODIS 1000 meter resolution
;    radiance data file (MOD021KM), extract a box of pixels from all MODIS
;    bands centered on a given latitude and longitude.
;
;    If the requested location is not found, no data are returned.
;
;    Both GSFC DAAC operational Level-1B and IMAPP direct broadcast Level-1B
;    data formats are recognized.
;
; CATEGORY:
;    MODIS utilities.
;
; CALLING SEQUENCE:
;    MODIS_EXTRACT, GEOFILE, RADFILE, CLAT, CLON
;
; INPUTS:
;    GEOFILE        Name of MODIS geolocation file
;    RADFILE        Name of MODIS 1000 meter resolution radiance file
;    CLAT           Latitude at center of region to extract (degrees N)
;    CLON           Longitude at center of region to extract (degrees E)
;
; OPTIONAL INPUTS:
;    None.
;
; KEYWORD PARAMETERS:
;    REFLECTANCE    If set, data for bands 1-19 and 26 are returned as
;                   TOA reflectances *with* solar zenith angle correction
;                   (default is to return radiances in units of
;                   Watts/meter^2/steradian/micron *without* solar zenith
;                   angle correction).
;    TEMPERATURE    If set, data for bands 20-25 and 27-36 are returned as
;                   brightness temperatures in units of Kelvin
;                   (default is to return radiances in units of
;                   Watts/meter^2/steradian/micron).
;    WIDTH          Width of box of pixels
;                   (default 3, must be an odd positive integer)
;    OUTFILE        Name of output file
;                   (default is to send output to screen)
;
;    
; OUTPUTS:
;    A text listing of the data extracted for the requested location.
;
;    The output fields are self-explanatory, with the exception of 'Land/Sea Mask',
;    which is explained below:
;    0: Shallow Ocean (Ocean <5k from coast OR <50m deep).
;    1: Land (not anything else).
;    2: Ocean Coastlines and Lake Shorelines.
;    3: Shallow Inland Water (Inland Water < 5km from shore OR < 50m deep).
;    4: Ephemeral (intermittent) Water.
;    5: Deep Inland Water (Inland water > 5km from shoreline AND > 50m deep).
;    6: Moderate or Continental Ocean (Ocean > 5km from coast AND > 50m deep AND < 500m deep).
;    7: Deep Ocean (Ocean > 500m deep).
;
; OPTIONAL OUTPUTS:
;    None.
;
; COMMON BLOCKS:
;    None
;
; SIDE EFFECTS:
;    None.
;
; RESTRICTIONS:
;    Requires IDL 5.0 or higher (square bracket array syntax).
;
;    Requires the following MODIS procedures by Liam.Gumley@ssec.wisc.edu:
;    MODIS_FILE_INFO     Get information about a MODIS HDF file
;    MODIS_LEVEL1B_READ  Read a single band from a MODIS HDF radiance image file
;
;    Requires the following HDF procedures by Liam.Gumley@ssec.wisc.edu:
;    HDF_SD_ATTINFO      Get information about an attribute
;    HDF_SD_ATTLIST      Get list of attribute names
;    HDF_SD_VARINFO      Get information about an SDS variable
;    HDF_SD_VARLIST      Get a list of SDS variable names
;    HDF_SD_VARREAD      Read an SDS variable
;    HDF_VD_VDATAINFO    Get information about a Vdata
;    HDF_VD_VDATALIST    Get list of Vdata names
;    HDF_VD_VDATAREAD    Read a Vdata field
;
;    Requires the following Planck procedures by Liam.Gumley@ssec.wisc.edu:
;    MODIS_BRIGHT        Compute MODIS brightness temperature
;    MODIS_PLANCK        Compute MODIS Planck radiance
;    BRIGHT_M            Compute brightness temperature (EOS radiance units)
;    BRITE_M             Compute brightness temperature (UW-SSEC radiance units)
;    PLANCK_M            Compute monochromatic Planck radiance (EOS units)
;    PLANC_M             Compute monochromatic Planck radiance (UW-SSEC units)
;
;    Requires the Solarsoft time library routines:
;    http://sohowww.nascom.nasa.gov/solarsoft/gen/idl/time
;
; EXAMPLES:
;
;; Extract a 3x3 box of pixels in Charlotte Harbor in SW Florida on 2000/12/14
;geofile = 'MOD03.A2000349.1610.002.2001003054620.hdf'
;radfile = 'MOD021KM.A2000349.1610.002.2001003080203.hdf'
;clat = 26.75
;clon = -82.11
;modis_extract, geofile, radfile, 26.75, -82.11, outfile='test.out'
;modis_extract, geofile, radfile, 26.75, -82.11, /temperature, outfile='test.out'
;
; MODIFICATION HISTORY:
; Liam.Gumley@ssec.wisc.edu
; http://cimss.ssec.wisc.edu/~gumley
; $Id: modis_extract.pro,v 1.1 2003/06/30 20:27:21 gumley Exp $
;
; Copyright (C) 1999, 2000 Liam E. Gumley
;
; This program is free software; you can redistribute it and/or
; modify it under the terms of the GNU General Public License
; as published by the Free Software Foundation; either version 2
; of the License, or (at your option) any later version.
;
; This program is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
; GNU General Public License for more details.
;
; You should have received a copy of the GNU General Public License
; along with this program; if not, write to the Free Software
; Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
;-

rcs_id = '$Id: modis_extract.pro,v 1.1 2003/06/30 20:27:21 gumley Exp $'

;- Check arguments
if (n_elements(geofile) eq 0) then message, 'Argument GEOFILE is undefined'
if (n_elements(radfile) eq 0) then message, 'Argument RADFILE is undefined'
if (n_elements(clat) eq 0)    then message, 'Argument CLAT is undefined'
if (n_elements(clon) eq 0)    then message, 'Argument CLON is undefined'
if (n_elements(width) eq 0)   then width = 3
if ((width mod 2) ne 1) then begin
  message, 'WIDTH must be an odd positive integer', /continue
  return
endif

;- Save start time
start_time = systime()

;- Print header on-screen
print, ''
print, rcs_id
print, 'Started at ', start_time
print, ''

;- Set success flag
success = 0

;-------------------------------------------------------------------------------
;- CHECK INPUT FILES
;-------------------------------------------------------------------------------

;- Check validity of geolocation file
geoinfo = modis_file_info(geofile, /quiet)
if (geoinfo.read eq 0) then begin
  message, 'Geolocation file cannot not be opened: ' + geofile, /continue
  return
endif
if (geoinfo.nvars eq 0) then begin
  message, 'Geolocation file contains no variables: ' + geofile, /continue
  return
endif
index = where(geoinfo.var_names eq 'Land/SeaMask', count)
if (count eq 0) then begin
  message, 'Geolocation file format is incorrect: ' + geofile, /continue
  return
endif

;- Check validity of radiance file
radinfo = modis_file_info(radfile, /quiet)
if (radinfo.read eq 0) then begin
  message, 'Radiance file cannot be opened: ' + radfile, /continue
  return
endif
index = where(radinfo.var_names eq 'EV_1KM_Emissive', count)
if (count eq 0) then begin
  message, 'Radiance file format is incorrect: ' + geofile, /continue
  return
endif

;- Check for matching start times in geolocation and radiance files
hdfid = hdf_sd_start(geofile)
hdf_sd_varread, hdfid, 'EV start time', geotime
hdf_sd_end, hdfid
hdfid = hdf_open(radfile)
hdf_vd_vdataread, hdfid, $
  'Level 1B Swath Metadata', 'EV Sector Start Time', radtime
hdf_close, hdfid
if (abs(geotime[0] - radtime[0]) gt 1.0e-6) then begin
  message, 'Start times of geolocation and radiance files do not match', /continue
  return
endif

;-------------------------------------------------------------------------------
;- READ GEOLOCATION DATA AND FIND REQUESTED LAT/LON
;-------------------------------------------------------------------------------

;- Get latitude, longitude, and land/sea mask arrays
hdfid = hdf_sd_start(geofile)
hdf_sd_varread, hdfid, 'Latitude', lat
hdf_sd_varread, hdfid, 'Longitude', lon
hdf_sd_varread, hdfid, 'Land/SeaMask', msk
hdf_sd_end, hdfid

;- Find the requested lat/lon
modis_locate, clat, clon, lat, lon, x, y, minrange=minrange
if ((x lt 0) or (y lt 0) or (minrange gt 5.0)) then begin
  message, 'Requested lat/lon was not found', /continue
  return
endif

;- Open output file if required
if (n_elements(outfile) eq 1) then begin
  openw, outlun, outfile, /get_lun, width=350
endif else begin
  outlun = -1
endelse

;- Print header to file
if (outlun gt 0) then begin
  printf, outlun, ''
  printf, outlun, rcs_id
  printf, outlun, 'Started at ', start_time
  printf, outlun, ''
endif

;- Print input parameters
printf, outlun, 'Geolocation file: ', geofile
printf, outlun, 'Radiance file   : ', radfile
printf, outlun, clat, format='("Requested latitude  (deg N): ", f10.4)'
printf, outlun, clon, format='("Requested longitude (deg E): ", f10.4)'
printf, outlun, ''
  
;- Print location of actual pixel
printf, outlun, lat[x, y], format='("Actual pixel latitude    (deg N): ", f10.4)'
printf, outlun, lon[x, y], format='("Actual pixel longitude   (deg E): ", f10.4)'
printf, outlun, minrange,  format='("Range to requested lat/lon  (km): ", f6.2)'
printf, outlun, x,         format='("Actual pixel number (zero-based): ", i5)'
printf, outlun, y,         format='("Actual line number  (zero-based): ", i5)'
printf, outlun, ''

;-------------------------------------------------------------------------------
;- GET ANCILLARY DATA AT REQUESTED LAT/LON
;-------------------------------------------------------------------------------

;- Get time code
hdfid = hdf_sd_start(geofile)
hdf_sd_varread, hdfid, 'EV start time', time, start=[y / 10L], count=[1]
hdf_sd_end, hdfid

;- Convert ECS toolkit TAI time to UTC time
ref = '1993-01-01T00:00:00.000Z'
utc = tai2utc(time + utc2tai(ref), /ecs)
printf, outlun, 'Start time of scan (UTC): ', utc
printf, outlun, ''

;- Get zenith and azimuth angles, surface elevation, and surface type
hdfid = hdf_sd_start(geofile)
hdf_sd_varread, hdfid, 'SolarZenith',   solzen, start=[x, y], count=[1, 1]
hdf_sd_varread, hdfid, 'SolarAzimuth',  solazm, start=[x, y], count=[1, 1]
hdf_sd_varread, hdfid, 'SensorZenith',  senzen, start=[x, y], count=[1, 1]
hdf_sd_varread, hdfid, 'SensorAzimuth', senazm, start=[x, y], count=[1, 1]
hdf_sd_varread, hdfid, 'Height',        height, start=[x, y], count=[1, 1]
hdf_sd_varread, hdfid, 'Land/SeaMask',  lansea, start=[x, y], count=[1, 1]
hdf_sd_end, hdfid
solzen = 0.01 * solzen
solazm = 0.01 * solazm
senzen = 0.01 * senzen
senazm = 0.01 * senazm
printf, outlun, solzen, format='("Solar zenith   (deg): ", f7.2)'
printf, outlun, solazm, format='("Solar azimuth  (deg): ", f7.2)'
printf, outlun, senzen, format='("Sensor zenith  (deg): ", f7.2)'
printf, outlun, senazm, format='("Sensor azimuth (deg): ", f7.2)'
printf, outlun, height, format='("Terrain height   (m): ", i7)'
printf, outlun, lansea, format='("Land/SeaMask   (0-7): ", i1)'
printf, outlun, ''

;- Print actual pixel and line numbers extracted
modis_level1b_read, radfile, 31, data, $
  area=[x - width / 2, y - width / 2, width, width], $
  start=start, count=count
printf, outlun, start[0], format='("Start pixel (zero-based): ", i6)'
printf, outlun, start[1], format='("Start line  (zero-based): ", i6)'
printf, outlun, count[0], format='("Number of pixels        : ", i6)'
printf, outlun, count[1], format='("Number of lines         : ", i6)'
printf, outlun, ''

;- Print lat/lons and land/sea mask for pixels extracted
fmt = strcompress('(' + string(count[0]) + 'f10.4)', /remove_all)
printf, outlun, 'Latitude (deg N)'
printf, outlun, lat[start[0] : start[0] + count[0] - 1, start[1] : start[1] + count[1] - 1], $
  format=fmt
printf, outlun, ''  
printf, outlun, 'Longitude (deg E)'
printf, outlun, lon[start[0] : start[0] + count[0] - 1, start[1] : start[1] + count[1] - 1], $
  format=fmt
printf, outlun, ''  
fmt = strcompress('(' + string(count[0]) + 'i3)', /remove_all)
printf, outlun, 'Land/Sea Mask'
printf, outlun, msk[start[0] : start[0] + count[0] - 1, start[1] : start[1] + count[1] - 1], $
  format='
printf, outlun, ''

;- Get solar zenith angles for extracted pixels
hdfid = hdf_sd_start(geofile)
hdf_sd_varread, hdfid, 'SolarZenith', zen, start=start[0:1], count=count[0:1]
hdf_sd_end, hdfid
zen = zen * 0.01

;- Get mirror side information for extracted pixels
hdfid = hdf_sd_start(geofile)
hdf_sd_varread, hdfid, 'Mirror side', mirror_side
hdf_sd_end, hdfid
nscans = n_elements(mirror_side)
mirror_side = rebin(mirror_side, 10L * nscans, /sample)
printf, outlun, 'Scan Mirror Side'
printf, outlun, mirror_side[start[1] : start[1] + count[1] - 1], format='(i3)'
printf, outlun, ''

;- Print detector numbers for extracted pixels
printf, outlun, 'Detector Number'
for line = start[1], start[1] + count[1] - 1 do begin
  printf, outlun, line mod 10L, format='(i3)'
endfor
printf, outlun, ''

;-------------------------------------------------------------------------------
;- GET FORMAT AND VERSION INFORMATION
;-------------------------------------------------------------------------------

;- Get format information from radiance file
hdfid = hdf_sd_start(radfile)
attlist = hdf_sd_attlist(hdfid, '', /global)
index = where(attlist.attnames eq 'CoreMetadata.0', count)
if (count eq 1) then begin
  fmt = 'DAAC operational'
  attinfo = hdf_sd_attinfo(hdfid, '', 'CoreMetadata.0', /global)
  pge_ver = get_metadata(attinfo.data, 'PGEVERSION')
endif else begin
  fmt = 'IMAPP direct broadcast'
  pge_ver = 'unknown'
endelse

;- Get calibration lookup table version information
attname = 'Reflective LUT Serial Number and Date of Last Change'
attinfo = hdf_sd_attinfo(hdfid, '', attname, /global)
srb_lut_ver = attinfo.data
attname = 'Emissive LUT Serial Number and Date of Last Change'
attinfo = hdf_sd_attinfo(hdfid, '', attname, /global)
teb_lut_ver = attinfo.data
hdf_sd_end, hdfid

;- Print format and version information
printf, outlun, 'Radiance file format  : ', fmt
printf, outlun, 'PGE version           : ', pge_ver
printf, outlun, 'Reflective LUT version: ', srb_lut_ver
printf, outlun, 'Emissive LUT version  : ', teb_lut_ver
printf, outlun, ''
    
;-------------------------------------------------------------------------------
;- EXTRACT AND PRINT RADIANCE DATA
;-------------------------------------------------------------------------------

;- Print output units
if keyword_set(reflectance) then begin
  printf, outlun, 'Units for bands 1-19 and 26   : Reflectance (solar zenith corrected)'
endif else begin
  printf, outlun, 'Units for bands 1-19 and 26   : Radiance (Watts/meter^2/steradian/micron)'
endelse
if keyword_set(temperature) then begin
  printf, outlun, 'Units for bands 20-36 (not 26): Brightness temperature (K)'
endif else begin
  printf, outlun, 'Units for bands 20-36 (not 26): Radiance (Watts/meter^2/steradian/micron)'
endelse
printf, outlun, ''
  
;- Extract radiances for all bands
for band = 1, 36 do begin

  ;- Get radiance data
  if ((band le 19) or (band eq 26)) then begin
  
    ;- Get a solar reflection band
    modis_level1b_read, radfile, band, data, $
      area=[x - width / 2, y - width / 2, width, width], $
      reflectance=keyword_set(reflectance), $
      valid_index=valid_index, valid_count=valid_count

    ;- Apply solar zenith angle correction if necessary
    if (keyword_set(reflectance) and (valid_count gt 0)) then $
      data[valid_index] = data[valid_index] / cos(!dtor * zen[valid_index])
      
  endif else begin

    ;- Get a thermal emission band
    modis_level1b_read, radfile, band, data, $
      area=[x - width / 2, y - width / 2, width, width], $
      temperature=keyword_set(temperature), $
      valid_index=valid_index, valid_count=valid_count

  endelse

  ;- Compute statistics
  if (valid_count ge 2) then begin
    result = moment(data[valid_index], sdev=sdev)
    mean = result[0]
  endif else begin
    mean = -1.0
    sdev = -1.0
  endelse

  ;- Print data for this band
  printf, outlun, band, $
    format='("Band: ", i2.2)'
  dims = size(data, /dimensions)
  fmt = strcompress('(' + string(dims[0]) + 'f8.3)', /remove_all)  
  printf, outlun, data, format=fmt
  printf, outlun, valid_count, mean, sdev, $
    format='("Nvalid: ", i4, ", Mean: ", f10.5, ", Std Dev: ", f10.5)'
  printf, outlun, ''

endfor

;-------------------------------------------------------------------------------
;- CLEAN UP
;-------------------------------------------------------------------------------

;- Close the output file if needed and set success flag
if (outlun gt 0) then free_lun, outlun
success = 1

END