The final result

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: ...

Get ArcPy and ArcGIS – Geospatial Analysis with Python now with the O’Reilly learning platform.

O’Reilly members experience books, live events, courses curated by job role, and more from O’Reilly and nearly 200 top publishers.