#1 first export from modelbuilder to calculate a heat risk index for …

…the region of Bonn

spatial-data-science/scripts/heat_risk_index.py

@@ -0,0 +1,160 @@
+import arcpy
+from arcpy.sa import *
+from sys import argv
+
+def initialize_arcpy():
+    """ Initialize arcpy settings and check necessary extensions. """
+    arcpy.env.overwriteOutput = True
+    arcpy.CheckOutExtension("3D")
+    arcpy.CheckOutExtension("spatial")
+    arcpy.CheckOutExtension("ImageExt")
+    arcpy.CheckOutExtension("ImageAnalyst")
+
+def generate_tessellation(output_feature_class, extent, size, spatial_ref):
+    """ Generate a tessellation grid. """
+    arcpy.management.GenerateTessellation(
+        Output_Feature_Class=output_feature_class,
+        Extent=extent,
+        Size=size,
+        Spatial_Reference=spatial_ref
+    )
+
+def spatial_join(target_features, join_features, out_feature_class):
+    """ Perform a spatial join between two feature sets. """
+    arcpy.analysis.SpatialJoin(
+        target_features=target_features,
+        join_features=join_features,
+        out_feature_class=out_feature_class,
+        match_option="LARGEST_OVERLAP"
+    )
+
+def zonal_statistics(in_zone_data, zone_field, in_value_raster, out_table, statistics_type="MAXIMUM"):
+    """ Calculate zonal statistics as a table. """
+    arcpy.sa.ZonalStatisticsAsTable(
+        in_zone_data,
+        zone_field,
+        in_value_raster,
+        out_table,
+        "DATA",
+        statistics_type
+    )
+
+def copy_raster(in_raster, out_rasterdataset, raster_format="TIFF"):
+    """ Copy a raster to a new dataset. """
+    arcpy.management.CopyRaster(in_raster, out_rasterdataset, format=raster_format)
+
+def reclassify_raster(in_raster, remap, out_raster):
+    """ Reclassify a raster based on value ranges. """
+    reclass_raster = arcpy.sa.Reclassify(in_raster, "Value", remap)
+    reclass_raster.save(out_raster)
+
+def join_field(in_data, in_field, join_table, join_field, fields):
+    """ Join fields from one table to another. """
+    arcpy.management.JoinField(
+        in_data=in_data,
+        in_field=in_field,
+        join_table=join_table,
+        join_field=join_field,
+        fields=fields
+    )
+
+def calculate_field(in_table, field, expression, field_type="FLOAT"):
+    """ Calculate a new field based on an expression. """
+    arcpy.management.CalculateField(
+        in_table=in_table,
+        field=field,
+        expression=expression,
+        field_type=field_type
+    )
+
+def standardize_field(in_table, fields, method="MIN-MAX", min_value=1, max_value=5):
+    """ Standardize fields using the given method. """
+    arcpy.management.StandardizeField(
+        in_table=in_table,
+        fields=fields,
+        method=method,
+        min_value=min_value,
+        max_value=max_value
+    )
+
+def delete_features(in_features):
+    """ Delete features from a feature class or layer. """
+    arcpy.management.DeleteFeatures(in_features=in_features)
+
+def hri_main(land_cover, zensus_2022):
+    """ Main function to execute the HRI analysis. """
+    initialize_arcpy()
+
+    # Generate tessellation
+    tessellation_output = "HRI_Hexagone"
+    generate_tessellation(
+        output_feature_class=tessellation_output, 
+        extent="781745.292120143 6556576.21979931 802689.19726414 6581479.0533047",
+        size="1500 SquareMeters",
+        spatial_ref="PROJCS[\"WGS_1984_Web_Mercator_Auxiliary_Sphere\",GEOGCS[\"GCS_WGS_1984\",...]"
+    )
+
+    # Spatial join
+    spatial_join_output = "HRI_Hexagone_SpatialJoin1"
+    spatial_join(tessellation_output, zensus_2022, spatial_join_output)
+
+    # Zonal statistics for surface temperature
+    surf_temp_max = "surf_temp_max"
+    zonal_statistics(spatial_join_output, "GRID_ID", "Multispectral Landsat.tif", surf_temp_max)
+
+    # Copy raster for tree canopy
+    copy_raster(land_cover, "tree_canopy.tif")
+
+    # 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")
+
+    # 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")
+
+    # Join fields for tree canopy statistics
+    join_field(spatial_join_output, "GRID_ID", tree_canopy_count, "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!")
+
+    # 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"])
+
+    # Copy raster for built-up area
+    copy_raster(land_cover, "built_up_area.tif")
+
+    # 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")
+
+    # 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")
+
+    # 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"])
+
+    # 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")
+
+    # Select layer by location
+    arcpy.management.SelectLayerByLocation(
+        in_layer=[spatial_join_output],
+        overlap_type="WITHIN",
+        select_features="Ortsteile_Bonn",
+        invert_spatial_relationship="INVERT"
+    )
+
+    # Delete features not in Bonn
+    delete_features(spatial_join_output)
+
+if __name__ == '__main__':
+    hri_main(*argv[1:])

spatial-data-science/scripts/unittest.py

@@ -0,0 +1,115 @@
+import unittest
+from hri_module import *
+import arcpy
+
+class TestHRIModule(unittest.TestCase):
+
+    def setUp(self):
+        # Setup code for test environment, mock data or objects can be initialized here
+        pass
+
+    def test_generate_tessellation(self):
+        output = "test_output_tessellation"
+        extent = "781745.292120143 6556576.21979931 802689.19726414 6581479.0533047"
+        size = "1500 SquareMeters"
+        spatial_ref = "WGS_1984"
+        try:
+            generate_tessellation(output, extent, size, spatial_ref)
+            success = True
+        except Exception as e:
+            success = False
+
+        self.assertTrue(success, "generate_tessellation failed.")
+
+    def test_spatial_join(self):
+        target_features = "target_features"
+        join_features = "join_features"
+        out_feature_class = "out_feature_class"
+        try:
+            spatial_join(target_features, join_features, out_feature_class)
+            success = True
+        except Exception as e:
+            success = False
+
+        self.assertTrue(success, "spatial_join failed.")
+
+    def test_zonal_statistics(self):
+        try:
+            zonal_statistics("zone_data", "GRID_ID", "value_raster", "out_table")
+            success = True
+        except Exception:
+            success = False
+
+        self.assertTrue(success, "zonal_statistics failed.")
+
+    def test_copy_raster(self):
+        try:
+            copy_raster("input_raster", "output_raster")
+            success = True
+        except Exception:
+            success = False
+
+        self.assertTrue(success, "copy_raster failed.")
+
+    def test_reclassify_raster(self):
+        try:
+            reclassify_raster("input_raster", "10 1;20 0;30 0", "output_raster")
+            success = True
+        except Exception:
+            success = False
+
+        self.assertTrue(success, "reclassify_raster failed.")
+
+    def test_join_field(self):
+        try:
+            join_field("in_data", "GRID_ID", "join_table", "join_field", ["field1", "field2"])
+            success = True
+        except Exception:
+            success = False
+
+        self.assertTrue(success, "join_field failed.")
+
+    def test_calculate_field(self):
+        try:
+            calculate_field("in_table", "field", "expression")
+            success = True
+        except Exception:
+            success = False
+
+        self.assertTrue(success, "calculate_field failed.")
+
+    def test_standardize_field(self):
+        try:
+            standardize_field("in_table", [["field1", "min_max1"], ["field2", "min_max2"]])
+            success = True
+        except Exception:
+            success = False
+
+        self.assertTrue(success, "standardize_field failed.")
+
+    def test_delete_features(self):
+        try:
+            delete_features("in_features")
+            success = True
+        except Exception:
+            success = False
+
+        self.assertTrue(success, "delete_features failed.")
+
+    def test_full_workflow(self):
+        """ Test the full workflow with all steps. This would need more elaborate mocking or test data. """
+        try:
+            # Running the full workflow would need valid test data
+            hri_main("Land_Cover_Living_Atlas", "Zensus_2022_Gitter_100m")
+            success = True
+        except Exception as e:
+            success = False
+
+        self.assertTrue(success, "hri_main workflow failed.")
+
+    def tearDown(self):
+        # Cleanup after tests, if necessary
+        pass
+
+if __name__ == '__main__':
+    unittest.main()