Don’t forget to Stretch! Using ENVI’s stretch tools to see things our eyes can’t.

Author: Jeff McKissick

Living in Boulder, we mountain people out here like to do a lot of physical activities whether it’s hiking, skiing, or yoga. Everyone knows the first thing you have to do before any physical activity is STRETCH! This also applies in ENVI as well! Over the past few months I have worked on various projects where, had I applied one of our stretches in ENVI first, I would have saved a lot of time for myself. This example today was a dataset of a large grass field in which the user was looking for an invasive species weed within this field.

You can see from the figure above that EVERYTHING LOOKS GREEN! How can you pick out a weed when everything looks like grass? With a little help from the customer, we were able to get access to a shapefile they provided that showed us areas in the scene that actually were the weed we were looking for. Still, even with these shapefiles everything looks the same color. This is where, before you start any of your preprocessing or classification workflows, you stretch!

ENVI has some really great stretch tools to choose from, but seeing them isn’t actually helping you know what they mean. For this example we used a few different linear percent stretches to help accentuate some of our features. What these percent stretches do is trim the X% of extreme values at the beginning and end of the histogram.

So for example, if you look at our three images with the histogram stretch plot shown, you can see in the first image with no stretch that our pixel values are 0-255 which is standard. If we look at our Linear 2% and 5% stretched images respectively you see the pixel values get trimmed on each end of each color band.

From here we were easily able to identify the invasive weed in our scene and compare it to the shapefiles provided for us so that we could run a classification workflow and extract the features that we wanted. Our shapefiles, not shown here, were all around the areas in the scene above that were a very dark green. These stretches allowed us to make more accurate ROIs  (Regions of Interest) for our classification which in turn gave us a more accurate result.


Comments (0) Number of views (1461) Article rating: 5.0

Categories: ENVI Blog | Imagery Speaks





Interactive Vector GUI

Author: Austin Coates

With holidays rapidly approaching, I thought it might be nice to present a light hearted demonstration of capabilities; something that you can delight friends and family with during cold winter nights.  This week’s post shows how to use a simple widget_window to create an interactive vector animation.  The code is broken up into four pieces VectorBlog, VectorBlog_Event, VectorBlogMouseEvent, and VectorBlogVectorBuilder. The VectorBlog creates the widget itself and initially displays the vector information.  VectorBlog_Event runs the events for the base widget.  This is mainly the closing of the widget. VectorBlogMouseEvent is the mechanism by which the mouse locations are recoded by the widget window. Last but not least VectorBlogVectorBuilder builds the two matrices containing the horizontal and vertical components of the vector.  When the following code is run the cursor location is recorded, a distance weighted vector diagram is created, and the resulting image is then displayed in the widget window.  The only tricky bits in this code relate to the fact that there are two event handlers, one for the base widget and one for the mouse events in the widget window.  The object reference for the image displayed must also be maintained and manipulated each time the cursor is move.


Figure 1: Example Vector Plot


; Build the vecort components
; This is distance weighting
pro VectorBlogVectorBuilder, dims = dims, pt = pt, u = u, v = v
  compile_opt IDL2

  ; Get all the x locations
  xs = findgen(dims[0])
  xs = xs#replicate(1,dims[1])

  ; Get all the y locations
  ys = findgen(dims[1])
  ys = transpose(ys# replicate(1,dims[0]))

  ; Build the horizontal components
  u = pt[0] - xs

  ; Build the vertical components
  v =  pt[1] - ys

; The event handler called when the mouse is moved.
FUNCTION VectorBlogMouseEvent, win, x, y, keymods

  ; Set the window

  ; Set the scale factor
  sc = .05

  ; Get the size of the window
  dims = win.DIMENSIONS * sc

  ; Get the vecotr information
  VectorBlogVectorBuilder, dims = dims, pt = [x*sc,y*sc], $
    u = u , v = v

  ; Build the vecotr image
  vector = vector(u,v,xrange=[0,dims[0]], yrange=[0,dims[1]], /current, AUTO_COLOR=1,  $

  ; Delete the old image

  ; Display the new image
  win.uvalue = vector


; The event handler for the widget
pro VectorBlog_Event, ev
  compile_opt IDL2


pro VectorBlog
  compile_opt IDL2

  ; Setup the widget
  tlb = widget_base()
  wwindow = widget_window(tlb, MOUSE_MOTION_HANDLER='VectorBlogMouseEvent', $
    XSIZE=500, YSIZE=500)
  widget_control, tlb, /REALIZE

  ; Set the scale factor
  sc = .05

  ; Get infomation about the widget window
  dims = w.DIMENSIONS * sc

  ; Build the vecotr infomation
  VectorBlogVectorBuilder, dims = dims, pt = [(dims[0])/2.,(dims[1])/2.], $
    u = u , v = v

  ; Display the vectors
  vector = vector(u,v,xrange=[0,dims[0]], yrange=[0,dims[1]], /current, AUTO_COLOR=1,  $

  ; Record the vector object
  w.uvalue = vector

  ; Start the task manager
  XMANAGER, 'VectorBlog', tlb


Comments (1) Number of views (1031) Article rating: No rating

Categories: IDL Blog | IDL Data Point





Maintaining Backward Compatibility in IDL 8.6 - Keep Calm And Read Your Release Notes

Author: Jim Pendleton

One of the "benefits" of being a Harris Geospatial Solutions insider is access to pre-release candidates of the commercial products developed by our engineering group.

We in Custom Software Solutions are sometimes the canaries in the coal mine, learning on occasion that an "undocumented feature" that we had used benignly or even to our advantage has been removed from the language in a newer release. Generally, these changes are justifiable. 

Where to Find Release Notes

Sometimes it's not enough to look at the documentation center's What's New help page to learn of all these changes. An additional source of information, generally prepared after the "What's New" documentation has gone to press, is located in a place other than your IDL or ENVI installation's documentation subdirectory.

This file includes information about supported platforms and potential backward-compatibility issues.

If you have received a Harris product installation DVD, check the info subdirectory on the DVD itself for the release notes files.

If your installation was downloaded from the Harris Geospatial Solutions Download and License Center, you or your site's designated license administrator will need to retrieve the release notes from a link that is separate from the product installer. 

Downloading Release Notes

After logging in, select the "Harris" link near the bottom of the web page, under "Browse My Software and Documentation".

Select the "IDL" link, that pops up in a new list in the "Product Lines" column.

Select the "IDL" link in the "Current Releases" tab.

On the "Product Download" page, select the appropriate item for the Release Notes document.

An Example Backward Compatibility Note

Recently, I discovered some of my routines were using an admittedly illegal syntax involving the "_REF_EXTRA" keyword passing mechanism.  The code wasn't so much illegal, as it was ignored. And one could argue the compiler should have complained about it from the time that _REF_EXTRA was added to the language. For example,

function MyRoutine, _REF_EXTRA=extra

    MySubroutine, _REF_EXTRA=extra


Can you spot the problem?

The _REF_EXTRA keyword is only intended to appear in the declaration of a function or procedure.Within the body of the code, you should always use the _EXTRA keyword when making calls to other routines.

In the form shown above, the code is basically ignored at execution time.  It serves no useful purpose.

Up until IDL 8.6, illegal use of the _REF_EXTRA syntax would simply be skipped by the compiler and interpreter. In the example above, MySubroutine would be called without any keywords, regardless of what was passed to MyRoutine.

In IDL 8.6, the compiler has been beefed up to complain about the invalid syntax. If you have code that fits this pattern your code will not compile. You may want to simply remove the flagged code because it has never been operational. Or you may want to change the syntax.

function MyRoutine, _REF_EXTRA=extra

    MySubroutine, _EXTRA=extra


Carefully consider the implications of changing the keyword, however. Modifying the syntax will also alter the behavior. You may end up modifying keywords on output that you hadn't intended to change!

Comments (0) Number of views (1221) Article rating: No rating

Categories: IDL Blog | IDL Data Point





Optimizing Max Kernel Operation in IDL

Author: Atle Borsholm

I found this optimization question on the comp.lang.idl-pvwave, and decided to give it a try. The question was to implement an algorithm that replaces every element in a 2-D array with its neighborhood maximum value. In this case the neighborhood size was 101x101 (i.e. -50 to +50), and the array size was 3200x3248. The nested FOR loop approach looks like the following code snippet.

  data = randomu(seed,3200,3248)
  dim = size(data,/dimension)
  nx = dim[0]
  ny = dim[1]
  t0 = tic('Nested FOR')
  result2 = data
  for i=0,nx-1 do begin
    for j=0,ny-1 do begin
      result2[i,j] = max(data[(i-50)>0:(i+50)<(nx-1),(j-50)>0:(j+50)<(ny-1)])

My first thought was to use the > operator which returns the maximum of 2 arguments. It operates on arrays, and in conjunction with the SHIFT function it serves to return the larger of 2 neighbors. The other trick here is that since we are looking for a 101x101 neighborhood maximum, we can use a combination of smaller neighborhood maxima as input in a structured way in order to achieve the exact 101x101 neighborhood size. The code that I ended up with after some trial and error was the following.

  t0 = tic('Iterative >')
  ; Using SHIFT and > in an iterative way
  padded = replicate(min(data),size(data,/dimension)+100)
  padded[50,50] = data
  tmp3 = shift(padded,1,0) > padded > shift(padded,-1,0)
  tmp9 = shift(tmp3,3,0) > tmp3 > shift(tmp3,-3,0)
  tmp27 = shift(tmp9,9,0) > tmp9 > shift(tmp9,-9,0)
  tmp81 = shift(tmp27,27,0) > tmp27 > shift(temporary(tmp27),-27,0)
  tmp99 = shift(tmp9,44,0) > temporary(tmp81) > shift(temporary(tmp9),-44,0)
  tmp101 = shift(tmp3,49,0) > temporary(tmp99) > shift(temporary(tmp3),-49,0)
  ; Same for Y-dim
  tmp3 = shift(tmp101,0,1) > tmp101 > shift(temporary(tmp101),0,-1)
  tmp9 = shift(tmp3,0,3) > tmp3 > shift(tmp3,0,-3)
  tmp27 = shift(tmp9,0,9) > tmp9 > shift(tmp9,0,-9)
  tmp81 = shift(tmp27,0,27) > tmp27 > shift(temporary(tmp27),0,-27)
  tmp99 = shift(tmp9,0,44) > temporary(tmp81) > shift(temporary(tmp9),0,-44)
  tmp101 = shift(tmp3,0,49) > temporary(tmp99) > shift(temporary(tmp3),0,-49)
  result1 = (temporary(tmp101))[50:50+nx-1,50:50+ny-1]

I didn’t say that optimized code always looks pretty, but the goal here is to run fast. Adding in some result comparison checking to make sure the results are equivalent.

  print, array_equal(result1,result2) ? 'Results are matching' : 'SOMETHING went wrong'

Finally, here are the results, which yielded an impressive amount of speed-up, the execution time went from 189.4 seconds to 1.4 seconds, and the results are identical:

  % Time elapsed Nested FOR: 189.39470 seconds.
  % Time elapsed Iterative >: 1.4241931 seconds.
  Results are matching

Comments (0) Number of views (1112) Article rating: No rating

Categories: IDL Blog | IDL Data Point





New in ENVI 5.4: ENVITask Returning a Virtual Raster

Author: Brian Griglak

Since the creation of custom ENVITasks in ENVI 5.2 SP1, there has been the requirement that your procedure must commit all output objects to disk.  There was the rule that the procedure wrapped by the task must have an input keyword telling it the filename to use to write the object to.  In the task template, you would have an output parameter that mapped to that keyword, and then during task loading the framework would magically create an input parameter for you, mapped to the same keyword with TYPE set to “ENVIURI”. 

The automatic creation of the input parameter and the internal correlation of the two parameters were done with the best of intentions, to simplify the process of creating custom tasks.  Alas the user feedback on this feature wasn’t always as rosy as we hoped.

So in ENVI 5.4 we’re shaking things up and giving you the task developer more control.  If you still use the “version” property in the task template, and have it set to “5.2.1”, “5.3”, or “5.3.1”, then you’ll get the old behavior.  But if you switch to using “schema” set to the value “envitask_3.0”, then a new set of rules apply to the procedure, and what you can do inside it.  In the new paradigm, your procedure will have separate keywords for the input filename and output object reference.  If you like, you can skip the filename keyword completely and return an object that hasn’t been tethered to disk at all.  This makes life much easier for types like ENVIGCPSet and ENVITiePointSet, but also allows for a procedure that constructs a complex virtual raster chain based on the other input parameters.

You might be asking why you would want a task to return a virtual raster.  With the new Classification Framework that is part of ENVI 5.4, you need to make sure that you prepare the data you want to run through a trained classifier the same way you prepared your training data.  One way to do this is to create a task that returned the preprocessed data in a consistent way.  If you can get away with not having to save that preprocessed attribute raster to disk, why not take advantage of the time and space advantages of using a virtual raster.

The following example is a modified version of the new Custom Task Tutorial in the ENVI 5.4 release.  That task wraps a procedure that uses a number of other tasks internally to perform all the preprocessing.  Here I’ve modified it to use virtual rasters and other ENVI API function calls to avoid ever writing a file to disk.

The code goes through a number of steps, which I will describe after showing it:

pro SentinelVegIndices_blog, INPUT_RASTER_10M=raster10m, $
                             INPUT_RASTER_20M=raster20m, $
  compile_opt idl2, hidden
  ; Get the spatial reference of the 10-meter raster
  spatialRef = raster10m.SPATIALREF
  coordSys = ENVICoordSys(COORD_SYS_CODE=spatialRef.COORD_SYS_CODE)
  ; Create a spatial grid definition
  grid = ENVIGridDefinition(coordSys, $
                            PIXEL_SIZE=spatialRef.PIXEL_SIZE, $
                            NCOLUMNS = raster10m.NCOLUMNS, $
                            NROWS = raster10m.NROWS, $
                            TIE_POINT_MAP=spatialRef.TIE_POINT_MAP, $
                            TIE_POINT_PIXEL = spatialRef.TIE_POINT_PIXEL)
  ; Regrid the 20-meter bands to 10 meters
  regriddedRaster = ENVISpatialGridRaster(raster20m, GRID_DEFINITION=grid)
  ; Create a layer stack
  layerStack = ENVIMetaspectralRaster([raster10m, regriddedRaster], $
  ; Compute image statistics
  stats = ENVIRasterStatistics(layerStack)
  ; Perform dark subtraction as an alternative to atmospheric correction
  bandRasters = ObjArr(layerStack.nBands)
  for i = 1, layerStack.nBands do begin
    expression = 'b' + i.ToString() + ' - ' + stats.Min[i-1].ToString()
    bandRasters[i-1] = ENVIPixelwiseBandMathRaster(layerStack, expression)
  bandStackRaster = ENVIMetaspectralRaster(bandRasters, $
  ; we need to put the wavelengths back into the band stack,
  ; they were removed by the band math
  metadata = ENVIRasterMetadata()
  metadata['wavelength'] = layerStack.Metadata['wavelength']
  metadata['wavelength units'] = layerStack.Metadata['wavelength units']
  correctedRaster = ENVIMetadataOverrideRaster(bandStackRaster, $
  ; Scale pixel values from 0 to 1
  gains = MAKE_ARRAY(layerStack.NBANDS, /FLOAT, VALUE=0.0001)
  offsets = FLTARR(layerStack.NBANDS)
  scaledRaster = ENVIGainOffsetRaster(correctedRaster, gains, offsets)
  ; Create vegetation indices
  indices = [ 'Enhanced Vegetation Index', $
              'Global Environmental Monitoring Index', $
              'Leaf Area Index', $
              'Plant Senescence Reflectance Index', $
              'Red Edge Normalized Difference Vegetation Index' ]
  outputRaster = ENVISpectralIndexRaster(scaledRaster, indices)

The first step is to upsample the 20m raster to 10m pixels, which in the tutorial is performed using the ENVIRegridRasterTask.  This can be done with the ENVISpatialGridRaster virtual raster, once we have constructed an ENVIGridDefinition with the appropriate mixing of properties form the 10m and 20m rasters.

Next the tutorial uses the ENVIBuildBandStackTask to build a metaspectral stack of all the 10m bands.  Here we use the ENVIMetaspectralRaster virtual raster, though we have to pass in the original 10m raster’s spatial reference to keep this raster registered on the map.

Next is the dark subtraction.  The tutorial uses the ENVIRasterStatisticsTask, but here we just use the API function ENVIRasterStatistics() to accomplish the same thing.

The band minima values are used to perform dark object subtraction.  The tutorial uses the ENVIDarkSubtractionCorrectionTask, which handles this as a single raster.  Here I had to build separate band math equations for each band, as the ENVIPixelwiseBandMathRaster virtual raster always returns a single band output, so I have to build an array of band math expressions and then metaspectrally stack the results, again passing in the spatial reference.

One pitfall of the band math is that it removes most of the spectral metadata, so I have to put the wavelength metadata back into the raster so the spectral index calculations select the correct bands.  I do this with ENVIMetadataOverrideRaster(), using a copy of the original metaspectral raster’s metadata values.

The raster is then scaled down by a factor of 10000.0 with ENVIGainOffsetRaster, to simulate the atmospheric correction better and produce spectral index values that are more accurate.  Lastly, the scaled raster is passed into ENVISpectralIndexRaster to calculate 5 different spectral indices.

Once this has been written, we can use an almost identical .task file to wrap the procedure.  The main difference is that we no longer need to specify the OUTPUT_URI parameter, and I had a slightly different procedure name.  This task’s output raster is never commited to disk, but it can be used as input to another task by calling ENVIRaster::Dehydrate() on it, which yields a 25KB JSON representation.

Here is the updated .task file:

  "name": "SentinelVegIndices_blog",
  "base_class": "ENVITaskFromProcedure",
  "routine": "SentinelVegIndices_blog",
  "display_name": "Compute Sentinel-2A Vegetation Indices",
  "description": "This task regrids the visible and near-infrared bands of a Sentinel-2A Level-1C 20-meter image to 10 meters. It creates a layer stack of all bands. It applies dark-object subtraction as a simple alternative to atmospheric correction. It scales the reflectance pixel values from 0 to 1, then computes a select group of vegetation indices.",
  "schema": "envitask_3.0",
  "parameters": [
    "name": "INPUT_RASTER_10M",
    "display_name": "Select a 10-meter image",
    "type": "ENVIRASTER",
    "direction": "input",
    "required": true,
    "description": "Select a Sentinel Level-1C 10-meter image."
    "name": "INPUT_RASTER_20M",
    "display_name": "Select a 20-meter image",
    "type": "ENVIRASTER",
    "direction": "input",
    "required": true,
    "description": "Select a Sentinel Level-1C 20-meter image."
    "name": "OUTPUT_RASTER",
    "display_name": "Output image",
    "type": "ENVIRASTER",
    "direction": "output",
    "required": true,
    "description": "This is a reference to the output raster."

Comments (0) Number of views (1364) Article rating: No rating

Categories: IDL Blog | IDL Data Point




«March 2017»





















© 2017 Exelis Visual Information Solutions, Inc., a subsidiary of Harris Corporation