This is the first of several posts on performing Kappa and/or chi-square statistics. These tests are commonly used to see whether remote sensing training site data matches up with that determined via image classification. This first script prepares for subsequent analysis.
'''
PrepareData.py
Purpose: Prepares data for Kappa and Chi-square analysis
Author:
Dan Patterson
Dept of Geography and Environmental Studies
Carleton University, Ottawa, Canada
Dan_Patterson@carleton.ca
Date created: Jan 2010
Last modified: June 2013
Requires:
obs_FC a shapefile of observed points with a classification
buffer_dist a buffer distance
expected_lyr a layer file of expected values, either a shapefile or a grid.
output_FC an output shapefile
'''
#Import the basic modules
import os, sys, string
import arcpy
obs_FC = sys.argv[1]
buffer_dist = float(sys.argv[2])
expected_lyr = sys.argv[3]
output_FC = sys.argv[4]
temp_FC = os.path.dirname(sys.argv[0]) + "/temp.shp"
temp2_FC = os.path.dirname(sys.argv[0]) + "/temp2.shp"
expected_type = arcpy.Describe(expected_lyr).datasetType
#
#Check to see if the expected data layer is a raster, if it is convert it to a shapefile
#Buffer the observed layer, then perform an intersection
if expected_type == "RasterDataset":
arcpy.AddMessage("\nConverting to temporary shapefile: " + str(expected_lyr))
arcpy.RasterToPolygon_conversion(expected_lyr, temp2_FC, "NO_SIMPLIFY", "VALUE")
arcpy.AddMessage("\nIntersecting layers to produce: " + str(output_FC))
arcpy.AddMessage("\nDeleting temporary files")
if buffer_dist > 0.0:
arcpy.Buffer_analysis(obs_FC, temp_FC, buffer_dist, "#", "#", "NONE")
arcpy.Intersect_analysis([temp_FC, temp2_FC], output_FC, "NO_FID")
arcpy.Delete_management(temp_FC)
else:
arcpy.Intersect_analysis([obs_FC, temp2_FC], output_FC, "NO_FID")
arcpy.Delete_management(temp2_FC)
else:
arcpy.AddMessage("\nIntersecting layers to produce: " + str(output_FC))
arcpy.AddMessage("\nDeleting temporary files")
if buffer_dist > 0.0:
arcpy.Buffer_analysis(obs_FC, temp_FC, buffer_dist, "#", "#", "NONE")
arcpy.Intersect_analysis([temp_FC, expected_lyr], output_FC, "NO_FID")
arcpy.Delete_management(temp_FC)
else:
arcpy.Intersect_analysis([obs_FC, expected_lyr], output_FC, "NO_FID")
#
#Add an area field and calculate area.
if buffer_dist > 0.0:
arcpy.AddField_management(output_FC, "Area", "DOUBLE", 16, 6)
arcpy.CalculateField_management(output_FC, "Area", "!shape.area!", "PYTHON_9.3")
arcpy.AddMessage("\nProcessing complete....\n")
No comments:
Post a Comment