|
| 1 | +{ |
| 2 | + "cells": [ |
| 3 | + { |
| 4 | + "cell_type": "markdown", |
| 5 | + "id": "b4c12bc8", |
| 6 | + "metadata": {}, |
| 7 | + "source": [ |
| 8 | + "## Group #3\n" |
| 9 | + ] |
| 10 | + }, |
| 11 | + { |
| 12 | + "cell_type": "markdown", |
| 13 | + "id": "e2750e5d", |
| 14 | + "metadata": {}, |
| 15 | + "source": [ |
| 16 | + "1. Get the 15 variables from this raster for all Peru departments polygons. This is the link where shapefiles are located. This is the link of the source raster. The values should be the percentage of district area cover by this specific Morphological Settlement Zone." |
| 17 | + ] |
| 18 | + }, |
| 19 | + { |
| 20 | + "cell_type": "code", |
| 21 | + "execution_count": null, |
| 22 | + "id": "5c9d4aa8", |
| 23 | + "metadata": {}, |
| 24 | + "outputs": [], |
| 25 | + "source": [ |
| 26 | + "# pip install rasterio\n", |
| 27 | + "# %pip install geopandas matplotlib shapely rasterio numpy pandas sklearn-xarray -q\n", |
| 28 | + "# !pip install ipywidgets\n", |
| 29 | + "# !jupyter labextension install @jupyter-widgets/jupyterlab-manager\n", |
| 30 | + "# pip install rasterstats\n", |
| 31 | + "# %pip install geopandas matplotlib shapely rasterio numpy pandas sklearn-xarray -q\n", |
| 32 | + "# %pip install git+https://github.com/jgrss/geowombat -q" |
| 33 | + ] |
| 34 | + }, |
| 35 | + { |
| 36 | + "cell_type": "code", |
| 37 | + "execution_count": null, |
| 38 | + "id": "4789905c", |
| 39 | + "metadata": {}, |
| 40 | + "outputs": [], |
| 41 | + "source": [ |
| 42 | + "\n", |
| 43 | + "import geopandas as gpd\n", |
| 44 | + "import rasterio\n", |
| 45 | + "from rasterio.merge import merge\n", |
| 46 | + "from rasterio.plot import show\n", |
| 47 | + "from shapely.geometry import mapping\n", |
| 48 | + "import rasterio\n", |
| 49 | + "from pathlib import Path\n", |
| 50 | + "import matplotlib.pyplot as plt\n", |
| 51 | + "import matplotlib.patheffects as pe\n", |
| 52 | + "import os\n", |
| 53 | + "from rasterio.mask import mask\n", |
| 54 | + "from rasterstats import zonal_stats\n", |
| 55 | + "import pandas as pd\n", |
| 56 | + "#import geowombat as gw ## Chequear gw en collab\n", |
| 57 | + "\n" |
| 58 | + ] |
| 59 | + }, |
| 60 | + { |
| 61 | + "cell_type": "code", |
| 62 | + "execution_count": null, |
| 63 | + "id": "a3813d15", |
| 64 | + "metadata": {}, |
| 65 | + "outputs": [], |
| 66 | + "source": [ |
| 67 | + "import geopandas as gpd\n", |
| 68 | + "\n", |
| 69 | + "# Ruta al shapefile\n", |
| 70 | + "shapefile_path = r'../../_data/INEI_LIMITE_DEPARTAMENTAL'\n", |
| 71 | + "\n", |
| 72 | + "# Leer el shapefile usando Geopandas\n", |
| 73 | + "departmentos = gpd.read_file(shapefile_path)\n", |
| 74 | + "departmentos = departments.to_crs('esri:54009')\n", |
| 75 | + "\n", |
| 76 | + "print(departmentos.head())" |
| 77 | + ] |
| 78 | + }, |
| 79 | + { |
| 80 | + "cell_type": "code", |
| 81 | + "execution_count": null, |
| 82 | + "id": "3c91017c", |
| 83 | + "metadata": {}, |
| 84 | + "outputs": [], |
| 85 | + "source": [ |
| 86 | + "from google.colab import drive\n", |
| 87 | + "drive.mount('/content/drive')" |
| 88 | + ] |
| 89 | + }, |
| 90 | + { |
| 91 | + "cell_type": "code", |
| 92 | + "execution_count": null, |
| 93 | + "id": "3ac57118", |
| 94 | + "metadata": {}, |
| 95 | + "outputs": [], |
| 96 | + "source": [ |
| 97 | + "# Rutas en collab\n", |
| 98 | + "\n", |
| 99 | + "raster = [\"drive/content/diplomado/assignment_10/tif_files/GHS_BUILT_C_MSZ_E2018_GLOBE_R2023A_54009_10_V1_0_R10_C11.tif\",\n", |
| 100 | + " \"drive/content/diplomado/assignment_10/tif_files/GHS_BUILT_C_MSZ_E2018_GLOBE_R2023A_54009_10_V1_0_R10_C12.tif\",\n", |
| 101 | + " \"drive/content/diplomado/assignment_10/tif_files/GHS_BUILT_C_MSZ_E2018_GLOBE_R2023A_54009_10_V1_0_R11_C11.tif\",\n", |
| 102 | + " \"drive/content/diplomado/assignment_10/tif_files/GHS_BUILT_C_MSZ_E2018_GLOBE_R2023A_54009_10_V1_0_R11_C12.tif\",\n", |
| 103 | + " \"drive/content/diplomado/assignment_10/tif_files/GHS_BUILT_C_MSZ_E2018_GLOBE_R2023A_54009_10_V1_0_R12_C11.tif\",\n", |
| 104 | + " \"drive/content/diplomado/assignment_10/tif_files/GHS_BUILT_C_MSZ_E2018_GLOBE_R2023A_54009_10_V1_0_R12_C12.tif\"\n", |
| 105 | + " ]\n", |
| 106 | + "\n", |
| 107 | + "statistics = []" |
| 108 | + ] |
| 109 | + }, |
| 110 | + { |
| 111 | + "cell_type": "code", |
| 112 | + "execution_count": null, |
| 113 | + "id": "7786c4c9", |
| 114 | + "metadata": {}, |
| 115 | + "outputs": [], |
| 116 | + "source": [ |
| 117 | + "# Iterar sobre cada archivo raster\n", |
| 118 | + "for raster_path in raster:\n", |
| 119 | + " # estadísticas zonales\n", |
| 120 | + " stats = zonal_stats(departmentos, raster_path, stats=\"count sum\", categorical=True, all_touched=True)\n", |
| 121 | + " # transforma DF\n", |
| 122 | + " stats_df = pd.DataFrame(stats)\n", |
| 123 | + " # Concatenar \n", |
| 124 | + " df = pd.concat([departmentos.reset_index(drop=True), stats_df], axis=1)\n", |
| 125 | + " statistics.append(df)" |
| 126 | + ] |
| 127 | + }, |
| 128 | + { |
| 129 | + "cell_type": "code", |
| 130 | + "execution_count": null, |
| 131 | + "id": "f5476429", |
| 132 | + "metadata": {}, |
| 133 | + "outputs": [], |
| 134 | + "source": [ |
| 135 | + "# Concatenar todos los DataFrames en all_stats en un único DataFrame\n", |
| 136 | + "final_df = pd.concat(statistics, ignore_index=True)\n", |
| 137 | + "\n", |
| 138 | + "print(final_df.head())\n", |
| 139 | + "\n", |
| 140 | + "#Reemplazar todos los NaN por 0\n", |
| 141 | + "final_df.fillna(0, inplace=True)" |
| 142 | + ] |
| 143 | + }, |
| 144 | + { |
| 145 | + "cell_type": "code", |
| 146 | + "execution_count": null, |
| 147 | + "id": "14984616", |
| 148 | + "metadata": {}, |
| 149 | + "outputs": [], |
| 150 | + "source": [] |
| 151 | + }, |
| 152 | + { |
| 153 | + "cell_type": "markdown", |
| 154 | + "id": "9724792d", |
| 155 | + "metadata": {}, |
| 156 | + "source": [ |
| 157 | + "2. Then you are going to generate choropleth map using folium for these 15 variables." |
| 158 | + ] |
| 159 | + }, |
| 160 | + { |
| 161 | + "cell_type": "code", |
| 162 | + "execution_count": null, |
| 163 | + "id": "a15a1cd1", |
| 164 | + "metadata": {}, |
| 165 | + "outputs": [], |
| 166 | + "source": [ |
| 167 | + "# Definir la área de pixeles (resolución del raster al cuadrado)\n", |
| 168 | + "spatial_resolution = 10 # meters\n", |
| 169 | + "pixel_area = spatial_resolution ** 2\n", |
| 170 | + "\n", |
| 171 | + "# Iterar por cada una de las 15 categorías MSZ\n", |
| 172 | + "for category in range(1, 16):\n", |
| 173 | + " # Crear una nueva columna para cada categoría MSZ\n", |
| 174 | + " column_name = f'MSZ_{category}_coverage'\n", |
| 175 | + " \n", |
| 176 | + " # Usar una función lambda con .apply para calcular el porcentaje de cobertura\n", |
| 177 | + " coverage_function = lambda row: ((row.get(str(category), 0) * pixel_area) / row['geometry'].area) * 100 if str(category) in row else 0\n", |
| 178 | + " \n", |
| 179 | + " # Aplicar función al DataFrame\n", |
| 180 | + " final_df[column_name] = final_df.apply(coverage_function, axis=1)" |
| 181 | + ] |
| 182 | + }, |
| 183 | + { |
| 184 | + "cell_type": "markdown", |
| 185 | + "id": "f2fbdab7", |
| 186 | + "metadata": {}, |
| 187 | + "source": [ |
| 188 | + "3. Save your html in the same folder of your JN. Name your HTML as your branch. This HTML should have all these layers. Please do not forget to use Layer Control.\n" |
| 189 | + ] |
| 190 | + }, |
| 191 | + { |
| 192 | + "cell_type": "code", |
| 193 | + "execution_count": null, |
| 194 | + "id": "c4c88d45", |
| 195 | + "metadata": {}, |
| 196 | + "outputs": [], |
| 197 | + "source": [ |
| 198 | + "# Crear un mapa centrado en una ubicación específica (puedes ajustar esto según tus datos)\n", |
| 199 | + "map_center = [-12.0757538, -76.9863174]\n", |
| 200 | + "mymap = folium.Map(location=map_center, zoom_start=5)\n", |
| 201 | + "\n", |
| 202 | + "# Agregar las capas al mapa\n", |
| 203 | + "for category in range(1, 16):\n", |
| 204 | + " layer_name = f'MSZ_{category}_coverage'\n", |
| 205 | + " geojson_data = final_df.to_json()\n", |
| 206 | + " \n", |
| 207 | + " folium.Choropleth(\n", |
| 208 | + " geo_data=geojson_data,\n", |
| 209 | + " data=final_df,\n", |
| 210 | + " columns=['geometry', layer_name],\n", |
| 211 | + " key_on='feature.properties.geometry',\n", |
| 212 | + " fill_color='YlOrRd', \n", |
| 213 | + " fill_opacity=0.7,\n", |
| 214 | + " line_opacity=0.2,\n", |
| 215 | + " legend_name=f'MSZ {category} Coverage'\n", |
| 216 | + " ).add_to(mymap)\n", |
| 217 | + "\n", |
| 218 | + "# Agregar el control de capas al mapa\n", |
| 219 | + "folium.LayerControl().add_to(mymap)\n", |
| 220 | + "\n", |
| 221 | + "# Guardar el mapa como un archivo HTML\n", |
| 222 | + "g3a10DFG_html = 'g3a10DFG.html' \n", |
| 223 | + "mymap.save(g3a10DFG_html)" |
| 224 | + ] |
| 225 | + } |
| 226 | + ], |
| 227 | + "metadata": { |
| 228 | + "kernelspec": { |
| 229 | + "display_name": "Python 3 (ipykernel)", |
| 230 | + "language": "python", |
| 231 | + "name": "python3" |
| 232 | + }, |
| 233 | + "language_info": { |
| 234 | + "codemirror_mode": { |
| 235 | + "name": "ipython", |
| 236 | + "version": 3 |
| 237 | + }, |
| 238 | + "file_extension": ".py", |
| 239 | + "mimetype": "text/x-python", |
| 240 | + "name": "python", |
| 241 | + "nbconvert_exporter": "python", |
| 242 | + "pygments_lexer": "ipython3", |
| 243 | + "version": "3.11.5" |
| 244 | + } |
| 245 | + }, |
| 246 | + "nbformat": 4, |
| 247 | + "nbformat_minor": 5 |
| 248 | +} |
0 commit comments