How to programmatically rasterize and georeference point spectra in ENVI

 

If a user has point spectra data with associated map or geographic coordinates (e.g., ASD FieldSpec data collected with a GPS attached), it is possible to build the spectra into a georeferenced image within ENVI. This Help Article discusses how to go about accomplishing this programmatically and includes sample data and code to get you started.

Rasterizing and georeferencing point spectra is a multi-staged process. Below are the required steps:

1) For this example, the map and spectral data for each sample must be placed in a delimited ASCII text file, formatted in a very specific way. Each line in the file corresponds to one spectra. From left to right, the values should be latitude, longitude, band 1, band 2, band 3,.... Below is an example of how the first line of such file should look like, which uses tab delimited spacing between values:

38.96808337   -120.0809865   230   230   229   .....

2) Using the commented sample code provided below, a user can read the map and spectral data into IDL arrays, from which two images are created. In this example, the first contains spectral data for each sample in BIP format and has 41 bands. Each pixel contains an entire spectrum. The second, also in BIP format, has two bands (one for latitude, one for longitude) and has exactly the same X and Y dimensions as the first image. Each pixel contains the exact map/geographic location of the corresponding pixel in the first image. This second image will be used to create a Geographic Lookup Table (GLT) so that each spectrum can be properly placed in a georeferenced image.

3) Once the two images are created, the user can go to Map->Georeference from Input Geometry->Build GLT to start the georeferencing process. Select the Longitude band for X Geometry and the Latitude band for Y Geometry. In the next screen, provide the correct input projection and select an output projection. In this example, the input projection is Geographic Lat/Lon. In the next screen, select an output pixel size for the georeferenced image. The smaller the number, the larger the resulting file. Setting the rotation to 0 is usually best. Save the GLT file to disk.

4) Once the GLT file is saved, go to Map->Georeference from Input Geometry->Georeference from GLT. Select the new GLT file as the Input Geometry File and the spectral data file as the Input Data file. Select a background value, if desired, and provide an output filename for the georeferenced image.

Code Example:

;This program creates a  button in ENVI to access 
;the rasterization program

pro rasterize_spec_data_define_buttons, buttonInfo
    compile_opt idl2

	ENVI_DEFINE_MENU_BUTTON, buttonInfo, VALUE = 'Point Spectra', $
	   /MENU, REF_VALUE = 'File', /SIBLING, POSITION=10

	ENVI_DEFINE_MENU_BUTTON, buttonInfo, VALUE = 'Create Spectra and Geo Images', $
	   EVENT_PRO='rasterize_spec_data', POSITION=0, ref_value = 'Point Spectra', uvalue='Create Spectra and Geo Images'

end


;This function returns the dimensions of the 
;most "square" rectangular array possible based 
;on a single inputted value.  In this example, the 
;value corresponds to the number of spectra
;to be rasterized.  

function get_dimensions, x, fact=f
	compile_opt idl2
	
    if (n_elements(f) eq 0) then begin
        f = 1L
        flag = 1
    endif else $
        flag = 0
    t = floor(sqrt(x))
    p = 2L
    
    while (p le t) do begin
        if ((x mod p) eq 0) then begin
            f = get_dimensions(x/p, fact=[f,p])
            break
        endif else p = p+1L+(p gt 2)
    endwhile
    
    if (p gt t) then $
        f = [f,x]
    
    if flag then begin
        f=f[1:*]
        nFactors =n_elements(f)
        if nFactors le 2 then $
            return, f
        fCheck = sqrt(x)/2
        f1 = f[nFactors-1]
        index = 0
        while f1 le fCheck do begin
            f1 = f1*f[index]
            index++
        endwhile
        return, [f1,x/f1]
    endif
    
    return, f

end


pro rasterize_spec_data, event
	compile_opt idl2

    ;Provide an ASCII input file formatted in the following way:
    ;- There is no header information--just data
    ;- Each row contains all relevant data for one spectrum
    ;- Each row is set up so that the data, from left to right, are 
    ;   latitude, longitude, band 1, band 2, band 3....
    ;The sample ASCII file provided with this code (sampledata.txt) 
    ;contains 41 bands of data. 
    filename = dialog_pickfile(title='Select ASCII input file')
    if filename eq '' then return
	
    ;Strip full path from input filename, use input file's basename for output files
    base_name = file_basename(filename)
	
    ;Provide output path for processed data	
    output_path = dialog_pickfile(title='Select output directory', /directory)
    if output_path eq '' then return

    ;Read the formatted ASCII data into an IDL array.  Values will be in double precision
    ;floating point.  If there is a header, set the read_skip keyword.
    envi_read_cols, filename, data

    ;Retrieve dimensions of the data array, especially the number of spectra (lines).
    dimensions = size(data, /dimensions)

    ;Pass retrieved number of spectra to secondary function designed to 
    ;figure out the dimensions of the most "square" rectangular image 
    ;array in which to hold the spectra such that each pixel contains an
    ;entire spectrum (BIP).  This is necessary because ENVI's GLT routines
    ;require input data organized into more than one row.  
    extents = get_dimensions(dimensions[1])

    ;Check to see if the returned image array dimensions are greater
    ;than one row.  If not, that means the number of spectra is a prime
    ;number and the GLT process will not work.  Adding or removing 
    ;one spectrum in the ASCII file should do the trick.
    if n_elements(extents) eq 1 then begin
        ok = dialog_message('Please provide a number of samples that is not a prime number')
        return
    endif

    ;Create container arrays for the geographic information and spectral
    ;data for each sample in the ASCII file.
    geo_array = dblarr(dimensions[1], 1, 2)
    spectra_array = dblarr(dimensions[1], 1, 41)

    ;Fill container arrays with geographic and spectral data for each
    ;sample.
    for i = 0, dimensions[1]-1 do begin
	geo_array[i,0,*] = transpose(data[0:1,i])
	spectra_array[i,0,*] = transpose(data[2:*, i])
    endfor

    ;Reform both container arrays based on the optimal image array 
    ;dimensions returned from the secondary function.
    geo_array_reform = reform(geo_array, extents[1], extents[0], 2)
    spectra_array_reform = reform(spectra_array, extents[1], extents[0], 41)

    ;Write out ENVI files containing geographic and spectral data for each
    ;sample.  The geo image has two bands (one for latitude, one for longitude).
    ;The spectral image has 41 bands.  Both images have the same X and Y
    ;dimensions.  So, for example, the pixel located at (5,3) in both images contains
    ;information about the same spectrum.  This is necessary for ENVI's GLT
    ;process.  
    envi_write_envi_file, geo_array_reform, bnames=['Latitude (Y)','Longitude (X)'], $
    out_name=output_path+path_sep()+base_name+'_geo.img'
    envi_write_envi_file, spectra_array_reform, $
    out_name=output_path+path_sep()+base_name+'_spectra.img'

    ;Reminder to user that they must now use the newly created geo file to create
    ;a geographic lookup table (GLT) and apply it to the spectral data in order to
    ;georeference the spectra.
    ok=dialog_message('You must now build a GLT file and use it to georeference the spectra image')

end

___________________________________________
Reviewed by BC on 09/17/2014