Hi everyone,
I have around 1000 raster data layer. Every single one represents a distribution range of a species.
In the attribute table of each raster data layer are the body measurements of the specimen (there is only
one row per attribute table).
I am looking for a method in which I could calculate for each pixel the mean body size of the overlaying rasters.
Is there something like that? It would be best if I could repeat that for severals columns in the attribute table.
Is there a more or less simple solution for this problem? Haven't found anything till now..
Many thanks,
Nicolas
Solved! Go to Solution.
It may take some time to work out the idea I have. I started on another method, but maybe I will have to adapt this based on the data provided. I will get back to you as soon as I have something to show.
Kind regards, Xander
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_geolst_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()
Hi Xander,
I don't have experience with python. So it will take a bit to figure it out for me. But after all problems I had with the model builder I think it is anyway a must to know python. I let you know when I run in problem thanks.
Many thanks,
Nicolas