From bd32308a71e4bde9e7f40218368ac8601631cc3e Mon Sep 17 00:00:00 2001 From: Jan Tschada Date: Fri, 27 Sep 2024 14:08:56 +0200 Subject: [PATCH] #1 Defined the Python module - named heatri - sample notebook --- .gitignore | 4 + .../Urban Heat Risk Index Sample.ipynb | 125 +++++++++++ spatial-data-science/scripts/README.md | 22 ++ .../scripts/build/lib/heatri/__init__.py | 20 ++ .../build/lib/heatri/heat_risk_index.py | 212 ++++++++++++++++++ .../scripts/heatri/__init__.py | 20 ++ .../scripts/heatri/heat_risk_index.py | 212 ++++++++++++++++++ spatial-data-science/scripts/pyproject.toml | 24 ++ 8 files changed, 639 insertions(+) create mode 100644 spatial-data-science/notebooks/Urban Heat Risk Index Sample.ipynb create mode 100644 spatial-data-science/scripts/README.md create mode 100644 spatial-data-science/scripts/build/lib/heatri/__init__.py create mode 100644 spatial-data-science/scripts/build/lib/heatri/heat_risk_index.py create mode 100644 spatial-data-science/scripts/heatri/__init__.py create mode 100644 spatial-data-science/scripts/heatri/heat_risk_index.py create mode 100644 spatial-data-science/scripts/pyproject.toml diff --git a/.gitignore b/.gitignore index 3ae0de1..3599364 100644 --- a/.gitignore +++ b/.gitignore @@ -32,6 +32,10 @@ bld/ [Oo]bj/ [Ll]og/ [Ll]ogs/ +*.egg-info/ + +# Jupyter notebooks checkpoints +.ipynb_checkpoints/ # Visual Studio 2015/2017 cache/options directory .vs/ diff --git a/spatial-data-science/notebooks/Urban Heat Risk Index Sample.ipynb b/spatial-data-science/notebooks/Urban Heat Risk Index Sample.ipynb new file mode 100644 index 0000000..a354b4e --- /dev/null +++ b/spatial-data-science/notebooks/Urban Heat Risk Index Sample.ipynb @@ -0,0 +1,125 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "id": "150b3605", + "metadata": {}, + "outputs": [], + "source": [ + "from arcgis.gis import GIS\n", + "from heatri import calculate_heat_risk_index" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "472f9ae8", + "metadata": {}, + "outputs": [], + "source": [ + "gis = GIS('home')" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "16bbe714", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "a0d5a5824fd6475b9255d5cd0189f3fe", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "MapView(layout=Layout(height='400px', width='100%'))" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "bonn_map = gis.map('Bonn, Germany')\n", + "bonn_map" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "16076048", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "790891.206928151" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bonn_map.extent['xmin']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ece266db", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/spatial-data-science/scripts/README.md b/spatial-data-science/scripts/README.md new file mode 100644 index 0000000..9c821c1 --- /dev/null +++ b/spatial-data-science/scripts/README.md @@ -0,0 +1,22 @@ +### Analysis + +Using ArcGIS, we apply spatial analysis techniques to identify patterns and correlations between urban features and temperature variations. This helps us pinpoint specific areas that are most susceptible to the UHI effect. + +The following screenshot shows the model to calculate the urban heat risk index based on landsat imagery, landcover and population data. + +![Screenshot ModelBuilder Heat Risk Index](https://raw.githubusercontent.com/EsriDE/urban-heat-risk-index/main/doc/img/HRI.svg) +*Screenshot: ModelBuilder Heat Risk Index* + +### Setup heat risk index + +Install the `heatri` module using pip: +``` +pip install --user . +``` + +```python + +from heatri import calculate_heat_risk_index + +# TODO +``` \ No newline at end of file diff --git a/spatial-data-science/scripts/build/lib/heatri/__init__.py b/spatial-data-science/scripts/build/lib/heatri/__init__.py new file mode 100644 index 0000000..0ff38a8 --- /dev/null +++ b/spatial-data-science/scripts/build/lib/heatri/__init__.py @@ -0,0 +1,20 @@ + +from .heat_risk_index import hri_main +from arcpy import Extent, SpatialReference + + +def calculate_heat_risk_index(x_min: float, y_min: float, x_max: float, y_max: float, wkid: int) -> str: + """ + Calculates the heat risk index for a defined geographic region. + + :return: The local path to the calculated heat risk index spatial bin features. + """ + + landsat_surf_temp = "https://landsat2.arcgis.com/arcgis/rest/services/Landsat/MS/ImageServer" + land_cover = "https://tiledimageservices.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services/European_Space_Agency_WorldCover_2021_Land_Cover_WGS84_7/ImageServer" + zensus_2022 = "https://services2.arcgis.com/jUpNdisbWqRpMo35/arcgis/rest/services/Zensus2022_grid_final/FeatureServer/0" + + spatial_reference = SpatialReference(wkid) + extent = Extent(XMin=x_min, YMin=y_min, XMax=x_max, YMax=y_max, spatial_reference=spatial_reference) + hri_spatial_bins = hri_main(landsat_surf_temp, land_cover, zensus_2022, extent, spatial_reference) + return hri_spatial_bins \ No newline at end of file diff --git a/spatial-data-science/scripts/build/lib/heatri/heat_risk_index.py b/spatial-data-science/scripts/build/lib/heatri/heat_risk_index.py new file mode 100644 index 0000000..4963784 --- /dev/null +++ b/spatial-data-science/scripts/build/lib/heatri/heat_risk_index.py @@ -0,0 +1,212 @@ +import arcpy +from arcpy.da import SearchCursor +from arcpy.sa import * +from arcpy import Extent +from arcpy import SpatialReference +import pandas as pd + + + +def initialize_arcpy(): + """ Initialize arcpy settings and check necessary extensions. """ + arcpy.env.overwriteOutput = True + arcpy.CheckOutExtension("3D") + arcpy.CheckOutExtension("spatial") + 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(landsat_surf_temp, land_cover, zensus_2022, extent, spatial_reference, workspace=None): + """ Main function to execute the HRI analysis. """ + initialize_arcpy() + + if workspace is None: + arcpy.env.workspace = arcpy.env.scratchGDB + else: + arcpy.env.workspace = workspace + + # Generate tessellation + tessellation_output = "HRI_Hexagone" + generate_tessellation( + output_feature_class=tessellation_output, + extent=extent, + size="1500 SquareMeters", + spatial_ref=spatial_reference + ) + + # Spatial join + spatial_join_output = "HRI_Spatial_Bins" + spatial_join(tessellation_output, zensus_2022, spatial_join_output) + surf_temp_layer = "Landsat Surface Temperature" + image_server_result = arcpy.management.MakeImageServerLayer(landsat_surf_temp, surf_temp_layer, processing_template="Band 10 Surface Temperature in Celsius", template=extent) + landsat_layer = image_server_result.getOutput(0) + definition_query = "((Best < 2000000) OR (Name LIKE 'Ov%')) AND (CloudCover < 0.05) AND (Month = 7 OR Month = 8)" + landsat_layer.definitionQuery = definition_query + + landsat_images = [] + with SearchCursor(landsat_layer, ["OBJECTID","AcquisitionDate"], where_clause=definition_query, spatial_filter=extent) as cursor: + for row in cursor: + landsat_images.append({"OBJECTID": row[0], "AcquisitionDate": row[1]}) + + landsat_images_df = pd.DataFrame(landsat_images) + landsat_images_sorted_df = landsat_images_df.sort_values(by="AcquisitionDate", ascending=False) + if landsat_images_sorted_df.empty: + raise ValueError("No valid Landsat Image found!") + + objectid = landsat_images_sorted_df["OBJECTID"][0] + landsat_layer.definitionQuery = f"OBJECTID={objectid}" + + landsat_surf_temp_output = "landsat_surf_temp" + with arcpy.EnvManager(extent=extent): + arcpy.management.CopyRaster( + in_raster=landsat_layer, + out_rasterdataset=landsat_surf_temp_output, + config_keyword="", + background_value=None, + nodata_value="", + onebit_to_eightbit="NONE", + colormap_to_RGB="NONE", + pixel_type="8_BIT_SIGNED", + scale_pixel_value="NONE", + RGB_to_Colormap="NONE", + format="TIFF", + transform="NONE", + process_as_multidimensional="CURRENT_SLICE", + build_multidimensional_transpose="NO_TRANSPOSE" + ) + + # 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) + + # Copy raster for tree canopy + tree_canopy_output = "tree_canopy" + with arcpy.EnvManager(extent=extent): + copy_raster(land_cover, tree_canopy_output) + + # Reclassify tree canopy raster + 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_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_table_output, "GRID_ID", ["SUM"]) + + # Calculate percentage tree cover and lacking + 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_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 + 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 + 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_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_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"]] + ) + 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") + gp_result = calculate_field(spatial_join_output, "HRI", field_expression, field_type="FLOAT") + return gp_result.getOutput(0) + + # 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) diff --git a/spatial-data-science/scripts/heatri/__init__.py b/spatial-data-science/scripts/heatri/__init__.py new file mode 100644 index 0000000..0ff38a8 --- /dev/null +++ b/spatial-data-science/scripts/heatri/__init__.py @@ -0,0 +1,20 @@ + +from .heat_risk_index import hri_main +from arcpy import Extent, SpatialReference + + +def calculate_heat_risk_index(x_min: float, y_min: float, x_max: float, y_max: float, wkid: int) -> str: + """ + Calculates the heat risk index for a defined geographic region. + + :return: The local path to the calculated heat risk index spatial bin features. + """ + + landsat_surf_temp = "https://landsat2.arcgis.com/arcgis/rest/services/Landsat/MS/ImageServer" + land_cover = "https://tiledimageservices.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services/European_Space_Agency_WorldCover_2021_Land_Cover_WGS84_7/ImageServer" + zensus_2022 = "https://services2.arcgis.com/jUpNdisbWqRpMo35/arcgis/rest/services/Zensus2022_grid_final/FeatureServer/0" + + spatial_reference = SpatialReference(wkid) + extent = Extent(XMin=x_min, YMin=y_min, XMax=x_max, YMax=y_max, spatial_reference=spatial_reference) + hri_spatial_bins = hri_main(landsat_surf_temp, land_cover, zensus_2022, extent, spatial_reference) + return hri_spatial_bins \ No newline at end of file diff --git a/spatial-data-science/scripts/heatri/heat_risk_index.py b/spatial-data-science/scripts/heatri/heat_risk_index.py new file mode 100644 index 0000000..4963784 --- /dev/null +++ b/spatial-data-science/scripts/heatri/heat_risk_index.py @@ -0,0 +1,212 @@ +import arcpy +from arcpy.da import SearchCursor +from arcpy.sa import * +from arcpy import Extent +from arcpy import SpatialReference +import pandas as pd + + + +def initialize_arcpy(): + """ Initialize arcpy settings and check necessary extensions. """ + arcpy.env.overwriteOutput = True + arcpy.CheckOutExtension("3D") + arcpy.CheckOutExtension("spatial") + 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(landsat_surf_temp, land_cover, zensus_2022, extent, spatial_reference, workspace=None): + """ Main function to execute the HRI analysis. """ + initialize_arcpy() + + if workspace is None: + arcpy.env.workspace = arcpy.env.scratchGDB + else: + arcpy.env.workspace = workspace + + # Generate tessellation + tessellation_output = "HRI_Hexagone" + generate_tessellation( + output_feature_class=tessellation_output, + extent=extent, + size="1500 SquareMeters", + spatial_ref=spatial_reference + ) + + # Spatial join + spatial_join_output = "HRI_Spatial_Bins" + spatial_join(tessellation_output, zensus_2022, spatial_join_output) + surf_temp_layer = "Landsat Surface Temperature" + image_server_result = arcpy.management.MakeImageServerLayer(landsat_surf_temp, surf_temp_layer, processing_template="Band 10 Surface Temperature in Celsius", template=extent) + landsat_layer = image_server_result.getOutput(0) + definition_query = "((Best < 2000000) OR (Name LIKE 'Ov%')) AND (CloudCover < 0.05) AND (Month = 7 OR Month = 8)" + landsat_layer.definitionQuery = definition_query + + landsat_images = [] + with SearchCursor(landsat_layer, ["OBJECTID","AcquisitionDate"], where_clause=definition_query, spatial_filter=extent) as cursor: + for row in cursor: + landsat_images.append({"OBJECTID": row[0], "AcquisitionDate": row[1]}) + + landsat_images_df = pd.DataFrame(landsat_images) + landsat_images_sorted_df = landsat_images_df.sort_values(by="AcquisitionDate", ascending=False) + if landsat_images_sorted_df.empty: + raise ValueError("No valid Landsat Image found!") + + objectid = landsat_images_sorted_df["OBJECTID"][0] + landsat_layer.definitionQuery = f"OBJECTID={objectid}" + + landsat_surf_temp_output = "landsat_surf_temp" + with arcpy.EnvManager(extent=extent): + arcpy.management.CopyRaster( + in_raster=landsat_layer, + out_rasterdataset=landsat_surf_temp_output, + config_keyword="", + background_value=None, + nodata_value="", + onebit_to_eightbit="NONE", + colormap_to_RGB="NONE", + pixel_type="8_BIT_SIGNED", + scale_pixel_value="NONE", + RGB_to_Colormap="NONE", + format="TIFF", + transform="NONE", + process_as_multidimensional="CURRENT_SLICE", + build_multidimensional_transpose="NO_TRANSPOSE" + ) + + # 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) + + # Copy raster for tree canopy + tree_canopy_output = "tree_canopy" + with arcpy.EnvManager(extent=extent): + copy_raster(land_cover, tree_canopy_output) + + # Reclassify tree canopy raster + 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_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_table_output, "GRID_ID", ["SUM"]) + + # Calculate percentage tree cover and lacking + 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_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 + 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 + 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_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_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"]] + ) + 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") + gp_result = calculate_field(spatial_join_output, "HRI", field_expression, field_type="FLOAT") + return gp_result.getOutput(0) + + # 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) diff --git a/spatial-data-science/scripts/pyproject.toml b/spatial-data-science/scripts/pyproject.toml new file mode 100644 index 0000000..e0f8fe5 --- /dev/null +++ b/spatial-data-science/scripts/pyproject.toml @@ -0,0 +1,24 @@ +[build-system] +requires = ['setuptools', 'setuptools-scm'] +build-backend = 'setuptools.build_meta' + +[tool.setuptools.packages.find] +include = ['heatri'] + +[project] +name = 'heatri' +version = '0.1a' +description = 'Calculates the heat risk index for a geographic region.' +readme = 'README.md' +requires-python = '>=3.11' +keywords = ['urban', 'heat'] +classifiers = [ + 'Programming Language :: Python :: 3', + 'License :: OSI Approved :: Apache Software License', + 'Operating System :: OS Independent', +] +dependencies = ['arcgis>=2.3'] + +[project.urls] +'Homepage' = 'https://www.esri.de' +'GitHub' = 'https://github.com/EsriDE/urban-heat-risk-index' \ No newline at end of file