**Integration with ArcGIS (ArcPy)** =================================== This section provides information about how to setup Geomatica's and ArcPy's python libraries so that they can be used together cleanly. Some basic error handling concepts are also exposed in this section for both Geomatica and ArcPy. At the bottom of this page, you will find a simple example that demonstrates how to combine ArcPy functions with Geomatica's python functions for creating powerful integrated workflows. `ArcPy Cookbook `_ Installing ArcGIS *64bit Background Geoprocesses* ------------------------------------------------- In order to use Geomatica's geospatial Python library with ArcPy, users, must install *ArcGIS' 64bit Background Geoprocesses*. This is required, because you cannot call Geomatica's 64bit modules with ArcGIS' native 32bit modules. Therefore, it is necessary to install ArcGIS' 64bit background geoprocesses. This section provides step-by-step instructions for installing the *ArcGIS 64bit Background Geoprocesses* **Steps:** **1.** Install ArcGIS 10.2.1 or later **2.** Acquire the ArcGIS for Desktop installation CD or the digital .iso file **Note:** You may need to install a virtual DVDrom (`Daemon Tools Lite `_) **3.** Click *Setup* next to *ArcGIS for Desktop Background Geoprocessing (64-bit)* .. image:: images/64bit_geoprocesses.png **4.** Follow the prompts and install with the default settings **5.** Verify your installation - Navigate to *C:\Python27* and make sure that the folder *ArcGISx6410.2* is found in the python27 directory .. image:: images/ArcPy64bit_verify.png Import ArcPy and Geomatica Python Libraries and verify configuration -------------------------------------------------------------------- This section provides instructions on how to verify that the *ArcGIS 64bit Background Geoprocessing* and Geomatica's Python library are properly installed and configured. **1.** Open a command prompt **2.** Navigate to and open the python executable in the *ArcGISx6410.2* directory, which will start python's interactive mode **Note:** When executing python scripts that call Geomatica's and ArcGIS' python libraries, it is required that you launch the script with the python executable in *C:\Python27\ArcGISx6410.2\python.exe* .. image:: images/ArcPy64bit_cmd.png **3.** Using Python's interactive mode, type the following to verify Geomatica's and ArcGIS' python libraries are correctly configured: .. code-block:: python #verify Geomatica's python library import pci #verify ArcGIS' python library import arcpy **Note:** If Geomatica's library is correctly configured, the statement, *"PCI Pluggable Framework environment successfully loaded."* will be printed to the terminal. If ArcPy's library is successfully configured nothing will be printed to the terminal. Both libraries will print error messages if there is a problem .. image:: images/ArcPy64bit_pci_verify.png Handling Geomatica and ArcPy Exceptions --------------------------------------- This example demonstrates how to handle exceptions that may result from incorrectly invoking a Geomatica function in Geomatica and ArcPy. `Click here `_ for more examples of error handling with ArcPy. .. code-block:: python from pci.pansharp2 import pansharp2 from pci.exceptions import PCIException import arcpy #Geomatic Python Try-except statement ms_file = r"c:\NewYorkk_multispectral.pix" #this file is incorrectly spelt and thus, does not exist ms_channels = [1,2,3,4] #pansharpen channels 1 to 4 (Red, green, blue and NIR) pan_file = r"c:\NewYork_panchromatic.pix" out_file = r"c:\NewYork_pansharpened.pix" try: pansharp2(ms_file, ms_channels, [], pan_file, [], out_file, [], "", [], "", "pix", "tiled512") except PCIException, e: print e except Exception, e: print e #ArcPy Try-except statement try: # Execute the Buffer tool # arcpy.Buffer_analysis("c:/transport/roads.shp", "c:/transport/roads_buffer.shp") except Exception as e: print e.message # If using this code within a script tool, AddError can be used to return messages # back to a script tool. If not, AddError will have no effect. arcpy.AddError(e.message) Combining Geomatica and ArcPy Functions --------------------------------------- This recipe uses Geomatica's library to automatically interpolate a DSM from LIDAR data. Then use ArcPy to calculate 3D visibility polygons from specific points of interest (POIs) .. code-block:: python import arcpy.mapping import locale import os import arcpy from pci.vdemingest import vdemingest from pci.vdemsetup import vdemsetup from pci.vdemint import vdemint from arcpy import env arcpy.CheckOutExtension("3D") #Licensing the ArcGIS 3D Analyst extension arcpy.env.overwriteOutput = True locale.setlocale( locale.LC_ALL, "" ) locale.setlocale( locale.LC_NUMERIC, "C" ) #********************************* #Start of Geomatica Code #********************************* #ingests the LIDAR point cloud (.las) into a working format (.pix) lidar_point_cloud = r"D:\PCI_Work\Projects\Open\Marketing\Python_Cookbook\tests\source_review\data_transform\lidar\lidar_cloud.las" lidar_return = "RETURNALL" #uses all lidar returns in the .las file pix_vector = r"D:\PCI_Work\Projects\Open\Marketing\Python_Cookbook\tests\source_review\data_transform\lidar\lidar_vector.pix" map_units = "UTM 17 D000" #vdemingest( lidar_point_cloud, "", "", "", lidar_return, [], "", # "", "", pix_vector, map_units, "", [] ) #setup the vector files for interpolation into a raster DEM index_file = r"D:\PCI_Work\Projects\Open\Marketing\Python_Cookbook\tests\source_review\data_transform\lidar\lidar_dem.txt" # output index file, used by VDEMINT dem_res = [1,1] #vdemsetup( lidar_point_cloud, index_file, "TIF", map_units, # dem_res, [], "", [], "" ) #interpolate the LIDAR point cloud maxiter = [20] # max 20 smoothing iterations memsize= [2048] # use 2048 MB (~2GB) of memory #vdemint( index_file, [1], [], pix_vector, [], [], "", maxiter, memsize ) output_lidar_dem = r"D:\PCI_Work\Projects\Open\Marketing\Python_Cookbook\tests\source_review\data_transform\lidar\lidar_dem.tif" #********************************* #Start of ArcGIS (ArcPy) Code #*********************************[ #Creating File Geodatabase out_dir = r"D:\PCI_Work\Projects\Open\Marketing\Python_Cookbook\tests\source_review\data_transform\lidar\arcpy_work\\" arcpy.CreateFileGDB_management(out_dir, 'pci_esri-solution.gdb') geodatabase_file = out_dir + os.sep + 'pci_esri-solution.gdb' #Set parameters and run visibility analysis vis_raster = geodatabase_file + os.sep + 'Visibility_raster' Analysis_type = "OBSERVERS" Import_Surface_Raster = output_lidar_dem Input_Observer_Feature = r"D:\PCI_Work\Projects\Open\Marketing\Python_Cookbook\tests\source_review\data_transform\lidar\points\poi.shp" Output_Raster = vis_raster AboGLevel = "" NODATA= [] arcpy.Visibility_3d(Import_Surface_Raster, Input_Observer_Feature, Output_Raster, AboGLevel, Analysis_type, "NODATA", "1", "FLAT_EARTH", "0.13", "5", "", "1.5", "", "100", "", "", "", "") print "Visibility Analysis complete!" # Convert Raster to Polygon in ArcGIS inRaster = vis_raster out_poly_folder = out_dir + os.sep + "polygons" outPolygons = out_poly_folder + os.sep + "visibility_poly.shp" field = "VALUE" # Execute RasterToPolygon arcpy.RasterToPolygon_conversion(inRaster, outPolygons, "NO_SIMPLIFY", field) print "Raster to Polygon Conversion Completed!"