;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; SIMAP.4.pro ; ; ; Changes introduced with SIMAP.4: ; 1. Introduced special inter- and extrapolation of geolocations ; for MODIS L1B, including direct broadcast. Thus, the ; georeferencing (mapping) is substantially improved, ; especially in the case of 1-km L1B and internal geolocation arrays. ; 2. Introduced dynamic offsets,ia,ja, to better fill holes in the ; mapped images. ; 3. Made changes allowing certain AIRS data to be mapped as well, ; including the visible L1B, AIRVBRAD. ; SIMAP can now stitch that AIRS data. ; 4. The CONVERT_COORD is vectorized to speed up the process. If the ; geolocation arrays are too big, the code will switch to a FOR-cycle ; to avoid possible memory problems. ; ; simap.2.1. adds a function SCALEGEO to scale integer geolocations ; to real values, if local attributes like "Slope" and "Intercept" ; to the geolocation arrays are present. ; Such is the case with MODOCQC/MYDOCQC, and the new ; joint atmosphere product, MODATML2/MYDATML2. ; ; simap.2.1. also has a slight change to the common block, that ; makes more variables be seen from the batch file. This was ; done to allow more sofisticated maps generations where more ; than three layers can be superimposed, e.g. to combine ; true color visible image with some science parameter. ; ; SIMAP2 has the row/column subset setup moved to the batch file. ; Thus, the spatial subset is controlled from the batch file now. ; ; Since more and more users have IDL5.5 where the GIF license ; was not acquired by RSI, the code check the IDL version to determine ; what graphic interface to use when it comes to write out a one-layer ; image file. If the IDL version is 5.3 and below, it will use GIF. Otherwise ; it will use PNG. ; ; The code now uses pixel-by-pixel filling of the map space. ; Thus, at the edges of the mapped swath, it handles better the ; larger pixel size there. It is activated using the same flag, ; intrpl=1. This approach is slow for large arrays, but gives ; more correct results. A new module is added to perform the ; pixel-by-pixel filling, MAPIT. ; ; This IDL code was developed at the Goddard DAAC Modis ; Data Support Team. Contact information: ; Andrey Savtchenko ; (asavtche@daac.gsfc.nasa.gov). ; Please, address all inquiries to ; help$daac.gsfc.nasa.gov. ; ; We will appreciate your credits and critics any time you use our code. ; The program is intended for MODIS swath, Level 1B and 2, data. ; ; The program is designed to be run from an IDL batch file. See for ; batch file examples, and how to run it, in the documentation, ; http://daac.gsfc.nasa.gov/MODIS/simap/ ; ; Most important features are: mapping and stitching of MODIS and AIRS L1B and L2 ; granules, saving mapped data in physical units, correction for ; invalid geolocations, and an option to subset a portion of the SDS and ; map it. COMMON SHARE1,xo,yo,NXM,NYM,NX,NY,inpath,outpath,$ map_limit,mappos,True_col,dname,zbuf,wid,$ center_lon,center_lat,mnlat,mxlat,sel,debug,$ jpgname,flcnt,xstart,xend,ystart,yend,$ xtitle,unit,zlog,ctab,proj,gallery,rgb,r,g,b,L1,$ latmin,latmax,lonmin,lonmax,intrpl,$ extgeo,sdsprint,automap,$ autointen,intens,slice,backgnd PRO GETDATA, sd_id, SDSNo, name, array, dims, chanid, fillvalue,slice chanid=0 sds_id = HDF_SD_SELECT(SD_ID,SDSNo) HDF_SD_GETDATA, SDS_ID, array HDF_SD_GETINFO, sds_id, LABEL=l, UNIT=u, FORMAT=f, $ COORDSYS=c,NDIMS=ndims, DIMS=dims, TYPE=ty,NAME=name ;find the id of the band_names !QUIET=1 dindex = HDF_SD_ATTRFIND(sds_id, 'band_names') !QUIET=0 ;---This is the AIRS VIS module------- ;---It will reformat the sds from 5D into 3D IF TOTAL(name EQ ['radiances','PrelimCldMapVis','PrelimNDVI']) AND $ ( (SIZE(array))[0] EQ 5 ) THEN BEGIN array=TRANSPOSE(array,[0,3,1,4,2]) dims=[dims[0]*dims[3],dims[1]*dims[4],dims[2]] array=REFORM(array,dims[0],dims[1],dims[2],/OVERWRITE) ndims=3 ENDIF ;if MODIS L1B Earth View SD is selected: IF dindex GT 0 THEN BEGIN HDF_SD_ATTRINFO, sds_id, dindex, DATA=chanid chanid=STR_SEP(chanid, ",") ENDIF ELSE BEGIN ;if cannot find band names, ;then simply use 1,2,... IF ndims GT 2 THEN $ chanid=STRING(INDGEN(dims[slice-1])+1) chanid=STRCOMPRESS(chanid,/REMOVE_ALL) ENDELSE ;---find the id of the fill value--- !QUIET=1 dindex = HDF_SD_ATTRFIND(sds_id, '_FillValue') IF dindex GT 0 THEN $ HDF_SD_ATTRINFO, sds_id, dindex, DATA=fillvalue !QUIET=0 ;------------------------------------ RETURN END FUNCTION SCALEGEO,geosd_id,SDSNo,geo ;The following will scale geolocations if they ;have "Slope" and "Intercept" information ;which is the case with MODOCQC/MYDOCQC ;and MODATML2/MYDATML2 sds_id = HDF_SD_SELECT(geosd_id,SDSNo) Slope=0 Intercept=0 !QUIET=1 sindex = HDF_SD_ATTRFIND(sds_id, 'Slope') !QUIET=0 IF sindex GE 0 THEN $ HDF_SD_ATTRINFO, sds_id,sindex, DATA=Slope !QUIET=1 sindex = HDF_SD_ATTRFIND(sds_id, 'scale_factor') !QUIET=0 IF sindex GE 0 THEN $ HDF_SD_ATTRINFO, sds_id,sindex, DATA=Slope !QUIET=1 index = HDF_SD_ATTRFIND(sds_id, 'Intercept') !QUIET=0 IF index GE 0 THEN $ HDF_SD_ATTRINFO, sds_id,index, DATA=Intercept !QUIET=1 index = HDF_SD_ATTRFIND(sds_id, 'add_offset') !QUIET=0 IF index GE 0 THEN $ HDF_SD_ATTRINFO, sds_id,index, DATA=Intercept IF NOT(((Slope[0] EQ 1) and (Intercept[0] eq 0)) or (Slope[0] eq 0)) THEN $ geo=Slope[0]*geo+Intercept[0] RETURN,geo END ;This function is intended for AIRS visible ;bands geolocations. Presently, to make things fast and ;simple, it will use the geolocations of the first selected ;layer of "radiances". These geolocations are interpolated here ;and reformated to a 2D array. FUNCTION AIRS_VIS_GEO,input,sel,dims nch=1 tmp=FLTARR(2,2) geo=FLTARR(dims[2]*8,dims[3]*9,nch) input=TRANSPOSE(input,[1,2,3,0]) input=SHIFT(input,2,0,0,0) ;test for International Date Line test=TOTAL(input GT 179.5)*TOTAL(input LT -179.5) ;if not crossing Inter. Date Line: IF(test EQ 0) THEN BEGIN FOR k=0,nch-1 DO BEGIN FOR i=0,dims[2]-1 DO BEGIN FOR j=0,dims[3]-1 DO BEGIN tmp[*,*]=input[*,i,j,sel[0].chans[k]-1] geo[i*8:(i+1)*8-1,j*9:(j+1)*9-1,k]=CONGRID(tmp,8,9,/INTERP,/MINUS_ONE) ENDFOR ENDFOR ENDFOR ENDIF ;if crossing Inter. Date Line: IF(test NE 0) THEN BEGIN FOR k=0,nch-1 DO BEGIN FOR i=0,dims[2]-1 DO BEGIN FOR j=0,dims[3]-1 DO BEGIN tmp[*,*]=input[*,i,j,sel[0].chans[k]-1] xcoor=COS(tmp*!DTOR) ycoor=SIN(tmp*!DTOR) xcoor=CONGRID(xcoor,8,9,/interp,/minus_one) ycoor=CONGRID(ycoor,8,9,/interp,/minus_one) lons=ATAN(ycoor,xcoor)*!RADEG geo[i*8:(i+1)*8-1,j*9:(j+1)*9-1,k]=lons ENDFOR ENDFOR ENDFOR ENDIF RETURN,geo END PRO SHOWIMAGE,rrange,grange,brange COMMON SHARE1 SET_PLOT,'Z' DEVICE,SET_RESOLUTION=[NX,NY],Z_BUFFERING=0,$ SET_CHARACTER_SIZE=[6,9] CASE proj OF '1': BEGIN MAP_SET,center_lat,center_lon,0,/CYLINDRICAL,$ /NOBORDER,XMARGIN=0,YMARGIN=0,$ LIMIT=map_limit,POSITION=mappos END '2': BEGIN MAP_SET,center_lat,center_lon,0,$ /stereo,XMARGIN=0,YMARGIN=0,$ /NOBORDER,LIMIT=map_limit,POSITION=mappos END ENDCASE temp=WHERE(rgb[*,*,0] EQ -999. ,count) IF(True_col EQ 3) THEN BEGIN ;scaling to bytes: IF(zlog) THEN BEGIN ;The following manages for negative values. ;Also, it does not allow for more than 8 ;orders of magnitude for the log plot. ;That must be more than enough, having in mind ;the S/N ratio of the data. red1=rgb[*,*,0]>1.e-20 green1=rgb[*,*,1]>1.e-20 blue1=rgb[*,*,2]>1.e-20 rrange[0]=rrange[0]>rrange[1]/1.e8 grange[0]=grange[0]>grange[1]/1.e8 brange[0]=brange[0]>brange[1]/1.e8 red1=BYTSCL(ALOG10(red1),$ MIN=ALOG10(rrange[0]),MAX=ALOG10(rrange[1]),TOP=252) green1=BYTSCL(ALOG10(green1),$ MIN=ALOG10(grange[0]),MAX=ALOG10(grange[1]),TOP=252) blue1=BYTSCL(ALOG10(blue1),$ MIN=ALOG10(brange[0]),MAX=ALOG10(brange[1]),TOP=252) ENDIF ELSE BEGIN red1=BYTSCL(rgb[*,*,0],MIN=rrange[0],MAX=rrange[1],TOP=252) green1=BYTSCL(rgb[*,*,1],MIN=grange[0],MAX=grange[1],TOP=252) blue1=BYTSCL(rgb[*,*,2],MIN=brange[0],MAX=brange[1],TOP=252) ENDELSE ;adding the background: IF(count GT 0) THEN BEGIN red1[temp]=backgnd green1[temp]=backgnd blue1[temp]=backgnd ENDIF L1=[[[red1]],[[green1]],[[blue1]]] ;L1 is exported ;in case combined L1-L2 images ;are desired ENDIF ELSE BEGIN ;begin if not TrueColor: IF(zlog) THEN BEGIN red1=rgb>1.e-20 rrange[0]=rrange[0]>rrange[1]/1.e8 ;in case autointens gave red1=BYTSCL(ALOG10(red1),$ ;negatives MIN=ALOG10(rrange[0]),MAX=ALOG10(rrange[1]),TOP=252) ENDIF ELSE BEGIN red1=BYTSCL(rgb,MIN=rrange[0],MAX=rrange[1],TOP=252) ENDELSE IF(count GT 0) THEN red1[temp]=backgnd L1=red1 ENDELSE ;background outside map ;must be different than lines and fills ERASE,255 ;these are color indexes of fills and lines filcol=254 lccol=0 ;print,"map continents...." IF(flcnt) THEN $ MAP_CONTINENTS,/HIRES,FILL_CONTINENTS=1,COLOR=filcol MAP_CONTINENTS,/COUNTRIES,/COASTS,COLOR=lccol MAP_GRID,/BOX_AXES,CHARSIZE=1.5,LABEL=1, $ LONLAB=map_limit[2]-.31,LATLAB=(map_limit[3]-.3),color=lccol ; red1[100:750,5:35]=0 ; XYOUTS,/device,xo+100,yo+10,"!3!17MODIS Data Support Team at GES DAAC!X",$ ; CHARSIZE=2,COLOR=filcol lines=TVRD() line=WHERE(lines EQ lccol,lcount) fill=WHERE(lines EQ filcol,fcount) zbuf=MAKE_ARRAY(NX,NY,3,/BYTE,VALUE=255) IF (True_col EQ 3) THEN BEGIN TV,TEMPORARY(red1),xo,yo zbuf[*,*,0]=TVRD() TV,TEMPORARY(green1),xo,yo zbuf[*,*,1]=TVRD() TV,TEMPORARY(blue1),xo,yo zbuf[*,*,2]=TVRD() kk=2 ENDIF ELSE BEGIN TV,TEMPORARY(red1),xo,yo COLSCALE,rrange zbuf=TVRD() kk=0 ENDELSE FOR k=0,kk DO BEGIN trash=zbuf[*,*,k] IF(lcount GT 0) THEN trash[line]=lccol IF(fcount GT 0) THEN trash[fill]=filcol zbuf[*,*,k]=trash ENDFOR ;load color palette and ensure: ;0 = black ;253 = blue-ish (backgnd) ;254 = fill of the continents ;255 = white loadct,ctab tvlct,r,g,b,/get r[[0,253,254,255]]=[0,180,180,255] g[[0,253,254,255]]=[0,246,247,255] b[[0,253,254,255]]=[0,240,121,255] tvlct,r,g,b RETURN END ;----This module will generate color scale legend----- PRO COLSCALE,rrange COMMON SHARE1 pos=[mappos[0],0.035,0.45,.045] fill=BYTE(BINDGEN(253)#$ (BYTARR(50)+1b)) nclx=(pos[2]-pos[0])*NX ncly=(pos[3]-pos[1])*NY fill=CONGRID(fill,nclx,ncly) TV,fill,pos[0],pos[1],/NORMAL IF(zlog) THEN BEGIN ;--log scale-- xticks=ALOG10(rrange[1]/rrange[0]) PLOT,rrange,[0,1],/NODATA,/NOERASE,$ POS=pos, $ XSTYLE=1,YSTYLE=4,COLOR=0,/XLOG, $ XTICKS=xticks,XMINOR=1,TICKLEN=1,$ YTICKS=1,YMINOR=1,$ CHARSIZE=1.5,XTICKFORMAT='(E10.0)' ENDIF ELSE BEGIN ;--linear scale-- PLOT,rrange,[0,1],/NODATA,/NOERASE,$ POS=pos, $ XSTYLE=1,YSTYLE=4,COLOR=0, $ XTICKS=5,XMINOR=1,TICKLEN=1,$ YTICKS=1,YMINOR=1,$ CHARSIZE=1.5,XTICKFORMAT='(F10.1)' ENDELSE IF(N_ELEMENTS(xtitle) EQ 0) THEN xtitle="" IF(N_ELEMENTS(unit) EQ 0) THEN unit="" title=" "+xtitle+", "+unit XYOUTS,pos[2],pos[1],title,/NORMAL,$ CHARSIZE=1.5,COLOR=0 RETURN END PRO WRITEJPEG, namejpg COMMON SHARE1 temp=SIZE(zbuf) IF (temp[0] LT 3) THEN kk=0 ELSE kk=3 WRITE_JPEG,namejpg,zbuf,QUALITY=100,TRUE=kk RETURN END ; This module is activated when intrpl=1 ; It is a way of mapping, i.e. moving of pixels ; from the satellite raster to the particular ; map projection. It handles swath edges very well, ; but slows down the processing. PRO MAPIT,map,xv,yv,temp,x_size,y_size,ia,ja COMMON SHARE1 FOR j=0,y_size-2 DO BEGIN FOR i=0,x_size-2 DO BEGIN x1=xv[i,j] x2=(xv[i+1,j+ja]-1)>0 y1=yv[i,j] y2=yv[i+ia,j+1] IF(ABS(x2-x1) LT 30) THEN $ map[x1x1,y1y1] = temp[i,j] ENDFOR ENDFOR RETURN END ;;;;;;;;;;;;;;;;;;;;;;;;M A I N M O D U L E;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; ; ; ; PRO simap, hdfname COMMON SHARE1 ;============================================================ shortname="" long_name="" units="" dname=!D.NAME ;default graphic device temp="" ;a temporal, empty variable geosd_id=-1 ;default sd_id of geolocation file fillvalue=-999. ;in case _FillValue cannot be found ;The code will automatically distinguish between L1B and L2 products ;and will switch the scaling equation. Terra and Aqua are considered. L1B=['MOD021KM','MOD02HKM','MOD02QKM','GSUB1','MOD02SSH','MYD02SSH',$ 'MYD021KM','MYD02HKM','MYD02QKM','RMT021KM','RMT02HKM','RMT02QKM',$ 'RMA021KM','RMA02HKM','RMA02QKM'] OL2=['MODOCL2','MODOCL2A','MODOCL2B','MOD28L2','MODOCQC','MOD28QC',$ 'MYDOCL2','MYDOCL2A','MYDOCL2B','MYD28L2','MYDOCQC','MYD28QC'] AL2=['MOD04_L2','MOD05_L2','MOD06_L2','MOD07_L2','MOD35_L2','MODATML2',$ 'MYD04_L2','MYD05_L2','MYD06_L2','MYD07_L2','MYD35_L2','MYDATML2'] ;============================================================= ;check for backslash at the end of the paths: IF STRMID(inpath,STRLEN(inpath)-1) NE "/" THEN $ inpath=inpath+"/" IF STRMID(outpath,STRLEN(outpath)-1) NE "/" THEN $ outpath=outpath+"/" ;this is the log file. ;will log there if debug=4 CLOSE,/ALL logname=outpath+'log.txt' temp=FINDFILE(logname,COUNT=count) IF(debug GT 0) THEN BEGIN debug=4 IF (not count) THEN OPENW,debug,logname $ ELSE OPENU,/APPEND,debug,logname ENDIF ;check options IF(TOTAL(True_col EQ [0,1,3]) NE 1) THEN BEGIN PRINTF,debug,"True_col must be either 0, 1, or 3" CLOSE,/ALL RETURN ENDIF IF(gallery and True_col EQ 0) THEN BEGIN PRINTF,debug,"gallery=1 and True_col=0 is not acceptable." PRINTF,debug," If you wish to make galleries, "+$ "set True_col=1 or 3" RETURN ENDIF IF(gallery EQ 0) THEN BEGIN automap=0 ENDIF ;check if the input argument is a single file name ;or the text file containing file names temp='' IF (hdfname EQ 'filelist') THEN BEGIN OPENR,5,inpath+hdfname filename='' WHILE NOT EOF(5) DO BEGIN READF,5,temp filename=filename+temp+' ' ENDWHILE CLOSE,5 filename=STRTRIM(filename,2) FileName=STRSPLIT(filename,' ',/EXTRACT) numfls=SIZE(FileName,/N_ELEMENTS)-1 ENDIF ELSE BEGIN filename=[hdfname] numfls=0 ENDELSE ;***********Here we start to cycle over HDF files**************** FOR fn=0,numfls DO BEGIN ;Inquire for geolocation file, if external ;geolocations are requested: IF (extgeo) THEN BEGIN char1=STRPOS(FileName[fn],".") IF(STRMID(FileName[fn],0,2) EQ "MO") THEN prefix="MOD03" $ ELSE prefix="MYD03" geofname=inpath+prefix+STRMID(FileName[fn],char1,15)+"*.hdf" geofname=FINDFILE(geofname,COUNT=gcount) IF gcount GT 0 THEN BEGIN geofname=geofname[0] ENDIF ELSE BEGIN PRINTF,debug,"Cannot find geolocation file",geofname CLOSE,/ALL RETURN ENDELSE ENDIF ;open HDF for READ, take file ID: FileHandle=HDF_OPEN(inpath+FileName[fn],/READ) ;check if open is successfull: IF FileHandle LT 0 THEN BEGIN PRINTF,debug,"Error open ", (inpath+FileName[fn]) CLOSE,/ALL RETURN ENDIF ;Extract the common part to use in the ; output files. ;Later, will add parameter name and ;"bin","jpg",or "png" pos=1000 FOR l=0,2 DO pos=STRPOS(FileName[fn],".",pos-1,/REVERSE_SEARCH) gname=STRMID(FileName[fn],0,pos)+"." PRINTF,debug PRINTF,debug,"=======================" PRINTF,debug,SYSTIME(0) PRINTF,debug,"Working on ",inpath+FileName[fn] ;Initialize Science Data (SD) interface, ;take the SD ID and assign to sd_id. ;It will point to the FileName[fn]. sd_id=HDF_SD_START(inpath+FileName[fn],/READ) ;Take the number of data sets (SDS) ; and assign to the variable datasets. HDF_SD_FILEINFO, sd_id, datasets, attributes CASE sdsprint OF '1': BEGIN longprint=0 sdsprint=1 END '2': BEGIN longprint=1 sdsprint=1 END '0': BEGIN longprint=0 sdsprint=0 END '4': goto,ex ELSE: goto,ex ENDCASE If (sdsprint) THEN BEGIN ;The following will list all SDS found in the file. ; There are "datasets" number of SDS in the file. PRINTF,debug," The following SDS's [and attributes if "+ $ "long print is active] are available:" ;=================start cycle SDS's to print on STDOUT device============= FOR I=0,datasets-1 DO begin sds_id=HDF_SD_SELECT(sd_id,I) HDF_SD_GETINFO,sds_id,DIMS=dims,NDIMS=ndims,$ NAME=name,NATTS=natts PRINTF,debug,FORMAT='(I0,".",3A0,4(I0,:,"x"))', $ I, "Short name: ", name , ", size: " , dims IF longprint NE 0 THEN BEGIN ;_____start cycle attributes for every SDS____ FOR J=0,natts-1 DO BEGIN HDF_SD_ATTRINFO,sds_id,J,NAME=name,DATA=attr_dat PRINTF,debug,FORMAT='(A0,": ",5(A0))',name,attr_dat attr_dat="" name="" ENDFOR ;______end cycle attributes_______ ENDIF ;end if longprint NE 0 ENDFOR ;============end cycle SDS's====================== ENDIF ;end if(sdsprint NE 0) gindex= HDF_SD_ATTRFIND(sd_id, "Number of Day mode scans") IF gindex NE -1 THEN BEGIN HDF_SD_ATTRINFO,sd_id, gindex, NAME=name, DATA=ndms PRINTF,debug,name,ndms ENDIF gindex= HDF_SD_ATTRFIND(sd_id, "Number of Night mode scans") IF gindex NE -1 THEN BEGIN HDF_SD_ATTRINFO,sd_id, gindex, NAME=name, DATA=nnms PRINTF,debug,name,nnms nnms=nnms[0] ndms=ndms[0] IF(nnms EQ 0) THEN DAYNIGHTFLAG='Day' IF(ndms EQ 0) THEN DAYNIGHTFLAG='Night' IF((ndms GT 0) AND (nnms GT 0)) THEN $ DAYNIGHTFLAG='Both' ENDIF gindex= HDF_SD_ATTRFIND(sd_id, "Number of Scans") IF gindex NE -1 THEN BEGIN HDF_SD_ATTRINFO,sd_id,gindex, NAME=name, DATA=nscan nscan=nscan[0] PRINTF,debug,name,nscan ENDIF gindex= HDF_SD_ATTRFIND(sd_id, "Max Earth View Frames") IF gindex NE -1 THEN BEGIN HDF_SD_ATTRINFO,sd_id, gindex, NAME=name, DATA=nsamp PRINTF,debug,name,nsamp ENDIF gindex= HDF_SD_ATTRFIND(sd_id, "Incomplete Scans") IF gindex NE -1 THEN BEGIN HDF_SD_ATTRINFO,sd_id, gindex, NAME=name, DATA=iscns PRINTF,debug,name,iscns ENDIF gindex= HDF_SD_ATTRFIND(sd_id, "CoreMetadata.0") IF(gindex LT 0) THEN gindex= HDF_SD_ATTRFIND(sd_id, "coremetadata") IF gindex NE -1 THEN BEGIN PRINTF,debug,"Query CoreMetadata for DAYNIGHTFLAG:" HDF_SD_ATTRINFO,sd_id, gindex, DATA=CoreMetadata temp=STRPOS(CoreMetadata,'DAYNIGHTFLAG') temp=STRPOS(CoreMetadata,'"',temp+15) temp=STRMID(CoreMetadata,temp+1,3) CASE temp OF 'Day': DAYNIGHTFLAG='Day' 'Nig': DAYNIGHTFLAG='Night' 'Bot': DAYNIGHTFLAG='Both' ELSE: PRINTF,debug,"DAYNIGHTFLAG not in CoreMetadata, or format is different" ENDCASE ;extract SHORTNAME: SHORTNAME='SHORTNAME' temp=STRPOS(CoreMetadata,SHORTNAME) temp=STRPOS(CoreMetadata,'"',temp+15) length=STRPOS(CoreMetadata,'"',temp+1)-temp-1 SHORTNAME=STRMID(CoreMetadata,temp+1,length) ;extract Platform Name: ;PLATFORM='ASSOCIATEDPLATFORMSHORTNAME' ;temp=STRPOS(CoreMetadata,PLATFORM) ;temp=STRPOS(CoreMetadata,"VALUE",temp+1) ;temp=STRPOS(CoreMetadata,'"',temp+1) ;length=STRPOS(CoreMetadata,'"',temp+1)-temp-1 ;PLATFORM=STRMID(CoreMetadata,temp+1,length) ENDIF ELSE BEGIN PRINTF,debug,"CoreMetadata.0 does not exist in this file!" PRINTF,debug,"Scaling equation cannot be chosen and data "+$ "will not be scaled to physical units." ENDELSE IF(N_ELEMENTS(DAYNIGHTFLAG) EQ 0) THEN BEGIN PRINTF,debug,'Cannot find DAYNIGHTFLAG. Arbitrary set to "Day"' DAYNIGHTFLAG="Day" ENDIF PRINTF,debug,"DAYNIGHTFLAG=",DAYNIGHTFLAG ;The index of the last selection with non-empty name SS=MAX(WHERE(sel.name NE '') ) ;Calculate number of channels for each ;selection, list selected SDS names and ;channels: PRINTF,debug,"Channel selection:" FOR I=0,SS DO BEGIN sel[i].numch=TOTAL(sel[i].chans NE 0) IF (sel[i].numch GT 0) THEN BEGIN sel[i].chans=STRTRIM(sel[i].chans,2) temp=sel[i].numch-1 PRINTF,debug,sel[i].name,", channel(s):",sel[i].chans[0:temp] ENDIF ENDFOR ;Total number of all selected channels: totselch=TOTAL(sel[0:SS].numch) temp=1 CASE temp OF '1': temp=1 '2': temp=0 '3': goto,ex ELSE: goto,ex ENDCASE ;-------------------------------------------------------------- ;If arrays in the selection are with different precision, ;this module will choose the highest one and it will ;be used to define the type of the array "image" ; type=BYTARR(SS+1) FOR I=0,SS DO BEGIN ;find the SDS offset number SDSNo = HDF_SD_NAMETOINDEX(sd_id, sel[i].name) IF (SDSNo NE -1) THEN BEGIN sds_id = HDF_SD_SELECT(SD_ID,SDSNo) HDF_SD_GETINFO, sds_id, TYPE=ty CASE ty OF "BYTE":type[i]=1 "INT":type[i]=2 "UINT":type[i]=12 "LONG":type[i]=3 "FLOAT":type[i]=4 "DOUBLE":type[i]=5 ELSE:type[i]=4 ENDCASE ENDIF ENDFOR type=MAX(type) ; ;-----------------------------------------------; ;=====================START TO RETRIEVE DATA=============; I0=0 FOR I=0,SS DO BEGIN slope=[1] ;default values in case these attributes offset=[0] ;are not found name=STRTRIM(sel[i].name,2) ;find the SDS offset number SDSNo = HDF_SD_NAMETOINDEX(sd_id, name) IF (SDSNo NE -1) THEN BEGIN PRINTF,debug,"Retrieve:",name GETDATA, sd_id, SDSNo, name, array,$ dims, chanid, fillvalue,slice ;---------------------------------------------------- ;For each SDS name, find attributes ; ;that scale data to physical units ; ; !QUIET=1 SDSNo = HDF_SD_NAMETOINDEX(sd_id, name) sds_id = HDF_SD_SELECT(SD_ID,SDSNo) ;for L1B: aindex = HDF_SD_ATTRFIND(sds_id, 'radiance_scales') IF aindex GE 0 THEN $ HDF_SD_ATTRINFO, sds_id, aindex, DATA=Slope aindex = HDF_SD_ATTRFIND(sds_id, 'radiance_offsets') IF aindex GE 0 THEN $ HDF_SD_ATTRINFO, sds_id, aindex, DATA=offset ;for Oceans aindex = HDF_SD_ATTRFIND(sds_id, 'Slope') IF aindex GE 0 THEN $ HDF_SD_ATTRINFO, sds_id, aindex, DATA=Slope aindex = HDF_SD_ATTRFIND(sds_id, 'Intercept') IF aindex GE 0 THEN $ HDF_SD_ATTRINFO, sds_id, aindex, DATA=offset ;for Atmosphere aindex = HDF_SD_ATTRFIND(sds_id, 'scale_factor') IF aindex GE 0 THEN $ HDF_SD_ATTRINFO, sds_id, aindex, DATA=Slope aindex = HDF_SD_ATTRFIND(sds_id, 'add_offset') IF aindex GE 0 THEN $ HDF_SD_ATTRINFO, sds_id, aindex, DATA=offset aindex = HDF_SD_ATTRFIND(sds_id, 'units') IF aindex GE 0 THEN $ HDF_SD_ATTRINFO, sds_id, aindex, DATA=units aindex = HDF_SD_ATTRFIND(sds_id, 'unit') IF aindex GE 0 THEN $ HDF_SD_ATTRINFO, sds_id, aindex, DATA=units aindex = HDF_SD_ATTRFIND(sds_id, 'Units') IF aindex GE 0 THEN $ HDF_SD_ATTRINFO, sds_id, aindex, DATA=units aindex = HDF_SD_ATTRFIND(sds_id, 'radiance_units') IF aindex GE 0 THEN $ HDF_SD_ATTRINFO, sds_id, aindex, DATA=units !QUIET=0 ;------------------------------------------------------ IF N_ELEMENTS(dims) GT 2 THEN $ dims=dims[WHERE(dims[slice-1] NE dims)] x_size=dims[0] y_size=dims[1] ;initialize arrays to keep all selected data ;and scaling coefficients: IF I0 EQ 0 THEN BEGIN image=MAKE_ARRAY(x_size,y_size,totselch,TYPE=type) pname=STRARR(totselch) radsc=FLTARR(totselch) radof=FLTARR(totselch) unita=STRARR(totselch) ENDIF ;sanity check: IF ((SIZE(array,/N_DIMENSIONS) EQ 2) AND $ (sel[i].numch GT 1)) THEN BEGIN PRINTF,debug,"Check sel[",i,"] - you selected"+$ " more channels than available" RETURN ENDIF IF SIZE(array,/N_DIMENSIONS) GT 2 THEN BEGIN ;cycle over channels: FOR k=0,sel[i].numch-1 DO BEGIN I2=WHERE(sel[i].chans[k] EQ chanid, Count) IF (Count EQ 0) THEN BEGIN PRINTF,debug,"Channel ",sel[i].chans[k],$ " is not found in this SDS." GOTO,ex ENDIF ;This segment will retrieve selected channels ;and corresponding scaling coefficients CASE slice OF 1: image[*,*,I0]=array[I2,*,*] 2: image[*,*,I0]=array[*,I2,*] 3: image[*,*,I0]=array[*,*,I2] ENDCASE radsc[I0]=slope[I2] radof[I0]=offset[I2] pname[I0]=sel[i].name+chanid[I2] unita[I0]=units I0=I0+1 ENDFOR ;end k-cycle (channels) ENDIF ; end if (size of array) GT 2 ;--if SDS is 2D only:-- IF SIZE(array,/N_DIMENSIONS) EQ 2 THEN BEGIN image[*,*,I0]=array radsc[I0]=slope[0] radof[I0]=offset[0] pname[I0]=sel[i].name unita[I0]=units I0=I0+1 ENDIF ;---------------------- ENDIF ELSE BEGIN ;--------end if (SDSNo NE -1)--------------- PRINTF,debug,"!!!" PRINTF,debug,"CANNOT RETRIEVE ",name," !!!" PRINTF,debug,"Review names in this selection." RETURN ENDELSE ENDFOR ; end SS-cycle (selections) ;################THIS IS A PLACE TO PUT YOUR CODE####### ; ;You can put your code here. All data, as it is in the HDF, and as you ;selected in the batch file (the structure "sel") is in the array ; "image". If you selected more than one channel, "image" is three ;dimensional. It is not mapped, and not converted to physical ;units at this point. If, e.g. image[*,*0] is your data, and ;image[*,*,1] your quality flags, you can make quality filtering ;at this point. See documentation, "TRICKS AND TIPS" ;rgb=FLOAT(image) ;rgb[*,*,0]=radsc[0]*image[*,*,0]+radof[0] ;rgb[*,*,1]=radsc[1]*image[*,*,1]+radof[1] ;################################################################## ;Some sanity checks: IF(I0 EQ 0) THEN BEGIN PRINTF,debug,"No SDS retrieved. Exit." GOTO,ex ENDIF IF((True_col EQ 3) AND $ (TOTAL(sel.numch) LT 3)) THEN BEGIN PRINTF,debug,"Requested is true color image "+$ "but selected data layers are less than 3" PRINTF,debug,"Select 3 data layers, "+$ "or set True_col=1 or 0" RETURN ENDIF ;--------------FIND GEOLOCATION FIELDS------------------------------ IF(extgeo) THEN BEGIN GeoFileHandle=HDF_OPEN(geofname,/READ) ;check if open is successfull: IF GeoFileHandle LT 0 THEN BEGIN PRINTF,debug,"Error open geolocation file" GOTO, ex ENDIF geosd_id=HDF_SD_START(geofname,/READ) ENDIF ELSE geosd_id=sd_id ;----Get Latitude---- name=['cornerlats','Latitude','latitude'] ;find the valid name and SDS offset number FOR i=0,N_ELEMENTS(name)-1 DO BEGIN SDSNo = HDF_SD_NAMETOINDEX(geosd_id, name[i]) IF (SDSNo GE 0) THEN BEGIN name=name[i] & BREAK ENDIF ENDFOR IF (SDSNo LT 0) THEN BEGIN PRINTF,debug,"Cannot find either ",name GOTO,ex ENDIF GETDATA, geosd_id, SDSNo, name, array, dims, chanid, fillvalue,slice lats=TEMPORARY(array) lats=SCALEGEO(geosd_id,SDSNo,lats) IF(name EQ "cornerlats") THEN lats=AIRS_VIS_GEO(lats,sel,dims) latfill=fillvalue ;----Get Longitude---- name=['cornerlons','Longitude','longitude'] ;find the SDS offset number FOR i=0,N_ELEMENTS(name)-1 DO BEGIN SDSNo = HDF_SD_NAMETOINDEX(geosd_id, name[i]) IF (SDSNo GE 0) THEN BEGIN name=name[i] & BREAK ENDIF ENDFOR IF (SDSNo LT 0) THEN BEGIN PRINTF,debug,"Cannot find either",name GOTO,ex ENDIF GETDATA, geosd_id, SDSNo, name, array, dims, chanid, fillvalue,slice lons=TEMPORARY(array) lons=SCALEGEO(geosd_id,SDSNo,lons) IF(name EQ "cornerlons") THEN lons=AIRS_VIS_GEO(lons,sel,dims) lonfill=fillvalue dims[0]=(SIZE(lons))[1] dims[1]=(SIZE(lons))[2] xscl=FLOAT(x_size)/dims[0] ;will use this later to yscl=FLOAT(y_size)/dims[1] ; resize geolocations IF(xscl GE 1) THEN xscl=ROUND(xscl) IF(yscl GE 1) THEN yscl=ROUND(yscl) ;-----------------end find and load geolocation fields---------- ;---Here we inquire on invalid geolocations and try to fix them---- IF TOTAL(ABS(lats) GT 90.) THEN BEGIN PRINT,"Invalid latitudes detected. Will attempt to fix." good_pxl=(ABS(lats) LT 90.) IF(TOTAL(good_pxl) EQ 0) THEN BEGIN PRINT,debug,"All lats are bad or missing!!!" GOTO,ex ENDIF iwant=INDGEN(dims[1]) FOR i=0,dims[0]-1 DO BEGIN mask=WHERE(good_pxl[i,*] GT 0,count) IF count GT 1 THEN BEGIN lats[i,*]=interpol(REFORM(lats[i,mask]),mask,iwant) ENDIF ENDFOR ENDIF ; endif total(abs(lats) gt 90) IF TOTAL(ABS(lons) GT 180.) THEN BEGIN PRINT,"Invalid longitudes detected. Will attempt to fix." good_pxl=(ABS(lons) LT 180.) IF(TOTAL(good_pxl) EQ 0) THEN BEGIN PRINT,debug,"All lons are bad or missing!!!" GOTO,ex ENDIF iwant=INDGEN(dims[1]) FOR i=0,dims[0]-1 DO BEGIN mask=WHERE(good_pxl[i,*] GT 0,count) IF count GT 1 THEN BEGIN lons[i,*]=interpol(REFORM(lons[i,mask]),mask,iwant) ENDIF ENDFOR ENDIF ; endif total(ABS(lons) gt 180.) ;--------------- end fixing invalid geolocations ------------------- ;--------SPATIAL SUBSETTING---------------------------------- ;x1,x2 and y1,y2 define column and row subset, respectively. ;"end" and "start" parameters come from the batch file. IF (xend EQ "") THEN x2=x_size-1 ELSE x2=xend < (x_size-1) IF (xstart EQ "") THEN x1=0 ELSE x1=xstart < (x_size-1) IF (yend EQ "") THEN y2=y_size-1 ELSE y2=yend < (y_size-1) IF (ystart EQ "") THEN y1=0 ELSE y1=ystart < (y_size-1) IF (y1 EQ y2) OR (x1 EQ x2) THEN BEGIN PRINTF,debug, "Subset indices larger than array size,",$ " or start=end." ENDIF ;---------------------------------------------------------------- ;The following block of code ensures that when a subset ;is requested, it is taken exactly at swath boundary (the ;along-track dimension, y) and exactly at a pixel where ;there is a geolocation (across track dimension,x) IF TOTAL(SHORTNAME EQ ["MOD021KM","MOD02HKM","MOD02QKM",$ "MYD021KM","MYD02HKM","MYD02QKM","GSUB1",$ "RMA021KM","RMA02HKM","RMA02QKM",$ "RMT021KM","RMT02HKM","RMT02QKM"]) THEN BEGIN sww=y_size/nscan ;swath width ENDIF ELSE BEGIN sww=yscl ENDELSE IF TOTAL(SHORTNAME EQ ["MOD021KM","MYD021KM","GSUB1","RMT021KM","RMA021KM"] AND $ extgeo EQ 0) THEN BEGIN g_off=2 ENDIF ELSE BEGIN g_off=0 ;geolocation offset ENDELSE y1=y1-(y1 MOD sww) y2=y2-(y2 MOD sww)+ sww-1 x1=x1-(x1 MOD xscl)+g_off x2=x2-(x2 MOD xscl)+g_off IF (y1 EQ y2) OR (x1 EQ x2) THEN BEGIN PRINTF,debug, "Identical start/end subset indices. Please, set start0)<(NXM-1) yv=(yv>0)<(NYM-1) ;Here we form the name of the BINARY ;and image file intended to accommodate global ;maps: IF(gallery EQ 0) THEN BEGIN IF(fn EQ 0) THEN BEGIN gname0=gname bfname=gname+"dat" jpgname=gname+"jpg" OPENW,2,outpath+bfname ENDIF ELSE OPENU,2,outpath+bfname ENDIF ; end if (gallery EQ 0 and fn EQ 0) IF(gallery and (True_col EQ 3)) THEN BEGIN bfname=gname+"dat" jpgname=gname+"jpg" OPENW,2,outpath+bfname ENDIF asar=ASSOC(2,FLTARR(NXM,NYM)) IF(autointen) THEN intens=FLTARR(2,totselch) FOR k=0,totselch-1 DO BEGIN ;=======SCALING EQUATION============= ; convert to physical units ;valid for L1B IF TOTAL(SHORTNAME EQ L1B) GT 0 THEN $ temp=radsc[k]*(image[*,*,k]-radof[k]) ;valid for atmosphere L2: IF TOTAL(SHORTNAME EQ AL2) GT 0 THEN $ temp=radsc[k]*(image[*,*,k]-radof[k]) ;valid for ocean L2 ;Oceans have power equation for ;cocco_conc_detach, ipar and par: IF TOTAL(SHORTNAME EQ OL2) GT 0 THEN BEGIN IF ( (pname[k] EQ "cocco_conc_detach") OR $ (pname[k] EQ "ipar") OR $ (pname[k] EQ "arp")) THEN BEGIN temp=10.^(radsc[k]*image[*,*,k]+radof[k]) ENDIF ELSE BEGIN temp=radsc[k]*image[*,*,k]+radof[k] ENDELSE ENDIF ;If cannot figure out the product, ;do not scale: IF TOTAL(SHORTNAME EQ [L1B,AL2,OL2]) EQ 0 THEN $ temp=image[*,*,k] ;==================================== ;--------autoset image brightness----------- IF(autointen) THEN BEGIN tmp=MEAN(temp) tmp1=STDEV(temp) mnm=MIN(temp) intens[*,k]=[mnm>tmp-tmp1,tmp+tmp1] ENDIF ;-------------------------------------------- ;==The case global map, no gallery== IF(gallery EQ 0) THEN BEGIN PRINTF,debug,"Map slice:",k ;mapping IF(fn EQ 0) THEN BEGIN map=REPLICATE(-999.,NXM,NYM) IF (intrpl) THEN MAPIT,map,xv,yv,temp,x_size,y_size,ia,ja $ ELSE map[xv,yv]=temp[*,*] asar[k]=map ENDIF ELSE BEGIN map=asar[k] IF (intrpl) THEN MAPIT,map,xv,yv,temp,x_size,y_size,ia,ja $ ELSE map[xv,yv]=temp[*,*] asar[k]=map ENDELSE ENDIF ;end if (gallery EQ 0) ;===================================== ;==The case no global map, yes gallery== IF(gallery) THEN BEGIN PRINTF,debug,"Map slice:",k IF(TRue_col EQ 1) THEN BEGIN bfname=gname+pname[k]+".dat" OPENW,2,outpath+bfname ENDIF ;end if (True_col eq 1) ;mapping map=REPLICATE(-999.,NXM,NYM) IF (intrpl) THEN MAPIT,map,xv,yv,temp,x_size,y_size,ia,ja $ ELSE map[xv,yv]=TEMPORARY(temp[*,*]) ENDIF ;end if gallery ;============================================ xtitle=pname[k] ;these are the titles for the gray unit=unita[k] ;scale. They may change for each ;k when doing gallery (e.g. ;L2 atmosphere) ;The following are bin and png file-writes ;intended for galleries. Don't move them ;from here. IF (gallery) THEN BEGIN IF(True_col eq 1) THEN BEGIN WRITEU,2,map CLOSE,2 rgb=map showimage,intens[*,k],intens[*,k],intens[*,k] IF(!VERSION.RELEASE GT 5.3) THEN BEGIN temp=outpath+gname+pname[k]+".png" WRITE_PNG,temp,zbuf,r,g,b ENDIF ELSE BEGIN temp=outpath+gname+pname[k]+".gif" WRITE_GIF,temp,zbuf,r,g,b ENDELSE PRINTF,debug,"Saved ",temp ENDIF ; end if true_col eq 1 IF(True_col EQ 3) THEN asar[k]=map ENDIF ;end if gallery ENDFOR ;--------end k-cycle (slices)---------------- CLOSE,2 ;close the output data binary file xv=0 yv=0 IF(gallery and True_col EQ 3) THEN BEGIN rgb=FLTARR(NXM,NYM,True_col) OPENR,2,outpath+bfname FOR k=0,True_col-1 DO rgb[*,*,k]=asar[k] CLOSE,2 SHOWIMAGE,intens[*,0],intens[*,1],intens[*,2] WRITEJPEG,outpath+jpgname PRINTF,debug,"Saved ",outpath+jpgname ENDIF ; end if gallery and true_col eq 3 ex: HDF_SD_END,sd_id IF (geosd_id NE sd_id) AND $ (geosd_id GT 0) THEN $ HDF_SD_END,geosd_id HDF_CLOSE,FileHandle ENDFOR ;********************end of fn cycle******************** ;***********THIS IS THE END OF HDF FILES CYCLING******** ;This reading from binary, and map show on the monitor ; will activate if True_col is not 0, and ;gallery=0, and debug=-1, i.e. when doing global map, ; gallery is turned off and debug mode is on. ; IF(True_col AND (gallery EQ 0) ) THEN BEGIN ;This will save RGB image ;of global map: IF (True_col EQ 3) THEN BEGIN rgb=FLTARR(NXM,NYM,True_col) OPENR,2,outpath+bfname FOR k=0,2 DO rgb[*,*,k]=asar[k] CLOSE,2 SHOWIMAGE,intens[*,0],intens[*,0>1],$ intens[*,0>2] WRITEJPEG,outpath+jpgname PRINTF,debug,"Saved ",outpath+jpgname ENDIF ;end if true_col eq 3 ;------------------------- ;This will save GIF images of ;all global maps (map per channel) IF (True_col EQ 1) THEN BEGIN rgb=FLTARR(NXM,NYM) OPENR,2,outpath+bfname FOR k=0,totselch-1 DO BEGIN rgb[*,*]=asar[k] xtitle=pname[k] unit=unita[k] SHOWIMAGE,intens[*,k],intens[*,0],$ intens[*,0] IF (!VERSION.RELEASE GT 5.3) THEN BEGIN temp=outpath+gname0+pname[k]+".png" WRITE_PNG,temp,zbuf,r,g,b ENDIF ELSE BEGIN temp=outpath+gname0+pname[k]+".gif" WRITE_GIF,temp,zbuf,r,g,b ENDELSE PRINTF,debug,"Saved ",temp ENDFOR CLOSE,2 ENDIF ;end if true_col eq 1 ;----------------------------- ENDIF IF(debug GT 0) THEN $ CLOSE,debug ;close log file CLOSE,5 ;close "filelist" SET_PLOT,dname RETURN END