From 6e1987a9b4d655ca64979885ae6806cf152af5b9 Mon Sep 17 00:00:00 2001 From: Jan Tschada Date: Tue, 24 Sep 2024 12:28:52 +0200 Subject: [PATCH] #1 Bunch of changes necessary - first of all we use strict the environment manager - nearly every copy raster must have a small extent not creating more than 4000 x 4000 pixels - the python script named unittest caused some sideeffects, cause numpy calculations tried to import this file and failed - calculate field had not a valid python expression, it was using Arcade --- .../scripts/heat_risk_index.py | 43 +++++++++++-------- .../{unittest.py => test_heat_risk_index.py} | 0 2 files changed, 25 insertions(+), 18 deletions(-) rename spatial-data-science/scripts/{unittest.py => test_heat_risk_index.py} (100%) diff --git a/spatial-data-science/scripts/heat_risk_index.py b/spatial-data-science/scripts/heat_risk_index.py index e2c2586..4a9d313 100644 --- a/spatial-data-science/scripts/heat_risk_index.py +++ b/spatial-data-science/scripts/heat_risk_index.py @@ -146,50 +146,57 @@ def hri_main(landsat_surf_temp, land_cover, zensus_2022, extent, spatial_referen # Zonal statistics for surface temperature as table surf_temp_table_output = "surf_temp_max_1" zonal_statistics(spatial_join_output, "GRID_ID", landsat_surf_temp_output, surf_temp_table_output) - return # Copy raster for tree canopy - copy_raster(land_cover, "tree_canopy.tif") + tree_canopy_output = "tree_canopy" + with arcpy.EnvManager(extent=extent): + copy_raster(land_cover, tree_canopy_output) # Reclassify tree canopy raster - reclassify_raster("tree_canopy.tif", "10 1;20 0;30 0;40 0;50 0;60 0;70 0;80 0;90 0;95 0;100 0", "reclass_tree_canopy.tif") + tree_canopy_reclassify_output = "reclass_tree_canopy" + reclassify_raster(tree_canopy_output, "10 1;20 0;30 0;40 0;50 0;60 0;70 0;80 0;90 0;95 0;100 0", tree_canopy_reclassify_output) # Zonal statistics for reclassified tree canopy - tree_canopy_count = "tree_canopy_count" - zonal_statistics(spatial_join_output, "GRID_ID", "reclass_tree_canopy.tif", tree_canopy_count, statistics_type="SUM") + tree_canopy_table_output = "tree_canopy_count" + zonal_statistics(spatial_join_output, "GRID_ID", tree_canopy_reclassify_output, tree_canopy_table_output, statistics_type="SUM") # Join fields for tree canopy statistics - join_field(spatial_join_output, "GRID_ID", tree_canopy_count, "GRID_ID", ["SUM"]) + join_field(spatial_join_output, "GRID_ID", tree_canopy_table_output, "GRID_ID", ["SUM"]) # Calculate percentage tree cover and lacking - calculate_field(tree_canopy_count, "PCT_Tree_Cover", "(!SUM! / !COUNT!) * 100") - calculate_field(tree_canopy_count, "PCT_Lacking", "100 - !PCT_Tree_Cover!") + calculate_field(tree_canopy_table_output, "PCT_Tree_Cover", "(!SUM! / !COUNT!) * 100") + calculate_field(tree_canopy_table_output, "PCT_Lacking", "100 - !PCT_Tree_Cover!") # Join the fields for surface temperature and tree canopy data - join_field(spatial_join_output, "GRID_ID", surf_temp_max, "GRID_ID", ["MAX"]) - join_field(spatial_join_output, "GRID_ID", tree_canopy_count, "GRID_ID", ["PCT_Tree_Cover", "PCT_Lacking"]) + join_field(spatial_join_output, "GRID_ID", surf_temp_table_output, "GRID_ID", ["MAX"]) + join_field(spatial_join_output, "GRID_ID", tree_canopy_table_output, "GRID_ID", ["PCT_Tree_Cover", "PCT_Lacking"]) # Copy raster for built-up area - copy_raster(land_cover, "built_up_area.tif") + built_up_area_output = "built_up_area" + with arcpy.EnvManager(extent=extent): + copy_raster(land_cover, built_up_area_output) # Reclassify built-up area raster - reclassify_raster("built_up_area.tif", "10 0;20 0;30 0;40 0;50 1;60 0;70 0;80 0;90 0;95 0;100 0", "reclass_built_up_area.tif") + built_up_area_reclassify_output = "reclass_built_up_area" + reclassify_raster(built_up_area_output, "10 0;20 0;30 0;40 0;50 1;60 0;70 0;80 0;90 0;95 0;100 0", built_up_area_reclassify_output) # Zonal statistics for reclassified built-up area - built_up_area_count = "built_up_area_count" - zonal_statistics(spatial_join_output, "GRID_ID", "reclass_built_up_area.tif", built_up_area_count, statistics_type="SUM") + built_up_area_table_output = "built_up_area_count" + zonal_statistics(spatial_join_output, "GRID_ID", built_up_area_reclassify_output, built_up_area_table_output, statistics_type="SUM") # Join the built-up area statistics - join_field(spatial_join_output, "GRID_ID", built_up_area_count, "GRID_ID", ["SUM"]) - calculate_field(built_up_area_count, "PCT_built_up_area", "(!SUM! / !COUNT!) * 100") - join_field(spatial_join_output, "GRID_ID", built_up_area_count, "GRID_ID", ["PCT_built_up_area"]) + join_field(spatial_join_output, "GRID_ID", built_up_area_table_output, "GRID_ID", ["SUM"]) + calculate_field(built_up_area_table_output, "PCT_built_up_area", "(!SUM! / !COUNT!) * 100") + join_field(spatial_join_output, "GRID_ID", built_up_area_table_output, "GRID_ID", ["PCT_built_up_area"]) # Final HRI calculation using standardized values standardize_field( spatial_join_output, [["Einwohner", "Einwohner_MIN_MAX"], ["MAX", "TEMP_MAX_MIN_MAX"], ["PCT_Lacking", "PCT_Lacking_MIN_MAX"]] ) - calculate_field(spatial_join_output, "HRI", "Sum($feature.TEMP_MAX_MIN_MAX, $feature.PCT_Lacking_MIN_MAX, $feature.Einwohner_MIN_MAX)", field_type="FLOAT") + field_expression = "!TEMP_MAX_MIN_MAX! + !PCT_Lacking_MIN_MAX! + !Einwohner_MIN_MAX!" + #calculate_field(spatial_join_output, "HRI", "Sum($feature.TEMP_MAX_MIN_MAX, $feature.PCT_Lacking_MIN_MAX, $feature.Einwohner_MIN_MAX)", field_type="FLOAT") + calculate_field(spatial_join_output, "HRI", field_expression, field_type="FLOAT") # Select layer by location # arcpy.management.SelectLayerByLocation( diff --git a/spatial-data-science/scripts/unittest.py b/spatial-data-science/scripts/test_heat_risk_index.py similarity index 100% rename from spatial-data-science/scripts/unittest.py rename to spatial-data-science/scripts/test_heat_risk_index.py