import arcpy
# Add polygons
arcpy.management.MakeFeatureLayer("C:/book/data/statisticalareas_demography2019.gdb/statisticalareas_demography2019", "stat")
# Add raster
arcpy.MakeRasterLayer_management("C:/book/output/carmel.tif", "carmel.tif")
# Subset
arcpy.management.SelectLayerByAttribute("stat", "NEW_SELECTION", '"SHEM_YISHUV" = \'חיפה\'')
arcpy.management.CopyFeatures("stat", "stat2")
# Extract
ZonalStatisticsAsTable("stat2", "STAT11", "carmel.tif", "C:/book/output/out.dbf", "DATA", "MEAN")
# Join table with polygons
arcpy.AddJoin_management("stat2", "STAT11", "out", "STAT11")
# Write
arcpy.conversion.FeatureClassToShapefile("stat2", "C:/book/output")
# Copy the layer to a new permanent feature class
arcpy.CopyFeatures_management(veg_joined_table, outFeature)