wiki:Invoking MGET tools programmatically
Last modified 5 years ago Last modified on 07/25/12 10:03:13

Invoking MGET tools programatically

MGET is, in essence, a library of functions. The functions are implemented in various programming languages, mainly Python, R, C++, and MATLAB, but are exposed uniformly via several widely-used interface technologies. Depending on your favored programming environment, you may invoke a given function:

  • From ArcGIS as a geoprocessing tool
  • From Python as a classmethod of a class
  • From many other languages as a method of a COM Automation class

This document illustrates all of the common invocation scenarios. If you do not see your scenario and require assistance, please email Jason Roberts (jason.roberts@duke.edu).

Each scenario performs the same task, conversion of a Scientific Data Set stored in HDF file to an ArcGIS raster, using a different programming language. The SDS, stored in 199001.s04m1pfv50-sst-16b.hdf, represents the global monthly sea surface temperature for January 1990 from the AVHRR sensor on the NOAA POES satellites, and is published by the NOAA NODC 4 km AVHRR Pathfinder Project. You can read more about this data in the example Converting HDFs to ArcGIS rasters.

For more information on MGET's internal architecture, please see our paper (Roberts et al., 2010).

MGET does not have a formal API yet

Although we explain how to call MGET's functions below and provide some API-style documentation elsewhere, MGET does not provide a formal API at this time. All objects and functions are subject to change across MGET releases, with no warning from us. We do not suggest you build anything critical on top of MGET unless you are willing to tolerate that. As we evolve MGET's design, we may eventually commit to a formal API. When that happens, we will provide some guarantee of stability, proper notification of changes, and so on.

Now, back to the example...

Invoking MGET from ArcGIS

The MGET functions that may be invoked from ArcGIS appear as geoprocessing tools in the Marine Geospatial Ecology Tools toolbox. You can invoke MGET tools just like any other geoprocessing tools: from geoprocessing models, from the ArcGIS command line prompt, and from geoprocessing scripts.

From geoprocessing models

source:/WikiFiles/Examples/InvokingMGETToolsProgrammatically/GeoprocessingModel.png

From the ArcGIS command line

To find out the command-line name of an MGET tool (e.g. HDFSDSToArcGISRaster_GeoEco), open the help page for the tool. For example: Convert SDS in HDF to ArcGIS Raster.

source:/WikiFiles/Examples/InvokingMGETToolsProgrammatically/ArcGISCommandLine.png

From geoprocessing scripts

Important: Due to a bug in ArcGIS 9.3 and 10, acknowledged by ESRI, the technique of invoking MGET tools through the ArcGIS geoprocessor that is described immediately below may cause Python development tools such as PythonWin and IDLE to crash. To avoid this, try the method described in the next section, Invoking MGET from Python.

Important: This may not work with ArcGIS 10.0 if you access the MGET toolbox by calling arcpy.ImportToolbox(). Instead, you need to import the toolbox the same way you did with 9.3, using the arcgisscripting module, as shown below. Please see #534 for more information.

To invoke MGET tools from your own geoprocessing scripts, initialize the geoprocessor object using the normal procedure, load the MGET toolbox, and invoke the tool. Here is a Python example:

# Before executing the code below, you must initialize the gp object using the normal procedure.

gp.AddToolbox('C:\\Program Files\\GeoEco\\ArcGISToolbox\\Marine Geospatial Ecology Tools.tbx')

gp.HDFSDSToArcGISRaster_GeoEco('C:\\temp14\\hdf\\199001.s04m1pfv50-sst-16b.hdf', 'C:\\temp14\\Pathfinder_Example_Raster\\sst199001',
                               'sst', -180, -90, 0.0439453125, 0, None, None, None, None,
                               "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]",
                               None, None, None, None, None, None, 'float(inputRaster) * 0.075 - 3.0', True)

MGET is implemented as a Python package and is installed in your Python installation directory, but the toolbox itself appears within the Windows C:\Program Files\GeoEco directory.

To find out the name of the MGET tool (e.g. HDFSDSToArcGISRaster_GeoEco), open the help page for the tool. For example: Convert SDS in HDF to ArcGIS Raster.

When you invoke the tool from Python through the ArcGIS geoprocessor, pass None or "#" for optional parameters that you do not want to specify. The tool will use the default values for the parameters. For more information on passing parameters to tools, see the ArcGIS documentation.

Due to how ArcGIS is designed, MGET log messages may not be reported when you invoke MGET tools from your own geoprocessing script. If you're writing your script in Python and wish to see the MGET log messages, you can invoke MGET directly through Python rather than indirectly through the geoprocessor. Please see below for more information.

Invoking MGET from Python

Python is the best language for programming with MGET because it provides complete access to all MGET functionality. Regardless of what language is used to implement a given MGET function, it will always be exposed as a Python classmethod. We try to also expose every function as an ArcGIS geoprocessing tool and a method of a COM Automation class, but sometimes technical issues prevent this. By using Python, you will avoid this problem.

MGET is implemented as a collection of classes in the GeoEco Python package. The usual way to call an MGET function from Python is to import the class that contains it and call the classmethod:

>>> from GeoEco.DataManagement.HDFs import HDF
>>> HDF.SDSToArcGISRaster('C:\\temp14\\hdf\\199001.s04m1pfv50-sst-16b.hdf', 'C:\\temp14\\Pathfinder_Example_Raster\\sst199001', 'sst', -180, -90, 0.0439453125, 0,
...                       coordinateSystem="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]",
...                       mapAlgebraExpression='float(inputRaster) * 0.075 - 3.0', buildPyramids=True, overwriteExisting=True)
>>> 

MGET includes complete Python Reference Documentation which documents every module (e.g. GeoEco.DataManagement.HDFs), class (e.g. HDF), and method (e.g. SDSToArcGISRaster).

MGET functions report their progress using Python's logging infrastructure. This is turned off by default, but you can turn it on by calling Logger.Initialize before calling other MGET functions:

>>> from GeoEco.Logging import Logger
>>> Logger.Initialize()
>>> 
>>> from GeoEco.DataManagement.HDFs import HDF
>>> HDF.SDSToArcGISRaster('C:\\temp14\\hdf\\199001.s04m1pfv50-sst-16b.hdf', 'C:\\temp14\\Pathfinder_Example_Raster\\sst199001', 'sst', -180, -90, 0.0439453125, 0,
...                       coordinateSystem="GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]",
...                       mapAlgebraExpression='float(inputRaster) * 0.075 - 3.0', buildPyramids=True, overwriteExisting=True)
2008-05-28 13:54:26,634 INFO Converting SDS "sst" in HDF file C:\temp14\hdf\199001.s04m1pfv50-sst-16b.hdf to ArcGIS raster C:\temp14\Pathfinder_Example_Raster\sst199001...
>>> 

You can write ArcGIS geoprocessing scripts using MGET. The following example converts a specified SST HDF from the NOAA NODC 4 km AVHRR Pathfinder Project to an ArcGIS raster and masks all pixels that do not meet a given minimum quality flag. The pixel quality data comes from a second HDF file.

source:/WikiFiles/Examples/InvokingMGETToolsProgrammatically/PathfinderExample.png

# Initialize MGET logging and ensure that log messages are reported to ArcGIS.

from GeoEco.Logging import Logger
Logger.Initialize(activateArcGISLogging=True)

# Catch and log any exception that is raised.

try:
    # Initialize the ArcGIS geoprocessor.

    from GeoEco.ArcGIS import GeoprocessorManager
    GeoprocessorManager.InitializeGeoprocessor()
    gp = GeoprocessorManager.GetWrappedGeoprocessor()

    # Get the input parameters.

    sstHDF = gp.GetParameterAsText(0)
    qualHDF = gp.GetParameterAsText(1)
    minQual = int(gp.GetParameterAsText(2))
    outputRaster = gp.GetParameterAsText(3)

    # Create a temporary directory to hold temporary rasters. This directory
    # will be automatically deleted when the script exits.

    from GeoEco.DataManagement.Directories import TemporaryDirectory
    tempDir = TemporaryDirectory()

    # Convert the input HDFs to rasters.

    import os
    from GeoEco.DataManagement.HDFs import HDF
    cs = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]"
    HDF.SDSToArcGISRaster(sstHDF, os.path.join(tempDir.Path, 'sst'), 'sst', -180, -90, 0.0439453125, 0, coordinateSystem=cs, mapAlgebraExpression='float(inputRaster) * 0.075 - 3.0')
    HDF.SDSToArcGISRaster(qualHDF, os.path.join(tempDir.Path, 'qual'), 'qual', -180, -90, 0.0439453125, 0, coordinateSystem=cs)

    # Mask cells of the SST raster that do not meet the minimum quality level.

    Logger.Info('Masking SST pixels where the quality flag is less than %i...' % minQual)

    gp.CheckOutExtension('spatial')
    gp.SetNull_sa(os.path.join(tempDir.Path, 'qual'), os.path.join(tempDir.Path, 'sst'), os.path.join(tempDir.Path, 'masked'), 'Value >= %i' % minQual)

    # Copy the masked raster to the output location.

    from GeoEco.DataManagement.ArcGISRasters import ArcGISRaster
    ArcGISRaster.Copy(os.path.join(tempDir.Path, 'masked'), outputRaster, overwriteExisting=bool(gp.OverwriteOutput))
                      
except:
    Logger.LogExceptionAsError()
    raise

Example output:

Executing (Convert and Mask Pathfinder SST HDF): Script1 C:\temp14\hdf\199001.s04m1pfv50-sst-16b.hdf C:\temp14\hdf\199001.m04m1pfv50-qual.hdf 7 C:\temp14\Rasters\sst199001
Start Time: Wed May 28 16:31:59 2008
Running script Script1...
Converting SDS "sst" in HDF file C:\temp14\hdf\199001.s04m1pfv50-sst-16b.hdf to ArcGIS raster C:\Temp\GeoEcoTemp_Jason\tmperitds\sst...
Converting SDS "qual" in HDF file C:\temp14\hdf\199001.m04m1pfv50-qual.hdf to ArcGIS raster C:\Temp\GeoEcoTemp_Jason\tmperitds\qual...
Masking SST pixels where the quality flag is less than 7...
Copying ArcGIS raster C:\Temp\GeoEcoTemp_Jason\tmperitds\masked to C:\temp14\Rasters\sst199001...
Completed script Script1...
Executed (Convert and Mask Pathfinder SST HDF) successfully.
End Time: Wed May 28 16:39:06 2008 (Elapsed Time: 7 minutes 7 seconds)

Some important notes about this example script:

  • When it calls Logger.Initialize, it sets activateArcGISLogging=True. This must be done for log messages to appear in the ArcGIS geoprocessing window.
  • It obtains the geoprocessor object by calling GeoprocessorManager.InitializeGeoprocessor and GeoprocessorManager.GetWrappedGeoprocessor. This instantiates the geoprocessor object and returns a reference to MGET's internal instance of it. Obtaining the geoprocessor this way provides several advantages:
    • MGET automatically instantiates the geoprocessor using win32com.client.Dispatch or arcgisscripting.create, depending on which version of ArcGIS is on the machine, allowing the script to run on both ArcGIS 9.1 and 9.2+ with no changes.
    • MGET automatically handles some incompatibilities between the win32com.client.Dispatch and the arcgisscripting.create geoprocessor.
    • MGET logs all calls to the geoprocessor (although you must enable verbose logging to see these messages).

But note that when MGET calls arcgisscripting.create, it always calls it like this: arcgisscripting.create(), never like this: arcgisscripting.create(9.3). Therefore, if you have ArcGIS 9.3 installed, the MGET geoprocessor will use the 9.2-style interface, not the 9.3-style interface. We chose this design to simplify the job of making MGET simultaneously compatible with 9.1, 9.2, and 9.3. So if you call the geoprocessor instantiated by MGET, you must also use the 9.2-style methods and properties of the geoprocessor.

  • It wraps all processing in a try block and calls Logger.LogExceptionAsError in the except block. If an exception is raised, even by a non-MGET function, this exception handler will log debugging information (you must enable verbose logging to see these messages).

  • It doesn't do any parameter validation. MGET functions perform a lot of parameter validation and try to report sensible errors. Where possible, you can rely on MGET functions to validate input parameters for you. This example script relies on HDF.SDSToArcGISRaster to validate the input HDF files and on ArcGISRaster.Copy to validate the output raster.

Also:

  • You cannot pass '#' to omit a parameter. Instead, just omit it from the function call. Of course, if you want to omit a parameter but specify one that comes after it, you should omit the first parameter but then specify all subsequent parameters as keyword arguments (see the Python documentation for more information).
  • When passing a parameter that is a list, do not provide a semicolon-delimited list like you do when calling tools through the geoprocessor, such as 'aaaaa;bbbbb;ccccc'. Instead, pass a Python list: ['aaaaa', 'bbbbb', 'ccccc'].

Invoking MGET through Windows COM Automation

Many MGET functions are exposed as methods of Windows COM Automation classes. In the example presented here, the SDSToArcGISRaster function is a method of the class that has ProgID GeoEco.HDF. Most programming languages that run on Windows provide some mechanism for instantiating COM Automation classes and invoking their methods. Here are some examples.

Important note for ArcGIS 9.3 users: Due to an apparent bug in ArcGIS 9.3, if your program uses COM Automation to invoke an MGET tool that itself invokes ArcGIS, your program may crash when it exits. Please see ticket #291 for more information. As far as I can tell, this problem will affect all of the examples shown below.

VBScript

As in the Python example above, this example initializes the MGET logging system before converting the HDF:

Set logger = WScript.CreateObject("GeoEco.Logger")
logger.Initialize()

Set hdf = WScript.CreateObject("GeoEco.HDF")
hdf.SDSToArcGISRaster "C:\\temp14\\hdf\\199001.s04m1pfv50-sst-16b.hdf", "C:\\temp14\\Pathfinder_Example_Raster\\sst199001", "sst", -180, -90, 0.0439453125, 0, False, False, False, False, "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", Missing, Missing, Missing, Missing, Missing, Missing, "float(inputRaster) * 0.075 - 3.0", True, True

You can see the script's output by running it from a CMD shell:

C:\temp14>cscript test1.vbs
Microsoft (R) Windows Script Host Version 5.6
Copyright (C) Microsoft Corporation 1996-2001. All rights reserved.

2008-05-28 17:29:06,605 INFO Converting SDS "sst" in HDF file C:\temp14\hdf\199001.s04m1pfv50-sst-16b.hdf to ArcGIS raster C:\temp14\Pathfinder_Example_Raster\sst199001...

C:\temp14>

To omit an optional parameter in VBScript, pass Missing for the parameter.

JScript

The JScript code is nearly identical to the VBScript code:

var logger = WScript.CreateObject("GeoEco.Logger");
logger.Initialize();

var hdf = WScript.CreateObject("GeoEco.HDF");
hdf.SDSToArcGISRaster("C:\\temp14\\hdf\\199001.s04m1pfv50-sst-16b.hdf", "C:\\temp14\\Pathfinder_Example_Raster\\sst199001", "sst", -180, -90, 0.0439453125, 0, false, false, false, false, "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", null, null, null, null, null, null, "float(inputRaster) * 0.075 - 3.0", true, true);

To omit an optional parameter in JScript, pass null for the parameter. The output looks the same if you run the script from a CMD shell.

R

First install the rcom package so you can call COM Automation components from R. Then:

> library(rcom)
> hdf <- comCreateObject("GeoEco.HDF")
> comInvoke(hdf, "SDSToArcGISRaster", "C:\\temp14\\hdf\\199001.s04m1pfv50-sst-16b.hdf", "C:\\temp14\\Pathfinder_Example_Raster\\sst199001", "sst", -180, -90, 0.0439453125, 0, FALSE, FALSE, FALSE, FALSE, "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", NULL, NULL, NULL, NULL, NULL, NULL, "float(inputRaster) * 0.075 - 3.0", TRUE, TRUE)
NULL
> 

This example does not initialize the MGET logging system because the rcom package does not capture stdout and echo it back to the R console. If you initialized the logging system, the messages would be reported but you would never see them.

To omit an optional parameter in R, pass NULL for the parameter.

MATLAB

In MATLAB, the actxserver function provides access to COM Automation servers:

>> hdf = actxserver('GeoEco.HDF');
>> hdf.SDSToArcGISRaster('C:\\temp14\\hdf\\199001.s04m1pfv50-sst-16b.hdf', 'C:\\temp14\\Pathfinder_Example_Raster\\sst199001', 'sst', -180, -90, 0.0439453125, 0, false, false, false, false, 'GEOGCS[''GCS_WGS_1984'',DATUM[''D_WGS_1984'',SPHEROID[''WGS_1984'',6378137.0,298.257223563]],PRIMEM[''Greenwich'',0.0],UNIT[''Degree'',0.0174532925199433]]', [], [], [], [], [], [], 'float(inputRaster) * 0.075 - 3.0', true, true);
>>

This example does not initialize the MGET logging system because the actxserver infrastructure does not capture stdout and echo it back to the MATLAB console. If you initialized the logging system, the messages would be reported but you would never see them.

To omit an optional parameter when calling COM Automation components from MATLAB, pass [] for the parameter.

Invoking MGET through Windows COM

The MGET Windows COM Automation classes are dual interface classes, allowing them to be instantiated through traditional COM. Through this mechanism, MGET functions can be invoked via a "vtable call" from early-bound languages that cannot easily call COM Automation interfaces. These languages include C++, C#, and Fortran. To facilitate this process, the MGET installer registers a COM type library named "GeoEco X.Y Type Library", where X and Y are the MGET major and minor version numbers.

C#

Before compiling this code in Microsoft Visual Studio, you must add a COM reference to the component named "GeoEco X.Y Type Library", where X and Y are the MGET version numbers.

using System;
using GeoEco;

namespace Example
{
    class Program
    {
        static void Main(string[] args)
        {
            Logger logger = new Logger();
            logger.Initialize(false, null);

            HDF hdf = new HDF();
            hdf.SDSToArcGISRaster("C:\\temp14\\hdf\\199001.s04m1pfv50-sst-16b.hdf", "C:\\temp14\\Pathfinder_Example_Raster\\sst199001",
                                  "sst", -180, -90, 0.0439453125, 0, Type.Missing, Type.Missing, Type.Missing, Type.Missing,
                                  "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]",
                                  Type.Missing, Type.Missing, Type.Missing, Type.Missing, Type.Missing, Type.Missing, "float(inputRaster) * 0.075 - 3.0", true, true);
        }
    }
}

One problem with C# is that you cannot omit parameter values. For example, in the call to logger.Initialize, false and null are explicitly passed, even though these are the default values for these parameters. When methods have many optional parameters, as is often the case with MGET, this can be very cumbersome.

This problem can be broken down into two issues: First off, it is annoying that you have to specify values for all of the parameters. Secondly, it is dangerous that you have to choose a value for a parameter that you wish to leave alone. This second issue can sometimes be addressed by passing Type.Missing as the value. The compiler will not always allow you to do this (for reasons I'm not going to go into). But when it does, the MGET function will use the default value for the parameter. This allows your code to remain agnostic about what default MGET uses, and protects you if MGET decides to change the default in the future.