X
PrevPrev Go to previous topic
NextNext Go to next topic
Last Post 13 Oct 2015 04:17 PM by  anon
Possible to ignore certain values in bandmath?
 9 Replies
Sort:
You are not authorized to post a reply.
Author Messages

anon



New Member


Posts:
New Member


--
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
    Basic Member


    --
    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:
    New Member


    --
    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
    Basic Member


    --
    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:
    New Member


    --
    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
    Basic Member


    --
    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:
    New Member


    --
    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:
    New Member


    --
    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
    Basic Member


    --
    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:
    New Member


    --
    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.
    You are not authorized to post a reply.