We produced a subset feature class of the bus stops that has the elevation values added as a field. This process could be repeated for the entire city, one neighborhood at a time, or it could be performed with the original elevation raster on the entire bus stops feature class to generate a value for each stop:
import arcpy arcpy.CheckOutExtension("Spatial") arcpy.env.overwriteOutput = True busStops = r'C:\Projects\PacktDB.gdb\SanFrancisco\Bus_Stops' sanFranciscoHoods = r'C:\Projects\SanFrancisco.gdb\SanFrancisco\SFFind_Neighborhoods' sfElevation = r'C:\Projects\SanFrancisco.gdb\sf_elevation' somaGeometry = [] sql = "name = 'South of Market'" with arcpy.da.SearchCursor(sanFranciscoHoods,['SHAPE@XY'],sql,None, True) as cursor: for row in cursor: somaGeometry.append(arcpy.Point(row[0][0],row[0][1])) somaElev = arcpy.sa.ExtractByPolygon(sfElevation, somaGeometry,"INSIDE") somaOutput = sfElevation.replace('sf_elevation','SOMA_elev') somaElev.save(somaOutput) print...