PRO MODIS_LEVEL1B_READ, FILENAME, BAND, IMAGE, $ RAW=RAW, REFLECTANCE=REFLECTANCE, TEMPERATURE=TEMPERATURE, $ CORRECTED=CORRECTED, $ AREA=AREA, UNITS=UNITS, PARAMETER=PARAMETER, $ SCANTIME=SCANTIME, MIRROR=MIRROR, LATITUDE=LATITUDE, LONGITUDE=LONGITUDE, $ VALID_INDEX=VALID_INDEX, VALID_COUNT=VALID_COUNT, VALID_RANGE=VALID_RANGE, $ INVALID_INDEX=INVALID_INDEX, INVALID_COUNT=INVALID_COUNT, RANGE=RANGE, $ MISSING=MISSING, START=START, COUNT=COUNT, SCALE_FACTOR=SCALE_FACTOR ;+ ; NAME: ; MODIS_LEVEL1B_READ ; ; PURPOSE: ; Read a single band from a MODIS Level 1B HDF product file at ; 1000, 500, or 250 m resolution. ; ; The output image is available in the following units (Radiance is default): ; RAW DATA: Raw data values as they appear in the HDF file ; RADIANCE: Radiance (Watts per square meter per steradian per micron) ; REFLECTANCE: Reflectance (Bands 1-19 and 26; *without* solar zenith correction) ; TEMPERATURE: Brightness Temperature (Bands 20-25 and 27-36, Kelvin) ; ; This procedure uses only HDF calls (it does *not* use HDF-EOS), ; and only reads from SDS and Vdata arrays (it does *not* read ECS metadata). ; ; CATEGORY: ; MODIS utilities. ; ; CALLING SEQUENCE: ; MODIS_LEVEL1B_READ, FILENAME, BAND, IMAGE ; ; INPUTS: ; FILENAME Name of MODIS Level 1B HDF file ; BAND Band number to be read ; (1-36 for 1000 m, 1-7 for 500 m, 1-2 for 250m) ; (Use 13, 13.5, 14, 14.5 for 13lo, 13hi, 14lo, 14hi) ; ; OPTIONAL INPUTS: ; None. ; ; KEYWORD PARAMETERS: ; RAW If set, image data are returned as raw HDF values ; (default is to return image data as radiance). ; REFLECTANCE If set, image data for bands 1-19 and 26 only are ; returned as reflectance *without* solar zenith angle correction ; (default is to return image data as radiance). ; TEMPERATURE If set, image data for bands 20-25 and 27-36 only are ; returned as brightness temperature ; (default is to return image data as radiance). ; AREA A four element vector specifying the area to be read, ; in the format [X0,Y0,NX,NY] ; (default is to read the entire image). ; UNITS On return, a string describing the image units. ; PARAMETER On return, a string describing the image (e.g. 'Radiance'). ; SCANTIME On return, a vector containing the start time of each scanline ; (SDPTK TAI seconds). ; MIRROR On return, a vector containing the mirror side for each scanline ; LATITUDE On return, an array containing the reduced resolution latitude ; data for the entire granule (degrees, every 5th line and pixel). ; LONGITUDE On return, an array containing the reduced resolution longitude ; data for the entire granule (degrees, every 5th line and pixel). ; VALID_INDEX On return, a vector containing the 1D indexes of pixels which ; are within the 'valid_range' attribute values. ; VALID_COUNT On return, the number of pixels which ; are within the 'valid_range' attribute values. ; INVALID_INDEX On return, a vector containing the 1D indexes of pixels which ; are not within the 'valid_range' attribute values. ; INVALID_COUNT On return, the number of pixels which ; are not within the 'valid_range' attribute values. ; RANGE On return, a 2-element vector containing the minimum and ; maximum data values within the 'valid range'. ; MISSING If set, defines the value to be inserted for missing data when returning ; the result in radiance, reflectance, or brightness temperature units ; (default is -1.0). ; ; OUTPUTS: ; IMAGE A two dimensional array of image data in the requested units. ; ; OPTIONAL OUTPUTS: ; None. ; ; COMMON BLOCKS: ; None ; ; SIDE EFFECTS: ; None. ; ; RESTRICTIONS: ; Requires IDL 5.0 or higher (square bracket array syntax). ; ; 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 ; BRIGHT_M Compute brightness temperature (EOS radiance units) ; ; EXAMPLES: ; ;; These examples require the IMDISP image display procedure, available from ;; http://cimss.ssec.wisc.edu/~gumley/imdisp.html ; ;; Read band 1 in radiance units from a 1000 m resolution file ;file = 'MOD021KM.A2000062.1020.002.2000066023928.hdf' ;modis_level1b_read, file, 1, band01 ;imdisp, band01, margin=0, order=1 ; ;; Read band 31 in temperature units from a 1000 m resolution file ;file = 'MOD021KM.A2000062.1020.002.2000066023928.hdf' ;modis_level1b_read, file, 31, band31, /temperature ;imdisp, band31, margin=0, order=1, range=[285, 320] ; ;; Read a band 1 subset in reflectance units from a 500 m resolution file ;file = 'MOD02HKM.A2000062.1020.002.2000066023928.hdf' ;modis_level1b_read, file, 1, band01, /reflectance, area=[1000, 1000, 700, 700] ;imdisp, band01, margin=0, order=1 ; ;; Read band 6 in reflectance units from a 1000 m resolution file, ;; and screen out invalid data when scaling ;file = 'MOD021KM.A2000062.1020.002.2000066023928.hdf' ;modis_level1b_read, file, 6, band06, /reflectance, valid_index=valid_index ;range = [min(band06[valid_index]), max(band06[valid_index])] ;imdisp, band06, margin=0, order=1, range=range ; ; MODIFICATION HISTORY: ; Liam.Gumley@ssec.wisc.edu ; http://cimss.ssec.wisc.edu/~gumley ; $Id: modis_level1b_read.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_level1b_read.pro,v 1.1 2003/06/30 20:27:21 gumley Exp $' ;------------------------------------------------------------------------------- ;- CHECK INPUT ;------------------------------------------------------------------------------- ;- Check arguments if (n_params() ne 3) then $ message, 'Usage: MODIS_LEVEL1B_READ, FILENAME, BAND, IMAGE' if (n_elements(filename) eq 0) then $ message, 'Argument FILENAME is undefined' if (n_elements(band) eq 0) then $ message, 'Argument BAND is undefined' if (arg_present(image) ne 1) then $ message, 'Argument IMAGE cannot be modified' if (n_elements(area) gt 0) then begin if (n_elements(area) ne 4) then $ message, 'AREA must be a 4-element vector of the form [X0,Y0,NX,NY]' endif if (n_elements(missing) eq 0) then missing = -1.0 ;------------------------------------------------------------------------------- ;- CHECK FOR VALID MODIS L1B HDF FILE, AND GET FILE TYPE (1km, 500 m, or 250 m) ;------------------------------------------------------------------------------- ;- Check that file exists if ((findfile(filename))[0] eq '') then $ message, 'FILENAME was not found => ' + filename ;- Get expanded filename openr, lun, filename, /get_lun fileinfo = fstat(lun) free_lun, lun ;- Check that file is HDF if (hdf_ishdf(fileinfo.name) ne 1) then $ message, 'FILENAME is not HDF => ' + fileinfo.name ;- Get list of SDS arrays sd_id = hdf_sd_start(fileinfo.name) varlist = hdf_sd_varlist(sd_id) hdf_sd_end, sd_id ;- Locate image arrays index = where(varlist.varnames eq 'EV_1KM_Emissive', count_1km) index = where(varlist.varnames eq 'EV_500_RefSB', count_500) index = where(varlist.varnames eq 'EV_250_RefSB', count_250) case 1 of (count_1km eq 1) : filetype = 'MOD021KM' (count_500 eq 1) : filetype = 'MOD02HKM' (count_250 eq 1) : filetype = 'MOD02QKM' else : message, 'FILENAME is not MODIS Level 1B HDF => ' + fileinfo.name endcase ;------------------------------------------------------------------------------- ;- CHECK BAND NUMBER, AND KEYWORDS WHICH DEPEND ON BAND NUMBER ;------------------------------------------------------------------------------- ;- Check band number case filetype of 'MOD021KM' : if (band lt 1) or (band gt 36) then $ message, 'BAND range is 1-36 for this MODIS type => ' + filetype 'MOD02HKM' : if (band lt 1) or (band gt 7) then $ message, 'BAND range is 1-7 for this MODIS type => ' + filetype 'MOD02QKM' : if (band lt 1) or (band gt 2) then $ message, 'BAND range is 1-2 for this MODIS type => ' + filetype endcase ;- Check for invalid request for reflectance units if ((band ge 20) and (band ne 26)) and keyword_set(reflectance) then $ message, 'REFLECTANCE units valid for bands 1-19, 26 only' ;- Check for invalid request for temperature units if ((band le 19) or (band eq 26)) and keyword_set(temperature) then $ message, 'TEMPERATURE units valid for bands 20-25, 27-36 only' ;------------------------------------------------------------------------------- ;- SET VARIABLE NAME FOR IMAGE DATA ;------------------------------------------------------------------------------- case filetype of ;- 1 km resolution data 'MOD021KM' : begin 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_Band26' (band ge 20 and band le 36) : sds_name = 'EV_1KM_Emissive' endcase end ;- 500 m resolution data 'MOD02HKM' : begin case 1 of (band ge 1 and band le 2) : sds_name = 'EV_250_Aggr500_RefSB' (band ge 3 and band le 7) : sds_name = 'EV_500_RefSB' endcase end ;- 250 m resolution data 'MOD02QKM' : sds_name = 'EV_250_RefSB' endcase ;------------------------------------------------------------------------------- ;- SET ATTRIBUTE NAMES FOR IMAGE DATA ;------------------------------------------------------------------------------- ;- Set names of scale, offset, and units attributes if keyword_set(reflectance) then begin scale_name = 'reflectance_scales' offset_name = 'reflectance_offsets' units_name = 'reflectance_units' parameter = 'Reflectance' endif else begin scale_name = 'radiance_scales' offset_name = 'radiance_offsets' units_name = 'radiance_units' parameter = 'Radiance' endelse ;------------------------------------------------------------------------------- ;- OPEN THE FILE IN SDS MODE ;------------------------------------------------------------------------------- sd_id = hdf_sd_start(fileinfo.name) ;------------------------------------------------------------------------------- ;- CHECK BAND ORDER AND GET BAND INDEX ;------------------------------------------------------------------------------- if (band ne 26) then begin ;- Get the list of band numbers for this SDS array att_info = hdf_sd_attinfo(sd_id, 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' endif else begin band_index = 0 endelse ;------------------------------------------------------------------------------- ;- READ IMAGE DATA ;------------------------------------------------------------------------------- ;- Get information about the image array varinfo = hdf_sd_varinfo(sd_id, sds_name) if (varinfo.name eq '') then $ message, 'Image array was not found: ' + sds_name npixels_across = varinfo.dims[0] npixels_along = varinfo.dims[1] ;- Set start and count values start = [0L, 0L, band_index] count = [npixels_across, npixels_along, 1L] if (band eq 26) then begin start = start[0 : 1] count = count[0 : 1] endif ;- Use AREA keyword if it was supplied if (n_elements(area) eq 4) then begin start[0] = (long(area[0]) > 0L) < (npixels_across - 1L) start[1] = (long(area[1]) > 0L) < (npixels_along - 1L) count[0] = (long(area[2]) > 1L) < (npixels_across - start[0]) count[1] = (long(area[3]) > 1L) < (npixels_along - start[1]) endif ;- Read the image array (hdf_sd_varread not used because of bug in IDL 5.1) var_id = hdf_sd_select(sd_id, hdf_sd_nametoindex(sd_id, sds_name)) hdf_sd_getdata, var_id, image, start=start, count=count hdf_sd_endaccess, var_id ;------------------------------------------------------------------------------- ;- READ IMAGE ATTRIBUTES ;------------------------------------------------------------------------------- ;- Read scale attribute att_info = hdf_sd_attinfo(sd_id, sds_name, scale_name) if (att_info.name eq '') then message, 'Attribute not found: ' + scale_name scale = att_info.data ;- Read offset attribute att_info = hdf_sd_attinfo(sd_id, sds_name, offset_name) if (att_info.name eq '') then message, 'Attribute not found: ' + offset_name offset = att_info.data ;- Read units attribute att_info = hdf_sd_attinfo(sd_id, sds_name, units_name) if (att_info.name eq '') then message, 'Attribute not found: ' + units_name units = att_info.data ;- Read valid range attribute valid_name = 'valid_range' att_info = hdf_sd_attinfo(sd_id, sds_name, valid_name) if (att_info.name eq '') then message, 'Attribute not found: ' + valid_name valid_range = att_info.data ;- Read latitude and longitude arrays if arg_present(latitude) then hdf_sd_varread, sd_id, 'Latitude', latitude if arg_present(longitude) then hdf_sd_varread, sd_id, 'Longitude', longitude ;------------------------------------------------------------------------------- ;- CLOSE THE FILE IN SDS MODE ;------------------------------------------------------------------------------- hdf_sd_end, sd_id ;------------------------------------------------------------------------------- ;- READ VDATA INFORMATION ;------------------------------------------------------------------------------- ;- Open the file in vdata mode hdfid = hdf_open(fileinfo.name) ;- Read scan start time (SDPTK TAI seconds) and mirror side vdataname = 'Level 1B Swath Metadata' hdf_vd_vdataread, hdfid, vdataname, 'EV Sector Start Time', scantime hdf_vd_vdataread, hdfid, vdataname, 'Mirror Side', mirror ;- Close the file in vdata mode hdf_close, hdfid ;- Resample and extract the scantime scantime = rebin(scantime, npixels_along, /sample) scantime = scantime[start[1]: start[1] + count[1] - 1] ;- Resample and extract the mirror side mirror = rebin(mirror, npixels_along, /sample) mirror = mirror[start[1]: start[1] + count[1] - 1] ;------------------------------------------------------------------------------- ;- PERFORM OPERATIONS IN SCALED INTEGER UNITS ;------------------------------------------------------------------------------- ;- Convert from unsigned short integer to long integer image = temporary(image) and 65535L valid_range = valid_range and 65535L ;- Get valid/invalid indexes and counts valid_check = (image ge valid_range[0]) and (image le valid_range[1]) invalid_index = where(valid_check eq 0B, invalid_count) if (arg_present(valid_index) or arg_present(valid_count)) then $ valid_index = where(valid_check eq 1B, valid_count) ;- If raw HDF values were requested, return to caller if keyword_set(raw) then begin units = 'Unsigned 16-bit integers' parameter = 'Raw HDF Values' if arg_present(range) then range = valid_range return endif ;- Replace saturated pixels in bands 1-5 with maximum scaled integer sat_count = 0 if (band le 5) then begin sat_check = (image eq 65533) or (image eq 65528) sat_index = where(sat_check eq 1B, sat_count) if (sat_count ge 1) then image[sat_index] = valid_range[1] endif ;------------------------------------------------------------------------------- ;- CONVERT IMAGE TO REQUESTED OUTPUT UNITS ;------------------------------------------------------------------------------- ;- Convert image to unscaled values image = scale[band_index] * (temporary(image) - offset[band_index]) valid_range = scale[band_index] * (valid_range - offset[band_index]) scale_factor = scale[band_index] ;- Convert radiance to brightness temperature for IR bands if keyword_set(temperature) then begin image = modis_bright(temporary(image), band, 1) valid_range = modis_bright(valid_range, band, 1) units = 'Kelvin' parameter = 'Brightness Temperature' endif ;- Set non-saturated invalid values to missing value if (invalid_count gt 0) then begin if (sat_count eq 0) then begin missing_index = invalid_index endif else begin missing_index = where(valid_check eq 0B and sat_check eq 0B, missing_count) endelse if (missing_index[0] ne -1L) then image[missing_index] = missing endif ;- Set non-finite values to missing value nonfinite_index = where(finite(image) ne 1, nonfinite_count) if (nonfinite_count gt 0) then image[temporary(nonfinite_index)] = missing ;- Get data range (min/max of non-missing image values) if arg_present(range) then begin nonmissing_index = where(image ne missing, nonmissing_count) if (nonmissing_count gt 0) then begin minval = min(image[temporary(nonmissing_index)], max=maxval) range = [minval, maxval] endif else begin range = [missing, missing] endelse endif END