pro misr_envi_level2brf,event ; This routine makes data from a MISR Level 2 data 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\Level2" AGPpath = "G:\data\misr\AGP" grid = "SubregParamsLnd" field = "LandBRF" bnames = ['Land BRF Blue','Land BRF Green','Land BRF Red','Land BRF NIR'] datafilepath = DIALOG_PICKFILE(PATH=DATApath,FILTER="*.hdf",TITLE="Select MISR Level 2 File") if (datafilepath eq '') then return datafile = FILE_BASENAME(datafilepath) ; note--extracting the path label in this manner depends on the MISR ; files retaining their original names; if the files have been renamed ; this code may need to be modified fparts = strsplit(datafile,'_',/EXTRACT) path = fparts[4] fid = EOS_GD_OPEN(datafilepath,/READ) ; get dimension information and projection parameters ; use the grid for the selected parameter so that appropriate ; information will be obtained gid = EOS_GD_ATTACH(fid,grid) status = EOS_GD_GRIDINFO(gid,xdimsize,ydimsize,upleft,lowright) status = EOS_GD_PROJINFO(gid,projcode,zonecode,spherecode,projparam) status = EOS_GD_FIELDINFO(gid,field,rank,dims,numtype,dimlist) case ydimsize of 8: pixsize = 70.4 * 1000 16: pixsize = 35.2 * 1000 32: pixsize = 17.6 * 1000 256: pixsize = 2.2 * 1000 512: pixsize = 1.1 * 1000 2048: pixsize = 275 else : endcase ; get the start and end blocks for the file status = read_global_attr(datafilepath,"Start_block",start_block) status = read_global_attr(datafilepath,"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(datafilepath,grid) ; 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 ; get the camera of interest label = "Choose a camera number from 1-9 (DF=1,CF=2,BF=3,AF=4,AN=5,AA=6,BA=7,CA=8,DA=9)" init_text = '5' get_text,init_text,camera,LABEL=label camera = camera - 1 ; get LandBRF Scale and Offset values status = read_grid_attr(datafilepath,grid,"Scale LandBRF",scale_landbrf) stauts = read_grid_attr(datafilepath,grid,"Offset LandBRF",offset_landbrf) status = read_grid_attr(datafilepath,grid,"Fill LandBRF",fill_landbrf) scale_landbrf = scale_landbrf[0] offset_landbrf = offset_landbrf[0] fill_landbrf = fill_landbrf[0] 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 for i=0,3 do begin start = [camera,i,0,0,start_block-1] ; read the data for one band edge = [1,1,dims[2],dims[3],num_blocks] status = EOS_GD_READFIELD(gid,field,bdata,START=start,EDGE=edge) bdata = reform(bdata) ; convert from raw values to scaled values, removing quality indicator ; in low order bit, multiplying by the scale factor, and adding the offset ; set bad values to 0 notgood = where(bdata ge fill_landbrf,count) print,"Not good count: ",count bdata = floor(bdata/2) bdata = bdata*scale_landbrf + offset_landbrf if (count gt 0) then begin bdata[notgood] = 0 endif ; 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_DETACH(gid) 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(datafilepath,'.hdf')+"_B"+strtrim(string(start_block),2)+"-"+$ strtrim(string(end_block),2) som_block_corner_meters,datafilepath,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