Have you ever wondered or thought about what the strengths of ERDAS IMAGINE Spatial Modeler are and what it can do? ERDAS IMAGINE Spatial Modeler on its own provides hundreds of functions, algorithms and analytical routines that can easily be chained together into models to solve complex Geospatial problems. On the other hand ERDAS IMAGINE Spatial modeler by design is a highly extensible tool. In this article we will explore one such example of extensibility and see how other Geospatial softwares' capability can be used within ERDAS IMAGINE Spatial Modeler. As an example we would bring in the functionality from ESRI’s Python site package - ArcPy into Spatial Modeler.
Before we get into the actual model and the specifics of the usage of Python Script operator, let’s try to get a sneak peek into the actual workflow that is being executed in the form of a Spatial Model. Using stereo pair images from a satellite as an input, Automatic Point Measurement and Triangulation are performed and then a surface is generated. Since the output surface that is generated is in the form of a Point Cloud i.e., a LAS file but ArcPy commands which are executed using the Python Script operator, expect a raster DEM as input, the surface generated (LAS file) will be converted into a raster DEM. In order to make sure the raster DEM is having contiguous data without any voids an external reference DEM is used to fill in the voids if any. Since the ArcPy commands expect a raster DEM to have a defined “NoData” value the same is defined using the operators available in the ERDAS IMAGINE Spatial Modeler. This raster DEM is now fed as input to the ArcPy commands to generate a stream network in the form of a shapefile. Once the stream network from ArcPy is generated, this is fed as input to ERDAS IMAGINE Spatial Modeler feature based operators and further analysis like Measure Length and Buffer after applying an Attribute Filter can be performed. So this workflow is an example of how ERDAS IMAGINE Spatial Modeler can be used not only to exploit other Geospatial softwares’ capability but also consume the outputs from them and carry on further analysis within ERDAS IMAGINE.
Dissecting the model and its inputs:
ERDAS IMAGINE Spatial Modeler has a Python Script operator (check out the link to learn how to configure and use the Python Script operator) which can be used to conveniently execute any python script. We will be using this operator in this example.
Please note that here we are using ArcGIS for Desktop 10.4 and ERDAS IMAGINE 2016 v16.1 to demonstrate the functionality. Currently the Python versions used by ArcPy and ERDAS IMAGINE Spatial Modeler are different. The Python Script operator in Spatial Modeler runs on python 3.4.x, ArcPy for desktop runs on python 2.7.x. So in order to tackle the Python version difference we are providing a workaround by explicitly providing the Python path that ArcPy uses along with the Python script containing the ArcPy commands along with its input parameters. This can be achieved by creating a function within a Python script and giving this script as an input to the Python Script operator in ERDAS IMAGINE Spatial Modeler. The following is an example of the function..
import os import subprocess ARCPYDIR='C:/Python27/ArcGIS10.4/pythonw.exe' def ArcGISWorkflow_ArcPy(ArcGISpythonScript_File,InputDEMFile_File,OutputDirectory_Dir,ArcOutput_File): subprocess.call (ARCPYDIR + ' ' + ArcGISpythonScript_File + ' ' + InputDEMFile_File + ' '+ OutputDirectory_Dir + ' '+ ArcOutput_File) return ArcOutput_File
In the above example “ARCPYDIR” would be the installation path of pythonw.exe (residing inside Python27 folder) installed by ArcGIS, “ArcGISpythonScript_File” is the path of the python script containing the ArcPy commands and “InputDEMFile_File”, “OutputDirectory_Dir” and “ArcOutput_File” are the inputs to the “ArcGISpythonScript_File”.
The following is an example of the “ArcGISpythonScript_File” which contains the actual ArcPy commands. The left side code snippet executes the ArcPy commands Fill, Flow Direction, Flow Accumulation, Raster Calculator, Stream Link, Stream Order and Stream to Feature in the same order by taking the output of one command as an input to the next command. In the end the output stream network shapefile would be from Stream to Feature command.
import arcpy from arcpy.sa import * from arcpy import env import os import shutil import sys #Set local variables strInputDEMFile=sys.argv strZParameter='1921' strOutputDirectory=sys.argv strStreamToFeatureOutput=sys.argv #Output file declaration strFillOutput=strOutputDirectory+'/Fill.img' strFlowDirectionOutput=strOutputDirectory+'/FlowDirection.img' strFlowAccumulationOutput=strOutputDirectory+'/FlowAccumulation.img' strRasterCalculatorOutput=strOutputDirectory+'/StreamNetwork.img' strStreamLinkOutput=strOutputDirectory+'/StreamLink.img' strStreamOrderOutput=strOutputDirectory+'/StreamOrder.img' strStreamToFeatureOutput=strOutputDirectory+'/StreamToFeature.img' #Call Fill command and produce output arcpy.CheckOutExtension("Spatial") OutFill = Fill(strInputDEMFile,strZParameter ) OutFill.save(strFillOutput) inSurfaceRaster = strFillOutput # Call Flow Direction command and produce an output outFlowDirection =FlowDirection(inSurfaceRaster,"NORMAL") outFlowDirection.save(strFlowDirectionOutput) # Call Flow Accumulation command and produce an output outFlowAccumulation = FlowAccumulation(strFlowDirectionOutput,"","FLOAT") outFlowAccumulation.save(strFlowAccumulationOutput) # Call Raster Calculation command and produce an output InputRaster=strFlowAccumulationOutput arcpy.gp.RasterCalculator('Con("'+InputRaster+'" > 500,1)',strRasterCalculatorOutput) # Call Stream Link command and produce an output outStreamLink = StreamLink(strRasterCalculatorOutput,strFlowDirectionOutput) outStreamLink.save(strStreamLinkOutput) # Call Stream Order command and produce an output outStreamLink = StreamOrder(strStreamLinkOutput,strFlowDirectionOutput,"STRAHLER") outStreamLink.save(strStreamOrderOutput) # Call Stream to Feature command and produce an output StreamToFeature(strStreamOrderOutput,strFlowDirectionOutput,strStreamToFeatureOutput,"SIMPLIFY")
Alternatively the intermediate output files from the above commands can be deleted using the code snippet shown below.
Delete the Flow accumulation raster as this is an intermediate file if os.path.exists(strRasterCalculatorOutput): if os.stat(strRasterCalculatorOutput).st_size>0: arcpy.management.Delete(strFlowAccumulationOutput) #Delete the Fill Output as this is an intermediate file if os.path.exists(strFlowDirectionOutput): if os.stat(strFlowDirectionOutput).st_size>0: arcpy.management.Delete(strFillOutput) #Delete the raster calculator output as this is an intermediate file if os.path.exists(strStreamLinkOutput): if os.stat(strStreamLinkOutput).st_size>0: arcpy.management.Delete(strRasterCalculatorOutput) #Delete the Stream Link output as this is an intermediate file if os.path.exists(strStreamOrderOutput): if os.stat(strStreamOrderOutput).st_size>0: arcpy.management.Delete(strStreamLinkOutput) #Delete the Stream Order output as this is an intermediate file if os.path.exists(strStreamToFeatureOutput): if os.stat(strStreamToFeatureOutput).st_size>0: arcpy.management.Delete(strStreamOrderOutput)
The following is the screenshot of the model. The Photogrammetry operators are used to generate a surface from a stereo pair after performing Automatic Point Measurement and Triangulation respectively. Then the Point Cloud operators are used to convert the surface into a Raster DEM and then using this DEM as one of input to the Python Script operator (the other inputs being the actual Python script containing the ArcPy commands to which the DEM is fed as input, output folder and the output stream network shapefile) a stream network in the form of a shapefile is generated. This stream network shapefile is fed as an input to the Feature operators within ERDAS IMAGINE and a buffered output is generated after filtering the streams which are greater than 500 meters.
The following is the screenshot showing the output from the ArcPy commands on the left and on the right is the output upon further more processing the ArcPy commands output using the ERDAS IMAGINE Feature operators.
In the near future, we hope to simplify the integration.