Extract Raster data per Polygon without extensions

Manu — Tue, 20 Dec 2022 20:19:39 GMT

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.

topic Extract Raster data per Polygon without extensions in Python Questions

Extract Raster data per Polygon without extensions