New in IDL 8.6: IDLTasks

Author: Brian Griglak

ENVITasks have been around for a few years now, but they didn’t help those of you who only have IDL and not ENVI as well.  In IDL 8.6 we are happy to announce that the ENVITask concept has been carried over to work in pure IDL without ENVI.  Many of the features are the same, but there are slight differences that I want to point out in this post.

The first and most obvious is how you construct an IDLTask.  Rather than launch ENVI and then call the ENVITask()function, you can call the IDLTask()function at any point in time.  As the help points out, you can pass in either a scalar string that is the name of the task you want or its filename, or you can pass in a Hash object that is the task definition (what you would get if you called JSON_Parse() on a task file).  If you pass in a string that isn’t a filename (relative or absolute), then it is treated as a task name.  The function will then look at !PATH and search for all .task files in each folder listed there.  It catalogs all the .task files it finds, but if there are multiple folders with the same base name only the first is recognized(just like how IDL handles multiple .pro or .sav files with the same base name found in !PATH).  The list of .task files is filtered to those with the correct IDLTask schema, which currently is only “idltask_1.0”.  This way we don’t accidentally pick up an ENVITask files and cause confusion.  If a.task file with the same base name as the requested task name is found, it is used as the task definition.  If no exact match is found, but partial matches exist, then helpful error messages are returned telling you about the name(s) that partially match, so you can correct your code.  I should point out that the current working directory (which can be retrieved by calling CD with the CURRENT keyword) is searched before any of the folders in !PATH, so that can affect the behavior of IDLTask().

The “idltask_1.0” task schema used for IDLTasks in IDL 8.6 is very similar to the “envitask_3.0” schema used by ENVITasks in ENVI 5.4.  The notable exception is that the TYPE property of your parameters won’t understand ENVI class types like ENVIRaster.  But all the basic datatypes available in IDL are supported by IDLTasks – strings, Booleans, and numbers, as well as List, Hash, OrderedHash, and Dictionary.

Another difference is how you interact with IDLTasks on GSF as opposed to ENVITasks.  The service endpoint for IDLTasks will be http://hostname:port/ese/services/IDL,while the ENVITasks use http://hostname:port/ese/services/ENVI.  The different endpoints are used to discriminate between the requests that should use the IDLTask() function vs the ENVITask() function to load the requested task.

Easy GSF deployment is one of the primary reasons you would want to build IDLTasks in the first place. If you have IDL functions or procedures that you are used to calling directly, then you are probably wondering why you would want to wrap them in an IDLTask.  As a C++ developer in a previous life, I appreciated the type safety that C++ requires, so I also appreciate the parameter validation that IDLTasks provide.  When developing your custom IDLTask, you will have to spend some time thinking about what the inputs and outputs are for your code, but once you do that you won’t need to worry about writing lots of input validation code, the IDLTask framework will take care of that for you.  The IDLTasks are also self-documenting like ENVITasks, so if someone else hands you a .task file and .sav file, you can load the task and then learn all about the parameter names, their types, cardinalities, and hopefully even descriptions. All of this information makes it possible to deploy your algorithms on GSF for running in the cloud, with all the same introspection capabilities over the REST endpoint.  Alternatively, you can set up batch processing using some sort of folder watch capability to spawn IDLTaskEngine instances to automatically run your code on each file that appears on your system.


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

Categories: IDL Blog | IDL Data Point





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 (418) 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 (658) 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 (604) 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 (908) Article rating: No rating

Categories: IDL Blog | IDL Data Point


12345678910 Last




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