Original code was generated by ArcGIS ModelBuilder on : 2022-08-17 09:10:52 Modified by SM for daily US data starting with tmean 5/11/2023 This code loops through Jan 2008 and calculates tract averages and appends results together to a file geodatabase table Added step to delete daily results after appending 5/15/23 Modified by NMR to handle varying month lengths and worked on converting hard code to soft code 5/18/2023 Deleted old cells and created a clean version 2/15/2024 Added annotations for SAM 10/22/2024 Additional edits to annotations This code reads the weather data as .bil files as received from PRISM
%%skip##### value definitions ####### Set the workspace to the geodatabase where you want to create the table## Work is saved to local computer which is faster and avoids server connectivity problemsarcpy.env.workspace =r"C:\Users\Public\PublicArcGISProjects\Weather_WD02Public\PRISMDAILYCT10.gdb"#Create variable for geodatabasePRISMDAILYCT10_gdb =r"C:\Users\Public\PublicArcGISProjects\Weather_WD02Public\PRISMDAILYCT10.gdb"#Originally tract boundaries from server used then tract boundaries copied to local computer because of problems with server connections#CT10_ContUS = "W:\\ResearchAndData\\UHC_Data\\Census2010\\Geodatabases\\Census2010_US.gdb\\AlbersContUS\\CT_2010_ContUSland"#tract geometry was repaired prior to running and we made sure to use Albers Equal Area USGS projectionCT10_ContUS = PRISMDAILYCT10_gdb+"\CT_2010_ContUSlandrepairAlbers"#can loop through multiple measures but for now will process one measure at a time by setting range as (0,1) below#Leah Schinasi HEAT project does not need vpdmax and vpdminmet = ["ppt", "tmean", "tdmean", "tmax", "tmin", "vpdmax", "vpdmin"]start_year =2008#we ran one year at a time but could loop through multiple yearsend_year =2008import arcpy##### functions #####def TractZonalStatistics(): # To allow overwriting outputs change overwriteOutput option to True arcpy.env.overwriteOutput =True#array list of weather variables#for multiple variables change rangefor m inrange(0,1):#looping through selected yearsfor yr inrange(start_year, end_year +1):#output table for tract averages of daily measures## Set the name of the table you want to create CT10PRISMDAILYYY ="CT10PRISMDAILY"+met[m]+str(yr)#template table previously saved to workspace specifies tract Id date and measure as text, value as double TEMPLATE ="CT10PRISMDAILYTEMPLATE"## Create the table using the CreateTable_management() function arcpy.CreateTable_management(out_path=PRISMDAILYCT10_gdb, out_name=CT10PRISMDAILYYY, template=[TEMPLATE])#looping through selected months mth=1## full yearfor mth inrange(1,13): #To handle months with different numbers of days assume months have 31 days and use step67 .Exists to determine if data are available for each day day =1for day inrange(1,32):# Process: Make Feature Layer (Make Feature Layer) (management) CT10_ContUS_Land ="CT10_ContUS_Land"#A feature layer is created from the tract spatial dataset arcpy.management.MakeFeatureLayer(in_features=CT10_ContUS, out_layer=CT10_ContUS_Land, workspace="", field_info="OBJECTID OBJECTID VISIBLE NONE;Shape Shape VISIBLE NONE;STATEFP10 STATEFP10 VISIBLE NONE;COUNTYFP10 COUNTYFP10 VISIBLE NONE;TRACTCE10 TRACTCE10 VISIBLE NONE;GEOID10 GEOID10 VISIBLE NONE;NAME10 NAME10 VISIBLE NONE;NAMELSAD10 NAMELSAD10 VISIBLE NONE;MTFCC10 MTFCC10 VISIBLE NONE;FUNCSTAT10 FUNCSTAT10 VISIBLE NONE;ALAND10 ALAND10 VISIBLE NONE;AWATER10 AWATER10 VISIBLE NONE;INTPTLAT10 INTPTLAT10 VISIBLE NONE;INTPTLON10 INTPTLON10 VISIBLE NONE;tractarea_sqm tractarea_sqm VISIBLE NONE;Shape_Length Shape_Length VISIBLE NONE;Shape_Area Shape_Area VISIBLE NONE")#The PRISM 800 m .bil file is assigned to PRISMraster PRISMraster ="T:\\UHC\\Schinasi_HEAT\\PRISM_Data\\"+met[m]+"\\"+str(yr)+"\\prism_"+met[m]+"_us_30s_"+str(yr)+str(mth).zfill(2)+str(day).zfill(2)+".bil"if arcpy.Exists(PRISMraster): valuerst = arcpy.Raster(PRISMraster)#NOTE: same object name "tmean" has been used for all variables. Name "tmean" this will be used for other variables like ppt, tmin etc#Parallel processing factor set to 14 make sure computer has at least 14 processors# Process: Zonal Statistics as Table (Zonal Statistics as Table) (sa) CT10_ContUS_tmean = PRISMDAILYCT10_gdb +"\\CT10_ContUS_"+met[m]+"_"+str(yr) +str(mth).zfill(2)+str(day).zfill(2)#30 m cell size from metropolitan Philadelphia land cover data should be Ok for ContUS since only cell size not extent used#Using the number 30 does not work - maybe it assumes 30 degrees not 30 m?#I got errors using arcpy.sa that Saima used (spatial analysis module)#When I tried with Model Builder arcpy.ia (image analysis module) was used so I switched to iawith arcpy.EnvManager(cellSize="W:\\ResearchAndData\\UHC_Data\\NLCD\\Geodatabases\\Land_Cover_2016\\NLCD16.gdb\\NLCD16_PRC_Clip", parallelProcessingFactor="14"): arcpy.ia.ZonalStatisticsAsTable(in_zone_data=CT10_ContUS_Land, zone_field="GEOID10", in_value_raster=valuerst, out_table=CT10_ContUS_tmean, ignore_nodata="DATA", statistics_type="MEAN", process_as_multidimensional="CURRENT_SLICE", percentile_values=[90], percentile_interpolation_type="AUTO_DETECT")#Reformat zonal statistics output append to C:\Users\Public keep tables in geodatabase format until end# I created date as text field for now the convert time field can be used to convert dates as text to date format# Process: Make Table View (Make Table View) (management) CT10_PRC_tmean_View ="CT10_ContUS_tmean_View" arcpy.management.MakeTableView(in_table=CT10_ContUS_tmean, out_view=CT10_PRC_tmean_View, where_clause="", workspace="", field_info="OBJECTID OBJECTID VISIBLE NONE;GEOID10 GEOID10 VISIBLE NONE;ZONE_CODE ZONE_CODE HIDDEN NONE;COUNT COUNT HIDDEN NONE;AREA AREA HIDDEN NONE;MEAN MEAN VISIBLE NONE")# Process: Add Field (Add Field) (management) CT10_PRC_tmean_2 = arcpy.management.AddField(in_table=CT10_ContUS_tmean, field_name="DATETXT", field_type="TEXT", field_precision=None, field_scale=None, field_length=10, field_alias="", field_is_nullable="NULLABLE", field_is_required="NON_REQUIRED", field_domain="")[0] CT10_PRC_tmean_3 = arcpy.management.AddField(in_table=CT10_PRC_tmean_2, field_name="MEASURE", field_type="TEXT", field_precision=None, field_scale=None, field_length=6, field_alias="", field_is_nullable="NULLABLE", field_is_required="NON_REQUIRED", field_domain="")[0]# Process: Calculate Field (Calculate Field) (management) CT10_PRC_tmean_4 = arcpy.management.CalculateField(in_table=CT10_PRC_tmean_3, field="DATETXT", expression='"'+str(mth).zfill(2) +"/"+str(day).zfill(2) +"/"+str(yr)+'"', expression_type="PYTHON3", code_block="", field_type="TEXT", enforce_domains="NO_ENFORCE_DOMAINS")[0] CT10_PRC_tmean_5 = arcpy.management.CalculateField(in_table=CT10_PRC_tmean_4, field="MEASURE", expression='"'+met[m]+'"', expression_type="PYTHON3", code_block="", field_type="TEXT", enforce_domains="NO_ENFORCE_DOMAINS")[0]# Process: Append (Append) (management) CT10PRISMDAILYYY_2 = arcpy.management.Append(inputs=[CT10_PRC_tmean_5], target=CT10PRISMDAILYYY, schema_type="NO_TEST", field_mapping='GEOID10 \"GEOID10\" true true false 11 Text 0 0,First,#,'+CT10_PRC_tmean_5+',GEOID10,0,11;DATETXT \"DATETXT\" true true false 10 Text 0 0,First,#,'+CT10_PRC_tmean_5+',DATETXT,0,10;MEASURE \"MEASURE\" true true false 6 Text 0 0,First,#,'+CT10_PRC_tmean_5+',MEASURE,0,6;VALUE \"VALUE\" true true false 8 Double 0 0,First,#,'+CT10_PRC_tmean_5+',MEAN,-1,-1', subtype ="", expression="")[0]#delete unneeded intermediate files arcpy.management.Delete(CT10_PRC_tmean_5)print("Day"+str(day) +"done") day = day +1 mth = mth +1 yr = yr +1print(met[m] +"done")##### Main Code ##### if__name__=='__main__':# Global Environment settingswith arcpy.EnvManager(scratchWorkspace=r"C:\Users\Public\PublicArcGISProjects\Weather_WD02Public\PRISMDAILYCT10.gdb", workspace=r"C:\Users\Public\PublicArcGISProjects\Weather_WD02Public\PRISMDAILYCT10.gdb"): TractZonalStatistics()print("Zonal Statistics Done")