Mfcttrf

If we have a raster, say of integer elevation data for a country,
and one polygon shapefile, say of 300 river basins in that country, with a unique name for each, how would we most easily get output like this for them all?

basinID, gridcellelev
a, 320
a, 321
a, 320
b, 17
b, 18
b, 19

The most burdensome way seems to be filtering/converting the single shapefile into 300 shapefiles,
clipping the raster 300 times into 300 uniqueID rasters, reading them back in, generating individual tables for each basin, then combining them all together.

On the other hand the ideal way seems to be skipping the file generation, not saving the xy data, and creating the same table from just using one raster and one shapefile - by iteratively selecting the cells within a basin, stamping them with the basinID, creating a table, losing the coordinates, and keep repeating and appending that table until the 300th basin.

I'm not looking for any statistics, just the raw data listing of the grid cell elevations that would have been part of some standard clip method. I believe the attribute table that comes with the clipped raster from ArcMap is the counts/frequency of the cells. That output table format works for me too.

I don't know to minimally reproduce a raster and a polygon shapefile,
so I'd just be grateful for any tips/libraries/functions/examples. If starting with R, here's a starting point:

library(tidyverse)
library(raster)
library(rgdal)
library(sf)

elev_raster <- raster("spain_elev_meters.tif") #integer raster
basins <- readOGR("spainbasins.shp", "spainbasins") %>% st_as_sf() #unique basin ID column: `basinID` 

I'd prefer to do everything in R, but happy to try any batch methods in ArcGIS as well (I have 10.6 with Spatial Analyst, but not Pro).

While this is not a complete answer, as I do not have your exact files etc., it should serve as a close basis for what you are trying to do.

setwd("D:/blah/")
library(velox)
library(raster)
library(rgdal)
library(sf)
library(tidyverse)
elev_raster <- raster("ElevationRaster.tif")
basins <- readOGR("basin.shp", "basin") %>% st_as_sf() 
basins_ID <- data.frame(c(1:length(basins$geometry)),as.character(basins$name),stringsAsFactors = F) # create a frame to link geometry ID to the basin name
colnames(basins_ID) <- c("ID","Name") # a bit of house-keeping to make the frame easier to understand
vlx_data <- velox(elev_raster) # velox rasters are super fast for extracting. Normal raster-package ones are sloooow.
extracted_data <- vlx_data$extract(basins,df=T,small=T) # extracting with velox is easy
colnames(extracted_data) <- c("ID","Elev") # house keeping again
merged_data <- merge(extracted_data,basins_ID,by="ID") # and now we merge the data tables to make nice looking dataset
head(merged_data)

This results in this nice dataFrame:

 ID Elev Name
1 1 31.050 a
2 1 14.550 a
3 1 19.840 a
4 1 63.226 a
5 1 86.181 a
6 1 77.979 a

If you're comfortable with arcpy, you can try this script.

Basically, it loop through each basin, clip the input raster and write in the output table all the values found in the clipped raster along with the name of the current basin.

import arcpy
from arcpy.sa import *

arcpy.env.overwriteOutput = True

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# CHANGE THESE TO FIT YOUR NAMES!!!
# COMPLETE PATH AND NAME OF YOUR INPUT RASTER
in_raster = r"C:Umbertocorsiumbe_arcgisproGDEM30m_clip02.tif"
# COMPLETE PATH AND NAME OF YOUR BASIN SHAPE
basins = r"C:UsersminorauDocumentsArcGISDefault.gdbbasins"
# NAME OF THE FIELD YOU WANT TO PRESERVE FROM THE BASIN ATTRIBUTE TABLE
basin_field = "name"
# OUTPUT PATH TO STORE YOU FINAL TABLE (I USED A GDB)
output_gdb = r"C:UsersminorauDocumentsArcGISDefault.gdb"
out_table = arcpy.CreateTable_management(output_gdb, "output_table")

arcpy.AddField_management(out_table, basin_field, "TEXT")
arcpy.AddField_management(out_table, "Value", "LONG")

basins_lyr = arcpy.MakeFeatureLayer_management(basins, "basins_lyr")

# loop through each basin
with arcpy.da.SearchCursor(basins, [basin_field]) as b_cur:
 for basin in b_cur:
 arcpy.SelectLayerByAttribute_management(basins_lyr, "NEW_SELECTION", """0 = '1'""".format(basin_field, basin[0]))
 raster_clip = ExtractByMask(in_raster, basins_lyr)
 ##Build attribute table for single band raster dataset
 ##Overwrite the existing attribute table file
 #arcpy.BuildRasterAttributeTable_management(raster_clip, "Overwrite")
 temp_result = arcpy.TableToTable_conversion(raster_clip, output_gdb, "output")
 with arcpy.da.SearchCursor(temp_result, ['Value']) as r_cur:
 for val in r_cur:
 with arcpy.da.InsertCursor(out_table, [basin_field, "Value"]) as i_cur:
 i_cur.insertRow((basin[0], val[0]))

Here is a printscreen of the results and the inputs I've used to test it:

How about "Convert Raster to Polygon", then do a "Spatial Join" of those new elevation polygons with your basin polygons. Do it so the centroids of the elevation polygons being inside of a basin polygon as the match method.