I am currently working on a Python Toolbox for ArcMap. This Toolbox needs to work without any extensions (Spatial Analyst, etc).
One step is to summarise (multi-band) raster values within polygons of a Shapefile or Feature Layer. Ideally, raster values are averaged within the polygon areas, taking into account that some raster pixels are only partly covered by a polygon.
The final result is a table with a Polygon-ID column and columns containing the average values for each raster band.
I tried the following (which is not the ideal attempt, since it does not account for fraction covered of raster cells):
import arcpy
arcpy.conversion.RasterToPoint(raster_path, "in_memory/pts")
arcpy.analysis.SpatialJoin(polygon_path, "in_memory/pts", "in_memory/pgn", "JOIN_ONE_TO_MANY")
which failed with some Runtime error: "RasterToPoint raise ExecuteError: ERROR 001143: The background server threw an exception" (translated).
Anyway, since this was not the ideal approach it might not even be worth solving the issue but instead, I hope somebody knows a better approach to this task. I found it trivial to code this in R or QGIS, so I hope there is also an easy and effective way to solve this basic task in an ArcMAP (basic) Python Toolbox.
Edit
I thought about iterating through the polygons of the Feature Class, creating points at raster cell centroids, selecting the points within the polygon and then using those to access raster values. The code I've written so far is as follows:
import os
import arcpy as ap
import numpy as np
import pandas as pd
import tempfile as tmp
from ast import literal_eval
ap.env.overwriteOutput = True
arcpy.env.addOutputsToMap = False
# Functions
## Extract raster information
def get_raster_info(raster):
vars = ["top", "bottom", "left", "x_size", "y_size", "n_col", "n_row", "wavelengths"]
fieldcodes = ["TOP", "BOTTOM", "LEFT", "CELLSIZEX", "CELLSIZEY", \
"COLUMNCOUNT", "ROWCOUNT", "WAVELENGTH"]
out = dict()
for var, code in zip(vars, fieldcodes):
value = ap.management.GetRasterProperties(p_in_raster, code)[0]
out[var] = literal_eval(value) if var != "wavelengths" else value
return out
## Extract raster values within polygons
class RasterInfo:
def __init__(self, raster):
self.info = get_raster_info(raster)
self.rows = np.array([self.info["top"] - n * self.info["y_size"] for n in range(1,
self.info["n_row"] + 1)]) if self.info["top"] > self.info["bottom"] else \
np.array([self.info["top"] + n * self.info["y_size"] for n in range(1,
self.info["n_row"] + 1)])
self.cols = np.array([self.info["left"] + n * self.info["x_size"] for n in range(1,
self.info["n_col"] + 1)])
self.sp_ref = ap.Describe(raster).spatialReference.exporttostring()
self.type = None
def select_cells_by_polygon(self, polygon):
### Get polygon information
ext = polygon.extent
xmin, xmax, ymin, ymax = ext.XMin, ext.XMax, ext.YMin, ext.YMax
tmp_pfeature = os.path.join("in_memory", "pft")
ap.management.CopyFeatures(polygon, tmp_pfeature)
tmp_poly = os.path.join("in_memory", "pgn")
ap.management.MakeFeatureLayer(tmp_pfeature, tmp_poly)
### Create points across polygon extent
xrange = np.where(np.logical_and(self.cols > xmin, self.cols < xmax))[0]
yrange = np.where(np.logical_and(self.rows > ymin, self.rows < ymax))[0]
indices = [[x, y] for x in xrange for y in yrange]
df = pd.DataFrame(data = indices, columns = ["c", "r"])
df["x"] = self.cols[df["c"]]
df["y"] = self.rows[df["r"]]
#df["OID"] = range(1, len(df) + 1)
with tmp.NamedTemporaryFile(delete = False) as tmp_table:
filename = tmp_table.name + ".txt"
df.to_csv(filename, sep = "\t", index = False)
tmp_table.close
tmp_plyr0 = os.path.join("in_memory", "plyr0")
tmp_0 = ap.MakeXYEventLayer_management(filename, "x", "y", tmp_plyr0, self.sp_ref)
# ArcGIS Pro: ap.management.XYTableToPoint(tmp_table, tmp_plyr, "x", "y")
### Copy to another layer (official ArcMap Website solution to missing OIDs)
tmp_plyr1 = os.path.join("in_memory", "plyr1")
ap.management.CopyFeatures(tmp_plyr0, tmp_plyr1)
tmp_plyrf = os.path.join("in_memory", "plyrf")
tmp_f = ap.management.MakeFeatureLayer(tmp_plyr1, tmp_plyrf)
### Select points within polygon
tmp_pts = os.path.join("in_memory", "pts")
selection = ap.management.SelectLayerByLocation(tmp_plyrf, overlap_type = "WITHIN",
select_features = tmp_poly, selection_type = "NEW_SELECTION")
ap.management.CopyFeatures(selection, tmp_pts)
out = ap.da.FeatureClassToNumPyArray(tmp_pts, ("c", "r"))
### Clear temporary files
ap.management.Delete(tmp_plyr0)
ap.management.Delete(tmp_plyr1)
ap.management.Delete(tmp_poly)
### Output
return out
With this code, I am able to extract column and row indices of the raster cells within a polygon.
However, this is not very efficient and it takes a lot of time to run.
