13 Oct 2015 04:17 PM |
|
I am creating some vegetative indices with bandmath and then I am going to make histograms of them later on, but I have values of 0.00 representing no data that I want to exclude from analysis. Is there a way to tell the bandmath experssion to do this?
For example, I am using this code to create an VI,
pro NDII
envi, /restore_base_save_files
envi_batch_init, log_file='batch.txt'
;set up path to input and output files
input_path = 'D:\Sheyenne\Atmospherically Corrected Landsat\no_tree_individual_scenes\Stacked'
;change to output directory where I want the processed files to be placed
cd, 'D:\Sheyenne\Atmospherically Corrected Landsat\Indices\Main\no_trees\NDII'
; search for data files in the specified directory
files = FILE_SEARCH(input_path,'*.tif', count=count)
IF (count EQ 0) THEN BEGIN
Print, 'No files were found to process'
ENDIF
for i=0, count-1 do begin
envi_open_file, files[i], r_fid=fid
if (fid eq -1) then begin
envi_batch_exit
return
endif
envi_file_query, fid, dims=dims, ns=ns, nl=nl, fname=filename
;set pos array for calculating NDVI
pos = [2,1] - 1
;set the output file names by stripping out the base name and appending 'ndvi.img'
out_name = file_basename(filename, '.tif') + '_NDII.tif'
;call the doit
ENVI_Doit, 'Math_Doit', $
FID = [fid, fid], $
DIMS = dims, $
POS = [3,4], $
EXP = '(float(b3)-float(b4))/(float(b3)+float(b4))', $
OUT_NAME = out_name
endfor
end
but I want to change all values of band 3 (b3) and band 4 (b4) that are equal to 63355.0 to NaN before I execute the band math.
|
|
|
|
Zachary Norman Basic Member
Posts:173  
14 Oct 2015 12:29 PM |
|
Hi Stefano,
I'm not sure how you would go about this with the ENVI classic API, but in the new API I have an idea that may work for you. The majority of the code can be found here in a help article that I wrote:
http://www.exelisvis.com/Support/Help...
You can replace the code for the tile iterator with something like:
; Iterate through tiles
FOREACH tile, all_tiles DO BEGIN
;get all spectral data from original input_raster
all_bands = input_raster.GetData(SUB_RECT=all_tiles.CURRENT_SUBRECT)
;choose our bands
b1 = all_bands[*,0,0] ;band 1
b2 = all_bands[*,0,1] ;band 2
;do the band math
data = (float(b1) - b2)/(b1 + b2)
;say band1 has bad data values set to zero so we can get these index values and
;replace the band_math raster with another value (-999 for this example). the
;number -999 then corresponds to DATA_IGNORE_VALUE in ENVI so we won't see it.
bad_index_locations = where(b1 eq -999)
data[bad_index_locations] = -999
;save the tile to the output raster
output_raster.SetTile, data, all_tiles
ENDFOREACH
This will perform your band math on an image tile and set the result of the band math to the original bad data value. Because band math is just adding arrays, the bad data value will probably change depending upon your bad data value. To get around this, we can get the data values in the first band that are bad and set the result from the band math to the same bad data value. To extend this, when you create an ENVIRaster object with the new ENVI API, you can specify the property DATA_IGNORE_VALUE which makes ENVI not display those pixels.
However, the pixels will still be used in other operations unless you create a mask raster to ignore processing on them.
|
|
|
|
Deleted User New Member
Posts:  
19 Oct 2015 02:25 PM |
|
This seems like it is what I want, I am trying to modify my code to use it but I am getting an error on the line:
all_tiles = (ENVIsubsetRaster(input_path,$
BANDS= [0])).CreateTileIterator(MODE='spectral')
which says Input raster is missing or invalid even though the folder it is being directed to has all the .tif files I want in it.
The exact code I am using is:
pro NDII2
envi, /restore_base_save_files
envi_batch_init, log_file='batch.txt'
;set up path to input and output files
input_path = 'H:\Sheyenne\Atmospherically Corrected Landsat\no_tree_individual_scenes\stacked_copy'
;change to output directory where I want the processed files to be placed
cd, 'H:\Sheyenne\Atmospherically Corrected Landsat\Indices\Main\no_trees\NDII'
; search for data files in the specified directory
files = FILE_SEARCH(input_path,'*.tif', count=count)
IF (count EQ 0) THEN BEGIN
Print, 'No files were found to process'
ENDIF
for i=0, count-1 do begin
envi_open_file, files[i], r_fid=fid
if (fid eq -1) then begin
envi_batch_exit
return
endif
envi_file_query, fid, dims=dims, ns=ns, nl=nl, fname=filename
;set pos array for calculating NDVI
pos = [2,1] - 1
;set the output file names by stripping out the base name and appending 'ndvi.img'
out_name = file_basename(filename, '.tif') + '_NDII.tif'
all_tiles = (ENVIsubsetRaster(input_path,$
BANDS= [0])).CreateTileIterator(MODE='spectral')
; Iterate through tiles
FOREACH tile, all_tiles DO BEGIN
;get all spectral data from original input_raster
all_bands = input_path.GetData(SUB_RECT=all_tiles.CURRENT_SUBRECT)
;choose our bands
b1 = all_bands[*,0,0] ;band 1
b2 = all_bands[*,0,1] ;band 2
;do the band math
data = (float(b1) - b2)/(b1 + b2)
;say band1 has bad data values set to zero so we can get these index values and
;replace the band_math raster with another value (-999 for this example). the
;number -999 then corresponds to DATA_IGNORE_VALUE in ENVI so we won't see it.
bad_index_locations = where(b1 eq 65535.0)
data[bad_index_locations] = -999
OUT_NAME = out_name
ENDFOREACH
endfor
end
|
|
|
|
Zachary Norman Basic Member
Posts:173  
20 Oct 2015 08:53 AM |
|
I bet your problem is coming from using the old ENVI API instead of the new API for ENVI. To show you how this is done, I added a foreach loop around my example from the help article which should do what you want. Here is the code:
foreach file, files do begin
;input raster that we will perform bandmath on
infile = file
infile_path = FILEPATH(infile, ROOT_DIR=e.ROOT_DIR, $
SUBDIRECTORY = ['data'])
;output raster filename that we will save to
;NDI = Normalized Difference Index
outfile = 'NDI_' + infile
outfile_path = thisdir + path_sep() + outfile
;we cannot write to the output file if it exists, so we need to delete it
if (file_which(thisdir, outfile) ne '') then FILE_DELETE, outfile_path
;Load the input raster
input_raster = e.OpenRaster(infile_path)
;Set up output raster to be the same as the input input_raster
;except for number of bands. Specifically, we need
;the spatial reference to be the same
output_raster = ENVIRaster(URI = outfile_path, $
NROWS=input_raster.NROWS, $
NCOLUMNS=input_raster.NCOLUMNS, $
NBANDS=1, $
SPATIALREF=input_raster.spatialref, $
DATA_TYPE=4, DATA_IGNORE_VALUE = -999)
;make a tile iterator of first band for input raster because our output
;raster will only have one band
all_tiles = (ENVIsubsetRaster(input_raster,$
BANDS= [0])).CreateTileIterator(MODE='spectral')
; Iterate through tiles
FOREACH tile, all_tiles DO BEGIN
;get all spectral data from original input_raster
all_bands = input_raster.GetData(SUB_RECT=all_tiles.CURRENT_SUBRECT)
;choose our bands
b1 = all_bands[*,0,0] ;band 1
b2 = all_bands[*,0,1] ;band 2
;do the band math
data = (float(b1) - b2)/(b1 + b2)
;say band1 has bad data values set to zero so we can get these index values and
;replace the band_math raster with another value (-999 for this example). the
;number -999 then corresponds to DATA_IGNORE_VALUE in ENVI so we won't see it.
bad_index_locations = where(b1 eq -999)
data[bad_index_locations] = -999
;save the tile to the output raster
output_raster.SetTile, data, all_tiles
ENDFOREACH
;save the output_raster
output_raster.save
endforeach
Another thing I though I would mention is your use of file search. It may be more useful to do something like this when searching for files:
cd, current = first_dir
cd, input_path
files = FILE_SEARCH('*.tif', count=count)
cd, first_dir
The reason I mention this is because the way you have it set up file_search will recursively search input_path. I forget this sometimes and it causes some problems in my code because I search other directories that are in that recursive search path. This then doubles some of the work that the program has to do because it is doing steps multiple times.
-Zach (VIS)
|
|
|
|
Deleted User New Member
Posts:  
21 Oct 2015 02:03 PM |
|
That works great, thank you! Its hard to find places to learn Envi IDL so you posting code and me working through what is going on was really helpful!
EDIT:
Even though Envi now recognizes they pixels as nodata, the value of -999 is still included when I look at the statistics. Looks like I will have to figure out a way to create masks for each of the 80 bands in my stacked images and apply the masks to the correct bands before calculating statistics
|
|
|
|
Zachary Norman Basic Member
Posts:173  
22 Oct 2015 09:57 AM |
|
Hi Stefano,
Glad to hear you got it to work! The new ENVI API is a little confusing when you first start learning about it. How is it that you are calculating your statistics for your images? I may be able to help give some ideas to filter those values out.
|
|
|
|
Deleted User New Member
Posts:  
22 Oct 2015 11:43 AM |
|
What I have been doing is to calculate statistics based on roi (I have 54 roi's), and I calculate them on a stacked image of a VI. I haven't been using IDL at all to do this so far but looks like I will need to. So for instance I have 80 landsat images in the stack with each band representing NDVI from different datesand the stats per ROI will give me a histogram for each ROI and all 80 bands as a text file. The problem with this is that the mask for each of the 80 images is different (they are clouds). Another option is to exclude the -999 from the histograms, which is what this article seems to do, http://www.exelisvis.com/Company/Pres... so I may look into that as well.
|
|
|
|
Deleted User New Member
Posts:  
22 Oct 2015 03:01 PM |
|
I have been playing around with the ENVIRasterStatistics function and it seems that values of -999 aren't incorporated in this. I am using this code:
e=ENVI()
File=Filepath('NDII_LT50300281984137PAC00.dat', Root_dir='H:\Sheyenne\Atmospherically Corrected Landsat\Indices\Main\no_trees\NDII')
Raster=e.OpenRaster(File)
stats = ENVIRasterStatistics(input_raster, /HISTOGRAMS, HISTOGRAM_BINSIZE=.2)
print, stats
which returns:
MIN: -0.44141799211502075
STDDEV: 0.12945765008883700
CORRELATION: !NULL
HISTOGRAMS: LIST
MAX: 0.51628601551055908
EIGENVALUES: !NULL
MEAN: -0.099569855931613105
NPIXELS: 868868
EIGENVECTORS: !NULL
COVARIANCE: !NULL
Since the minimum value is not -999 it must be registering it as no data. From what I have read the Histogram in this is within a dictionary as a list, would you know how to plot the histogram stored within by chance? Or better yet, If i could save the list as a csv file from the dictionary I can do what I want in python which I am more comfortable with.
|
|
|
|
Zachary Norman Basic Member
Posts:173  
22 Oct 2015 03:50 PM |
|
Hi Stefano,
Here is a good example of how to get the actual histogram information. It will generate plots of the histograms for each band of the input image. The variable bin_count contains the information for the frequency in each bin. The bin locations are represented by xbins, which need to be created by hand because they don't get created by default.
compile_opt idl2
ireset, /no_prompt
e=envi(/headless)
; Launch the application
e = ENVI(/headless)
; Create ENVIRaster statistics
file = FILEPATH('qb_boulder_msi', ROOT_DIR=e.ROOT_DIR, $
SUBDIRECTORY = ['data'])
raster = e.OpenRaster(file)
; Return the statistics
stats = ENVIRasterStatistics(raster, HISTOGRAM_BINSIZE = 2)
;histograms
histograms = stats.histograms
for i=0, n_elements(histograms)-1 do begin
band_hist = histograms[i]
xbins = band_hist.MIN + $
findgen(band_hist.BINCOUNT, INCREMENT = band_hist.BINSIZE)
bin_count = band_hist.COUNTS
hist = plot(xbins, bin_count, /test, $
WINDOW_TITLE = 'Band ' + strtrim(band_hist.BAND,2) + ' Histogram')
endfor
I would be a little cautious with the bin size for the histograms and/or specifying the HISTOGRAM keyword. I was having some crashes when I used both like you were, but I think it may have been because my bin size was way too small.
|
|
|
|
Deleted User New Member
Posts:  
22 Oct 2015 06:12 PM |
|
Perfect, thanks again! I think I will just subset the data based on the location of the corner pixels for each roi to create the stats per roi and I should be good to go then.
|
|
|
|