This is what I came up with until now (the code is still messy, but it threw out some results that look promising).
I have attached the geodatabases used for and create in the process.
Let me know if you run into problems.
I changed the name of 1 raster that wasn't listed in the csv to make it match the list of raster names in the csv (although that is wrong).
#-------------------------------------------------------------------------------
# Name: analyzeSpecies.py
# Purpose: raster overlay attribute analysis
#
# Author: Xander
#
# Created: 10-02-2015
#-------------------------------------------------------------------------------
def main():
import arcpy
import os
arcpy.env.overwriteOutput = True
# input settings
ws_ras = r"D:\Xander\GeoNet\SpeciesRaster\gdb\rasters.gdb"
ws_pol = r"D:\Xander\GeoNet\SpeciesRaster\gdb\polygons.gdb"
tbl_species = r"D:\Xander\GeoNet\SpeciesRaster\gdb\tables.gdb\SpeciesData"
fc_out = r"D:\Xander\GeoNet\SpeciesRaster\gdb\result.gdb\result"
# field names
fld_raster = "Name_SP"
fld_TW = "TW_SMUNDI_1"
fld_BL = "BL_SMUNDI_1"
fld_FL = "FL_SMUNDI_1"
lst_flds = [fld_TW, fld_BL, fld_FL]
# internal names
fcname_merge = "merged"
fcname_area = "aoi"
fcname_intersect = "intersect"
# path to tmp fc's
fc_merged = os.path.join(ws_pol, fcname_merge)
fc_area = os.path.join(ws_pol, fcname_area)
fc_intersect = os.path.join(ws_pol, fcname_intersect)
# dct with fields, input and output
dct_flds = {}
for fldname in lst_flds:
dct_flds[fldname] = "{0}_Mean".format(fldname[:2])
# create dictionary from species tabel (raster name vs stats)
flds_species = [fld_raster]
flds_species.extend(lst_flds)
dct_SpecieStats = {r[0]: [r[1], r[2], r[3]] for r in arcpy.da.SearchCursor(tbl_species, flds_species)}
# list polygons fc's
lst_polfcs = []
# list rasters and batch convert rasters to polygons
arcpy.env.workspace = ws_ras
lst_ras = arcpy.ListRasters()
cnt = 0
for ras in lst_ras:
cnt += 1
raster = os.path.join(ws_ras, ras)
fc_pol = os.path.join(ws_pol, ras)
print " - RasterToPolygon_conversion: {0}".format(ras)
# raster to polygon
arcpy.RasterToPolygon_conversion(raster, fc_pol, simplify="NO_SIMPLIFY")
# add fields to polygons
print " - AddDoubleField"
for fldname in lst_flds:
AddDoubleField(fc_pol, fldname)
# update values in polygon fc
print " - UpdateCursor(fc_pol..."
with arcpy.da.UpdateCursor(fc_pol, lst_flds) as curs:
for row in curs:
if ras in dct_SpecieStats:
lst_val = dct_SpecieStats[ras]
row[0] = correctValue(lst_val[0])
row[1] = correctValue(lst_val[1])
row[2] = correctValue(lst_val[2])
curs.updateRow(row)
else:
row[0] = -1
row[1] = -1
row[2] = -1
curs.updateRow(row)
# add to list to append or merge
lst_polfcs.append(fc_pol)
# sr
if cnt == 1:
sr = arcpy.Describe(fc_pol).spatialReference
# merge and append featureclasses
if cnt % 50 == 0:
if cnt == 50:
# merge
arcpy.Merge_management(lst_polfcs, fc_merged)
lst_polfcs = []
else:
# append
arcpy.Append_management(lst_polfcs, fc_merged, "NO_TEST")
lst_polfcs = []
# append the last fc's
print " - after loop, append or merge"
if cnt < 50:
# there is no merged fc, so use merge instead of append
if len(lst_polfcs) > 0:
arcpy.Merge_management(lst_polfcs, fc_merged)
lst_polfcs = []
else:
if len(lst_polfcs) > 0:
arcpy.Append_management(lst_polfcs, fc_merged, "NO_TEST")
lst_polfcs = []
# create area fc (extent of merged fc)
print "createRectangleExtendFC"
rectangle = createRectangleExtendFC(fc_merged, sr)
arcpy.CopyFeatures_management([rectangle], fc_area)
# execute intersect
print "Intersect_analysis"
arcpy.Intersect_analysis([fc_merged, fc_area], fc_intersect)
# create empty output fc
print "CreateFeatureclass_management"
fc_ws, fc_name = os.path.split(fc_out)
arcpy.CreateFeatureclass_management(fc_ws, fc_name, "POLYGON", spatial_reference=sr)
# add output fields
print "AddDoubleField output"
for fldname in lst_flds:
fldname_stat = dct_flds[fldname]
AddDoubleField(fc_out, fldname_stat)
# only use gridcode = 1
where = "{0} = {1}".format(arcpy.AddFieldDelimiters(fc_intersect, "gridcode"), 1)
# start analysis of unique polygons
print "start analysis of unique polygons"
dct_geo = {}
dct_att = {}
flds_stat = ["OID@", "SHAPE@"]
flds_stat.extend(lst_flds)
with arcpy.da.SearchCursor(fc_intersect, flds_stat, where_clause=where) as curs:
for row in curs:
oid = row[0]
polygon = row[1]
lst_val = [row[2], row[3], row[4]]
dct_att[oid] = lst_val
j = polygon.JSON
if j in dct_geo:
lst_oid = dct_geo
lst_oid.append(oid)
dct_geo = lst_oid
else:
dct_geo = [oid]
del row, curs
# calc stats
print "calculate statistics"
dct_res = {}
for i in range(0, 3):
for j, lst_oid in dct_geo.items():
lst_val = [dct_att[oid] for oid in lst_oid]
tot = sum(lst_val)
cnt = len(lst_val)
mean = tot / float(cnt)
oid_first = lst_oid[0]
if oid_first in dct_res:
lst_mean = dct_res[oid_first]
lst_mean.append(mean)
else:
lst_mean = [mean]
dct_res[oid_first] = lst_mean
# insert values
print "insert values into output fc"
flds_out = ["SHAPE@"]
for fldname in lst_flds:
flds_out.append(dct_flds[fldname])
# insert cursor
with arcpy.da.InsertCursor(fc_out, flds_out) as curs_out:
# search cursor
with arcpy.da.SearchCursor(fc_intersect, ("OID@", "SHAPE@")) as curs:
for row in curs:
oid = row[0]
if oid in dct_res:
polygon = row[1]
lst_mean = dct_res[oid]
curs_out.insertRow((polygon, lst_mean[0], lst_mean[1], lst_mean[2], ))
print "ready..."
def AddDoubleField(fc_pol, fldname):
if len(arcpy.ListFields(fc_pol, wild_card=fldname))==0:
arcpy.AddField_management(fc_pol, fldname, "DOUBLE")
def correctValue(val):
if type(val) == float:
return val
else:
return float(val)
def createRectangleExtendFC(fc, sr):
ext = arcpy.Describe(fc).extent
props = ["lowerLeft", "upperLeft", "upperRight", "lowerRight", "lowerLeft"]
arrPnts = arcpy.Array()
for crnr in props:
arrPnts.add(eval("ext.{0}".format(crnr)))
return arcpy.Polygon(arrPnts, sr)
if __name__ == '__main__':
main()
