pro SL_alignment_v5 ; SL_alignment_v5, 21 July 2020 ; This version written by Daniel Politte ; (Based on April 2017 version by Valerie Fox, SL_alignment_v2.pro) ; ; Reads in an L cube and an S cube, and then combines them into one. ; The L cube does NOT need to be reversed (i.e., wavelengths made to be ; ascending) before this script is run. ; ; Unlike previous versions, this one does not second-guess the map projection ; info in the inputs' ENVI headers. This means that, unless coregistration has ; already been done on the inputs, the S and L portions will not line up right. ; ; The output only includes data in spatial locations where both S and L data are ; available. If either spectrometer's data is missing, we fill that pixel with the ; no-data value. ; Acquire filenames from the user infileL = dialog_pickfile(title="Select L data", /MUST_EXIST) if ((strlen(infileL) eq 0)) then begin print, 'Operation Canceled... Returning' return endif infileS = dialog_pickfile(title="Select S data", /MUST_EXIST) if ((strlen(infileS) eq 0)) then begin print, 'Operation Canceled... Returning' return endif outfile = dialog_pickfile(/WRITE, /OVERWRITE_PROMPT, DEFAULT_EXTENSION="img") if ((strlen(outfile) eq 0)) then begin print, 'Operation Canceled... Returning' return endif ; End acquiring filenames e = ENVI() ;Read in files and determine dimensions print, 'initializing' Ldata_in= e.OpenRaster(infileL) if Ldata_in.metadata.hasTag('wavelength') then begin ; We want the L data's to be arranged with ascending wavelengths waves = Ldata_in.metadata['wavelength'] if waves[0] > waves[1] then begin ; It's indeed descending. Reverse it! reverse_raster_bands, Ldata_in, Ldata endif else begin ; It's already ascending. Leave it be. Ldata = Ldata_in endelse endif else begin ; This must be parameter maps or something. We acutally don't want to reverse the bands, then. Ldata = Ldata_in endelse L_ns = Ldata.ncolumns L_nl = Ldata.nrows L_nb = Ldata.nbands L_spatref = Ldata.spatialref Sdata = e.OpenRaster(infileS) S_ns = Sdata.ncolumns S_nl = Sdata.nrows S_nb = Sdata.nbands s_spatref = Sdata.spatialref print, 'L columns, rows, bands: ', [L_ns, L_nl, L_nb] print, 'S columns, rows, bands: ', [S_ns, S_nl, S_nb] print, 'stacking' ;Define the spatial reference for the stacked file. ;S data is forced to conform to L data. ; ;Example Coordinate System info ; ENVISTANDARDRASTERSPATIALREF <339023> ;COORD_SYS_CODE = 0 ;COORD_SYS_STR = 'PROJCS["Mars Equirectangular Default",GEOGCS["GCS_Unknown",DATUM["D_Unknown",SPHEROID["S_Unknown",3396190.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Equidistant_Cylindrical"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],UNIT["Meter",1.0]]' ;PIXEL_SIZE = 12.000000, 12.000000 ;ROTATION = 0.00000000 ;TIE_POINT_MAP = 8138379.1, -277881.29 ;TIE_POINT_PIXEL = 0.50000000, 0.50000000 ; coordsys_code = L_spatref.coord_sys_code coordsys_string = L_spatref.coord_sys_str Coordsys = ENVICoordSys(coord_sys_str=coordsys_string) print, 'Coordinate System: ', coordsys pixelsize = L_spatref.Pixel_Size tiePointMap = L_spatref.Tie_Point_Map tiePointPixel = L_spatref.Tie_Point_Pixel ;Now define the grid: Output_Grid = ENVIGridDefinition(Coordsys, PIXEL_SIZE=pixelsize, TIE_POINT_MAP=tiePointMap, TIE_POINT_PIXEL=tiePointPixel, NROWS=L_nl, NCOLUMNS=L_ns) ;Create the grid for both input rasters L_layers_spatialgridraster = ENVISpatialGridRaster(Ldata, grid_definition=output_grid) S_layers_spatialgridraster = ENVISpatialGridRaster(Sdata, grid_definition=output_grid) ;Now finally stack them LayerStack = ENVIMetaspectralRaster([S_layers_spatialgridraster,L_layers_spatialgridraster]) ; Extract the data cube = LayerStack.getData() fillvalue = LayerStack.metadata.hasTag('data ignore value') ? LayerStack.metadata['data ignore value'] : 0.0 cubesize = size(cube) ; Blank out areas where EITHER of our two input cubes have no data foo = L_layers_spatialgridraster.getdata(pixel_state=l_nodata,bands=[0]) foo = S_layers_spatialgridraster.getdata(pixel_state=s_nodata,bands=[0]) for r=0,cubesize[1]-1 do begin ;loop over rows for c=0,cubesize[2]-1 do begin ;loop over columns if l_nodata[r,c] or s_nodata[r,c] then begin cube[r,c,*] = fillvalue endif endfor endfor ; create a raster object with the modified data & save it to file LayerStack_cropped = ENVIRaster(cube, uri=outfile, metadata=LayerStack.metadata, spatialref=LayerStack.spatialref) LayerStack_cropped.save print, 'done!' end pro reverse_raster_bands, raster_in, raster_out compile_opt idl2 e = ENVI() inmeta = raster_in.metadata indata = raster_in.getData() outmeta = ENVIRasterMetadata.Hydrate(inmeta.Dehydrate()) ; Reverse the standard bandwise properties of the raster. ; WARNING: will NOT work with any other band-wise properties there might be in the header. if outmeta.hasTag('band names') then begin outmeta['band names'] = reverse(outmeta['band names']) endif if outmeta.hasTag('wavelength') then begin outmeta['wavelength'] = reverse(outmeta['wavelength']) endif if outmeta.hasTag('fwhm') then begin outmeta['fwhm'] = reverse(outmeta['fwhm']) endif if outmeta.hasTag('bbl') then begin outmeta['bbl'] = reverse(outmeta['bbl']) endif ; Just kill the default bands if outmeta.hasTag('default bands') then begin outmeta.removeItem, 'default bands' endif ; Reverse the raster's data, and put the new raster together outdata = reverse(indata, 3) newFile = e.getTemporaryFilename() ; Just a temp file raster_out = ENVIRaster(outdata, metadata=outmeta, uri=newFile, spatialref=raster_in.spatialref) raster_out.save end