datos_espaciales.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "datos_espaciales.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "private_outputs": true,
      "authorship_tag": "ABX9TyOHVGeEmibjOBZK/BQdRIPw",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "gpuClass": "standard"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/sporella/datos_geoespaciales_python/blob/main/datos_espaciales.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Instalar librerías necesarias\n",
        "\n",
        "En Google Colab ya se encuentran instaladas algunas librerías por defecto, pero para trabajar con datos geoespaciales debemos instalar algunas más.\n",
        "\n",
        "Podemos utilizar el instalador de paquetes `pip` al igual como se haría en la terminal de tu propio computador. \n",
        "\n",
        "> Para utilizar códigos que normalmente escribiríamos en la terminal, debemos anteponer `!`\n",
        "\n",
        "Instalaremos:\n",
        "\n",
        "- `pygeos` herramientas para análisis geoespaciales\n",
        "- `folium` para visualizaciones dinámicas (incluir `-U` para instalar la última versión)\n",
        "- `geopandas` para trabajar con datos vectoriales\n",
        "- `mapclassify` contiene clasificaciones para mapas de cloropletas"
      ],
      "metadata": {
        "id": "Qbwye-xHxOsW"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "! pip install pygeos\n",
        "! pip install folium -U\n",
        "! pip install geopandas\n",
        "! pip install mapclassify\n"
      ],
      "metadata": {
        "id": "MvOhbtZtoyqC"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Introducción\n",
        "## Primero debemos importar las librerías necesarias\n",
        "En python podemos importar las librerías con un alias, para geopandas se utiliza normalmente `gpd`.\n",
        "\n",
        "`geopandas`nos permite trabajar con datos vectoriales (puntos, líneas, polígonos, etc) en un formato tabular que nos permite realizar procesamiento de datos con la librería `pandas`."
      ],
      "metadata": {
        "id": "82C8Uuhd0IJH"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "zLnFClGcnPvr"
      },
      "outputs": [],
      "source": [
        "import geopandas as gpd\n",
        "import pandas as pd\n",
        "import folium"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Leer datos\n",
        "\n",
        "En este primer ejemplo trabajaremos con datos que trae geopandas. Son datos de Natural Earth con información de países (polígonos) y ciudades capitales (puntos).\n"
      ],
      "metadata": {
        "id": "U12QrE25P5oN"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "paises = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
        "paises.head()"
      ],
      "metadata": {
        "id": "PspYiJ-FoOHw"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "capitales = gpd.read_file(gpd.datasets.get_path(\"naturalearth_cities\"))\n",
        "capitales.head()"
      ],
      "metadata": {
        "id": "fQWR4_Wm8Gu9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Gráfico estático\n",
        "\n",
        "Podemos hacer un gráfico estático con el método plot() que utiliza la sintaxis de `matplotlib`"
      ],
      "metadata": {
        "id": "Vp3uFgdPQaRj"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "paises.plot();\n",
        "paises.plot(\"continent\");\n",
        "paises.plot(\"pop_est\");"
      ],
      "metadata": {
        "id": "BDrXn9YRPeRD"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "capitales.plot();"
      ],
      "metadata": {
        "id": "0oPqJFc5Pyp0"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Gráfico dinámico\n",
        "\n",
        "En sus versiones más recientes, geopandas nos permite usar la librería `folium`para hacer visualizaciones dinámicas mediante el método `.explore()`"
      ],
      "metadata": {
        "id": "aDBU96hvQ1Nd"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "paises.explore()"
      ],
      "metadata": {
        "id": "cpcW2S91Q0yZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "capitales.explore()"
      ],
      "metadata": {
        "id": "dAE5dgqSP-cA"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Proyección espacial\n",
        "\n",
        "Como vimos antes, es muy importante conocer la proyección en la que se encuentran nuestros datos, o incluso si es que tienen o no una proyección.\n",
        "\n",
        "Trabajar con datos en proyecciones diferentes nos lleva a incurrir en errores que no siempre son fáciles de detectar, por lo que siempre es necesario chequear las proyecciones.\n"
      ],
      "metadata": {
        "id": "-CMYwobsHe6y"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Proyección actual\n",
        "\n",
        "Podemos consultar la proyección de un set de datos utilizando `.crs`"
      ],
      "metadata": {
        "id": "RXbFeOSmUH1V"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "paises.crs"
      ],
      "metadata": {
        "id": "8C_iGvH0HmCR"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Cambiar proyección\n",
        "\n",
        "Para cambiar la proyección, utilizamos el método `to_crs()` e incluimos el número del epsg de destino.\n",
        "\n",
        "Podemos consultar diferentes epsg en https://epsg.io/"
      ],
      "metadata": {
        "id": "XnI9lT-hRydW"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "paises_psm = paises.to_crs(3857)\n",
        "paises_psm.crs"
      ],
      "metadata": {
        "id": "iTJ2RsevRvJO"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "paises.head()"
      ],
      "metadata": {
        "id": "HH-jD46XSLCl"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "paises_psm.head()"
      ],
      "metadata": {
        "id": "LzhsHiHESOFa"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "paises.plot(figsize = [10, 10]);"
      ],
      "metadata": {
        "id": "mF8hFpOoTiXN"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "paises_psm.plot(figsize = [10, 10]);"
      ],
      "metadata": {
        "id": "pFZ4y_4KSZz1"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "paises_psm.explore()"
      ],
      "metadata": {
        "id": "YsnFvR8QT3sB"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Existen proyecciones específicas para ciertos lugares del mundo, como por ejemplo, la Antártica"
      ],
      "metadata": {
        "id": "C4yqx-vGLUkn"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "antartica = paises[paises.continent == \"Antarctica\"]\n",
        "antartica.plot()"
      ],
      "metadata": {
        "id": "VK09kjwSSkpp"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "antartica_polar = antartica.to_crs(3031) # 5481\n",
        "antartica_polar.plot()\n"
      ],
      "metadata": {
        "id": "7ifwNNNHS9L5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "antartica_polar = antartica.to_crs(5481) # 5481\n",
        "antartica_polar.plot()\n"
      ],
      "metadata": {
        "id": "KlHLTXmfU6kU"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Trabajando con nuestros propios datos\n",
        "\n",
        "Por lo general, tendremos tres fuentes principales desde las cuales podemos obtener datos geospaciales:\n",
        "\n",
        "- **Archivos csv con coordenadas** de ubicaciones de algún tópico de interés\n",
        "\n",
        "- **Capas de servicios públicos o que estén en internet**. Son capas de polígonos, puntos o líneas que alguien más ha disponibilizado y que por lo general se encuentran en formato Shapefile (GIS), KML (GoogleEarth) o geojson.\n",
        "\n",
        "- **Capas hechas por nosotros mismos** es la forma menos común de obtener los datos, ya que debemos tener experiencia en el uso de Sistemas de Información Geográfica."
      ],
      "metadata": {
        "id": "8T_rMNfiUMIi"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Algunas fuentes de datos en Internet\n",
        "\n",
        "- [geonode CEDEUS](http://datos.cedeus.cl/)\n",
        "\n",
        "- [IDE Chile](https://www.ide.cl/)\n",
        "\n",
        "- [geodatos INE](https://www.ine.cl/herramientas/portal-de-mapas/geodatos-abiertos)\n",
        "\n",
        "- [Observatorio Ciudades UC](https://www.ine.cl/herramientas/portal-de-mapas/geodatos-abiertos)"
      ],
      "metadata": {
        "id": "xIlOSw0EUUEE"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Cargar archivos\n",
        "\n",
        "Geopandas nos permite cargar archivos geoespaciales de diferentes tipos, como .`shp.`, `.geojson`, `csv`, etc.\n"
      ],
      "metadata": {
        "id": "858JiiRpADHG"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Archivo shp\n",
        "\n",
        "Para cargar un archivo shp debemos poner la ruta del archivo shp., en el mismo directorio deben estar los otros archivos que componen el shapefile(.cpg, .dbf, .prj, .shx)"
      ],
      "metadata": {
        "id": "e2bRV5B__KqS"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "comunas_valparaiso = gpd.read_file(\"comunas_valparaiso.shp\")\n",
        "comunas_valparaiso.plot(figsize = [8,8])\n",
        "comunas_valparaiso.explore()"
      ],
      "metadata": {
        "id": "KkUFgtI8DlpG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Archivo Geojson\n",
        "\n",
        "Podemos leerlo desde el almacenamiento local o desde una url."
      ],
      "metadata": {
        "id": "ruD6IjI2_kgy"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "transporte = gpd.read_file(\"http://datos.cedeus.cl/geoserver/wfs?srsName=EPSG%3A4326&typename=geonode%3Atransporte_valpo&outputFormat=json&version=1.0.0&service=WFS&request=GetFeature\")\n",
        "#transporte[\"route_desc\"] = transporte.route_desc.replace(\"\", \"Sin Información\")\n",
        "transporte.explore()\n",
        "transporte.head()"
      ],
      "metadata": {
        "id": "S-1fLUhFAjyZ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "transporte.route_desc.unique()"
      ],
      "metadata": {
        "id": "ElkkrHlCm0Py"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "transporte[\"route_desc\"] = transporte.route_desc.replace(\"\", \"Sin información\")\n",
        "transporte"
      ],
      "metadata": {
        "id": "HvvEyKlMnb_9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "transporte.explore(\"route_desc\")"
      ],
      "metadata": {
        "id": "caaQiQ5Zm-wV"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Archivos csv\n",
        "\n",
        "Cargamos primero con pandas `pd.read_csv()` y luego convertimos con `gpd.GeoDataFrame()`"
      ],
      "metadata": {
        "id": "jeLe0iioG50E"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "monumentos = pd.read_csv(\"monumentos.csv\")\n",
        "monumentos.head()"
      ],
      "metadata": {
        "id": "mZ5UnsTOJcmk"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "type(monumentos)"
      ],
      "metadata": {
        "id": "alCr0tL3oM7Z"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "monumentos = gpd.GeoDataFrame(monumentos, geometry = gpd.points_from_xy(monumentos.lon, monumentos.lat))\n",
        "monumentos.explore()"
      ],
      "metadata": {
        "id": "8VdT0GrTA78v"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "En este caso es importante chequear la información de proyección geoespacial con `.crs` para darnos cuenta que debemos asignar el sistema de coordenadas al realizar la conversión:"
      ],
      "metadata": {
        "id": "e8Ikp2MKHSWh"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "print(monumentos.crs)"
      ],
      "metadata": {
        "id": "tUQ07BRTHfpe"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "monumentos = gpd.GeoDataFrame(monumentos, geometry = gpd.points_from_xy(monumentos.lon, monumentos.lat), crs = 4326) # inlcuir el argumento crs\n",
        "monumentos.explore()"
      ],
      "metadata": {
        "id": "-1Y3Pd15HkK9"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Cruce de información\n",
        "\n",
        "Nos gustaría generar un mapa solo para la región de Valparaíso pero nuestros puntos de monumentos nacionales no poseen información de región, comuna y provincia.\n",
        "\n",
        "En este caso, podemos utilizar un \"join espacial\". Para esto, lo primero que debemos hacer es confirmar que ambas capas tengan el mismo sistema de proyección:"
      ],
      "metadata": {
        "id": "sMSR25JKJ1BF"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "comunas_valparaiso.crs == monumentos.crs"
      ],
      "metadata": {
        "id": "-UBLgMlBKWBH"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "comunas_valparaiso.crs"
      ],
      "metadata": {
        "id": "YrKKjBgsqHGq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "monumentos.crs"
      ],
      "metadata": {
        "id": "ws01mLB7qKzD"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Entonces tenemos que convertir las capas al mismo sistema de coordenadas. Utilizaremos una proyección en metros (la de las comunas), así que transformaremos los monumentos a psm (3857) "
      ],
      "metadata": {
        "id": "-uxxXiQuKeYY"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "monumentos_psm = monumentos.to_crs(comunas_valparaiso.crs)\n",
        "comunas_valparaiso.crs == monumentos_psm.crs"
      ],
      "metadata": {
        "id": "u1-OR1FjKikR"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "monumentos.sjoin(comunas_valparaiso).shape"
      ],
      "metadata": {
        "id": "lpivVfwhqhzP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Y ahora podremos hacer el join"
      ],
      "metadata": {
        "id": "_Zam7N6LLG6l"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "monumentos_psm.shape\n",
        "comunas_valparaiso.shape"
      ],
      "metadata": {
        "id": "PO1jZBRSrFqb"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "monumentos_psm.sjoin(comunas_valparaiso).shape"
      ],
      "metadata": {
        "id": "S17trhy6q0Va"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "monumentos_valparaiso = monumentos_psm.sjoin(comunas_valparaiso[[\"Region\",\t\"Comuna\",\t\"Provincia\", \"geometry\"]])\n",
        "monumentos_valparaiso.explore()"
      ],
      "metadata": {
        "id": "k9iJfCiXLGKm"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "monumentos_valparaiso.explore(\"Comuna\", tiles = \"cartoDBPositron\")"
      ],
      "metadata": {
        "id": "7w7X-NtqMUSe"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "transporte.crs\n",
        "transporte_psm = transporte.to_crs(3857)"
      ],
      "metadata": {
        "id": "L_FOaB7isose"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "monumentos_valparaiso[\"distancia_km\"] = monumentos_valparaiso.geometry.apply(lambda x: transporte_psm.distance(x).min()/1000)\n",
        "monumentos_valparaiso.explore(\"distancia_km\", cmap = \"YlGnBu_r\")"
      ],
      "metadata": {
        "id": "BLNk7wXmNho0"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Mejorando nuestra visualización"
      ],
      "metadata": {
        "id": "J1fXGpK--Yf_"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "En esta sección trabajaremos de forma más libre, poniéndo énfasis en la importancia de consultar la documentación para obtener resultados más personalizados."
      ],
      "metadata": {
        "id": "cis8o0rz-poT"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "m = comunas_valparaiso.explore(tiles = \"CartoDB positron\", column = \"Provincia\", \n",
        "                               cmap = [\"#00e3d4\",\"#9f0f08\",\"#ac9bff\",\"#ecce42\",\"#4a308c\",\"#007735\",\"#ff97bf\"],\n",
        "                               tooltip = False,\n",
        "                               popup = [\"Comuna\", \"Provincia\"])\n",
        "monumentos_valparaiso.explore(m = m, marker_type = \"circle_marker\", marker_kwds = {\"color\" : \"#f5e9f7\"})\n",
        "transporte_psm.explore(m = m, column = \"route_desc\")\n",
        "folium.LayerControl().add_to(m)\n",
        "m"
      ],
      "metadata": {
        "id": "c_umvu9r-dJQ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "comunas_valparaiso.Provincia.nunique()"
      ],
      "metadata": {
        "id": "mdhviMlSugqh"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "m.save(\"mapa_valparaiso.html\")"
      ],
      "metadata": {
        "id": "P2EUvaUwzhph"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}