This is an exemplary code for using ESRI ArcGIS python packages in scripts. Similar scripts are possible with the use of QGIS libraries or similar.
As a starter, you should always use the model builder from ArcGIS to set the general structure of your code, including environment settings. Afterwards, you add the path to arcpy to your environment and code for automatization.
The following example shows a setup to produce cross-validation based goodness of fit measures on a simple IDW scheme.
NOTES:
(1) The following code can only be run on a machine that has an installation of ArcGIS. There is no preview in Jupyter Notebooks.
(2) Sometimes, it is a bit tricky to find out if files are assigned as relative or absolute paths. Be sure to have that in mind.
(3) If you want the data sets, needed for this example, contact me: b.mueller@iggf.geo.uni-muenchen.de
In [6]:
# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------------
# myIDW.py
# Created on: 2017-06-23 12:23:24.00000
# (generated by ArcGIS/ModelBuilder)
# Description:
# ---------------------------------------------------------------------------
# make arcpy module available
import sys
# this is depending on your installation
sys.path.extend(['C:\\Program Files (x86)\\ArcGIS\\Desktop10.5\\bin'])
sys.path.extend(['C:\\Program Files (x86)\\ArcGIS\\Desktop10.5\\arcpy'])
sys.path.extend(['C:\\Program Files (x86)\\ArcGIS\\Desktop10.5\\ArcToolbox\\Scripts'])
import arcpy
# additional packages
import os
import shapefile
import itertools
# make the spatial analyst available (or any other extension)
arcpy.CheckOutExtension("Spatial")
# set YOUR working directory (makes it easier to read data)
os.chdir("C:\\Users\\b.mueller\\Desktop\\Python ArcGIS Example\\ArcGIS\\")
# optional: option for overwriting files by arcpy functions
arcpy.env.overwriteOutput = 1
# Local variables:
# Input:
myplaces = "myplaces.shp"
# Output:
sm_inter = "C:\\Users\\b.mueller\\Desktop\\Python ArcGIS Example\\ArcGIS\\sm_inter"
# Set Geoprocessing environments (done by model builder if set)
# the raster, you want the results to be snapped to
arcpy.env.snapRaster = "mysnap"
# the extent, you want to be covered
arcpy.env.extent = "5,86599881313816 47,2703623267204 15,0377433356746 55,0573747014351"
# read shapefile:
M=shapefile.Reader(myplaces)
# get a list of combinations for either validation or interpolation (here: validation)
# first, get all tuples that make use of n records from the shapefile
n=1
All_Comb=list(itertools.permutations(range(M.numRecords),n))
# then, get rid of all redundant tuples
All_Comb=list(set([tuple(sorted(x)) for x in All_Comb]))
# NOTE: Which n you can use depends strongly on memory and data size! (workaround: make these lines more efficient)
# now, loop through the combinations
for combs in All_Comb:
# set the actual combination as a validation subset (or interpolation, as you like)
Valid_combs=[M.records()[i] for i in combs]
Valid_shape=[M.shapes()[i] for i in combs]
# calculate the left over combinations for interpolation (or validation, s.a.)
left_combs=range(M.numRecords)
for c in combs:
left_combs.remove(c)
# set the left overs as interpolation subset
Inter_combs=[M.records()[i] for i in left_combs]
Inter_shape=[M.shapes()[i] for i in left_combs]
# define unique filenames for the two subsets
filename_of_dbf_V = "C:\\Users\\b.mueller\\Desktop\\Python ArcGIS Example\\ArcGIS\\sub\\SHP_V_" + \
"".join([str(c) for c in list(combs)])
filename_of_dbf_I = "C:\\Users\\b.mueller\\Desktop\\Python ArcGIS Example\\ArcGIS\\sub\\SHP_I_" + \
"".join([str(c) for c in list(combs)])
# initialize the writer for validation data
M_valid_w=shapefile.Writer(M.shapeType)
# field types can be reused from original data...
for f in M.fields:
M_valid_w.field(f[0],f[1],f[2],f[3])
# ...as well as the bounding box
M_valid_w.bbox = lambda: M.bbox
# now write the coordinates into the data set...
for P in Valid_shape:
a,b = P.points[0]
M_valid_w.point(a,b)
# ... and the actual values
# (if there's more than ID and value, you'll need to add more "letters" (a,b,c,etc...))
for R in Valid_combs:
a,b = R
M_valid_w.record(a,b)
# finally, write the validation data
M_valid_w.save(filename_of_dbf_V)
# do the same for the interpolation
# PRO-TIPP: make all this a function of M and the subset samples
M_inter_w=shapefile.Writer(M.shapeType)
for f in M.fields:
M_inter_w.field(f[0],f[1],f[2],f[3])
M_inter_w.bbox = lambda: M.bbox
for P in Inter_shape:
a,b = P.points[0]
M_inter_w.point(a,b)
for R in Inter_combs:
a,b = R
M_inter_w.record(a,b)
M_inter_w.save(filename_of_dbf_I)
# filenames for the produced outputs
# the interpolated raster
filename_of_dbf_raster = "C:\\Users\\b.mueller\\Desktop\\Python ArcGIS Example\\ArcGIS\\sub\\R_" + \
"".join([str(c) for c in list(combs)])
# the zonal statistic files
filename_of_dbf_ZonalSt = "C:\\Users\\b.mueller\\Desktop\\Python ArcGIS Example\\ArcGIS\\sub\\Z_" + \
"".join([str(c) for c in list(combs)])
# the resulting excel file which can be used for rmse calculation, aggregation, etc.
filename_of_dbf_xls = "C:\\Users\\b.mueller\\Desktop\\Python ArcGIS Example\\ArcGIS\\sub\\X_" + \
"".join([str(c) for c in list(combs)]) + ".xls"
# Process: IDW
arcpy.gp.Idw_sa(filename_of_dbf_I + ".shp", "sm", filename_of_dbf_raster, "0,01", "5", "VARIABLE 12", "")
# Process: Zonal Statistics as Table
arcpy.gp.ZonalStatisticsAsTable_sa(filename_of_dbf_V + ".shp", "CID", filename_of_dbf_raster,
filename_of_dbf_ZonalSt, "DATA", "MEAN")
# Process: Table to CSV
arcpy.TableToExcel_conversion(filename_of_dbf_ZonalSt,filename_of_dbf_xls, "NAME", "CODE")
# finally, give back the license (for security and multi-use reasons)
arcpy.CheckInExtension("Spatial")