{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Data layer" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using angles for meter grid.\n", "Using angles for meter grid.\n", "Using angles for meter grid.\n", "Using angles for meter grid.\n", "Using angles for meter grid.\n" ] } ], "source": [ "import json\n", "import io\n", "import os\n", "import re\n", "import itertools as iter\n", "import numpy as np\n", "\n", "# data\n", "import pandas as pd\n", "import geopandas as gpd\n", "import shapely\n", "import duckdb\n", "import h3\n", "\n", "import scalenav.data as snd\n", "import scalenav.scale_nav as sn\n", "from scalenav.plotting import cmap\n", "# from scipy import io as io\n", "# import nctoolkit as nc\n", "# import xarray as xr\n", "# import rioxarray as rx\n", "import glob\n", "import ibis as ib\n", "from ibis import _\n", "import ibis.selectors as s\n", "ib.options.interactive = True\n", "\n", "# plots\n", "\n", "# from datashader import transfer_functions as tf, reductions as rd\n", "import pypalettes as pypal\n", "import pydeck as pdk\n", "from seaborn import color_palette" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "conn = ib.connect(\"duckdb://\")\n", "conn.raw_sql(\"install spatial; load spatial\")" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'/Users/cenv1069/Documents/scale_nav/source/notebook'" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%pwd" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "processed_data_folder = \"../../outputs/\"" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.7238241365054198" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# trying an interpolation method for the right resolution from a given grid size\n", "dist = 0\n", "\n", "alpha = 5/np.log(1_000)\n", "\n", "A = 13 + alpha*np.log(10)\n", "\n", "A \n", "alpha\n", "# 13+alpha*np.log(10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setup" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'10': {'h3_res': 14, 'nn': []},\n", " '100': {'h3_res': 12, 'nn': []},\n", " '1000': {'h3_res': 11, 'nn': []},\n", " '5000': {'h3_res': 10, 'nn': []},\n", " '10000': {'h3_res': 8, 'nn': []}}" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "extra_res = 0\n", "grid_params = []\n", "res_params = []\n", "\n", "if dist > 0 : \n", " # extra_res=round(A-np.log(alpha*dist))\n", "\n", " grid_params = [10,100,1000,5000,10000,dist]\n", " grid_params.sort()\n", " \n", " # res_params = [13,12,11,10,8,extra_res]\n", " res_params = [round(A-alpha*np.log(size)) for size in grid_params]\n", " # res_params.sort(reverse=True)\n", "\n", "else:\n", " grid_params = [10,100,1000,5000,10000]\n", " res_params = [14,12,11,10,8]\n", "\n", "res_params\n", "\n", "rast_to_h3 = { str(size): {\"h3_res\" : res, \"nn\" : [] } for (size,res) in zip(grid_params,res_params)}\n", "rast_to_h3" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using angles for meter grid.\n", "Using angles for meter grid.\n", "Using angles for meter grid.\n", "Using angles for meter grid.\n", "Using angles for meter grid.\n", "Using base resolution: 14\n" ] } ], "source": [ "grid_param = 10\n", "\n", "rast_to_h3 = snd.rast_to_h3_map(x=8,y=12,ref=\"m\",dist=dist)\n", "\n", "base_res = rast_to_h3[str(grid_param)][\"h3_res\"]\n", "neighbs = rast_to_h3[str(grid_param)][\"nn\"]\n", "\n", "print(\"Using base resolution: \",base_res)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(0, 0),\n", " (2, 1),\n", " (1, -1),\n", " (1, 1),\n", " (0, 2),\n", " (-1, 0),\n", " (0, -2),\n", " (-2, -1),\n", " (1, 2),\n", " (2, 0),\n", " (0, 1),\n", " (-1, -1),\n", " (1, -2),\n", " (3, 1),\n", " (1, 0),\n", " (2, 2),\n", " (-1, -2),\n", " (0, -1)]" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "neighbs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### GHSL MSZ data set\n", "\n", "could be moved to another notebook." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# close up london-oxford\n", "# limits = [-1.486,51.138,0.352,51.882]\n", "\n", "# wide europe window\n", "# limits = [-14.06,48.34,31.95,59.27]\n", "\n", "\n", "# Kano Nigeria\n", "limits = [8.46085,11.96741,8.53552,12.01733]\n", "\n", "# Lagos\n", "# limits = [3.0624,6.3699,3.5925,6.6592]" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "# src_file = \"outputs/kummu.parquet\"\n", "src_file = \"new_process.parquet\"\n", "\n", "layer = conn.read_parquet(f\"{processed_data_folder}{src_file}\",table_name=\"layer\")" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "msz_categories = {\n", " \"1\" : \"Open spaces, low vegetation\",\n", " \"2\" : \"Open spaces, medium vegetation\",\n", " \"3\" : \"Open spaces, high vegetation\",\n", " \"4\" : \"Open spaces, water surfaces\",\n", " \"5\" : \"Open spaces, road surfaces\",\n", " \"11\" : \"Built spaces, residential building height <= 3m\",\n", " \"12\" : \"Built spaces, residential building 3m < height <= 6m\",\n", " \"13\" : \"Built spaces, residential building 6m < height <= 15m\",\n", " \"14\" : \"Built spaces, residential building 15m < height <= 30m\",\n", " \"15\" : \"Built spaces, residential building 30m < height\",\n", " \"21\" : \"Built spaces, non-residential building height <= 3m\",\n", " \"22\" : \"Built spaces, non-residential building 3m < height <= 6m\",\n", " \"23\" : \"Built spaces, non-residential building 6m < height <= 15m\",\n", " \"24\" : \"Built spaces, non-residential building 15m < height <= 30m\",\n", " \"25\" : \"Built spaces, non-residential building 30m < height\",\n", "}" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['layer']\n" ] }, { "data": { "text/html": [ "
┏━━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┓\n",
       "┃ lon       ┃ lat       ┃ band_var ┃\n",
       "┡━━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━┩\n",
       "│ float32   │ float32   │ float32  │\n",
       "├───────────┼───────────┼──────────┤\n",
       "│ -0.322563 │ 16.243849 │      1.0 │\n",
       "│ -0.322461 │ 16.243849 │      1.0 │\n",
       "│ -0.322359 │ 16.243849 │      1.0 │\n",
       "│ -0.322563 │ 16.243767 │      1.0 │\n",
       "│ -0.322461 │ 16.243767 │     11.0 │\n",
       "│ -0.322359 │ 16.243767 │     11.0 │\n",
       "│ -0.322563 │ 16.243685 │      1.0 │\n",
       "│ -0.322461 │ 16.243685 │     11.0 │\n",
       "│ -0.322358 │ 16.243685 │     11.0 │\n",
       "│ -0.322153 │ 16.242781 │      2.0 │\n",
       "└───────────┴───────────┴──────────┘\n",
       "
\n" ], "text/plain": [ "┏━━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━┓\n", "┃\u001b[1m \u001b[0m\u001b[1mlon\u001b[0m\u001b[1m \u001b[0m\u001b[1m \u001b[0m┃\u001b[1m \u001b[0m\u001b[1mlat\u001b[0m\u001b[1m \u001b[0m\u001b[1m \u001b[0m┃\u001b[1m \u001b[0m\u001b[1mband_var\u001b[0m\u001b[1m \u001b[0m┃\n", "┡━━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━┩\n", "│ \u001b[2mfloat32\u001b[0m │ \u001b[2mfloat32\u001b[0m │ \u001b[2mfloat32\u001b[0m │\n", "├───────────┼───────────┼──────────┤\n", "│ \u001b[1;36m-0.322563\u001b[0m │ \u001b[1;36m16.243849\u001b[0m │ \u001b[1;36m1.0\u001b[0m │\n", "│ \u001b[1;36m-0.322461\u001b[0m │ \u001b[1;36m16.243849\u001b[0m │ \u001b[1;36m1.0\u001b[0m │\n", "│ \u001b[1;36m-0.322359\u001b[0m │ \u001b[1;36m16.243849\u001b[0m │ \u001b[1;36m1.0\u001b[0m │\n", "│ \u001b[1;36m-0.322563\u001b[0m │ \u001b[1;36m16.243767\u001b[0m │ \u001b[1;36m1.0\u001b[0m │\n", "│ \u001b[1;36m-0.322461\u001b[0m │ \u001b[1;36m16.243767\u001b[0m │ \u001b[1;36m11.0\u001b[0m │\n", "│ \u001b[1;36m-0.322359\u001b[0m │ \u001b[1;36m16.243767\u001b[0m │ \u001b[1;36m11.0\u001b[0m │\n", "│ \u001b[1;36m-0.322563\u001b[0m │ \u001b[1;36m16.243685\u001b[0m │ \u001b[1;36m1.0\u001b[0m │\n", "│ \u001b[1;36m-0.322461\u001b[0m │ \u001b[1;36m16.243685\u001b[0m │ \u001b[1;36m11.0\u001b[0m │\n", "│ \u001b[1;36m-0.322358\u001b[0m │ \u001b[1;36m16.243685\u001b[0m │ \u001b[1;36m11.0\u001b[0m │\n", "│ \u001b[1;36m-0.322153\u001b[0m │ \u001b[1;36m16.242781\u001b[0m │ \u001b[1;36m2.0\u001b[0m │\n", "└───────────┴───────────┴──────────┘" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print(conn.list_tables())\n", "conn.sql(\"Select * from layer limit 10\")" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "# conn.raw_sql(\"SUMMARIZE layer;\")\n", "# layer.describe()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fine scale data with CL tool" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "# the coarse resolution file\n", "# ! python -m scalenav.rast_converter /Users/cenv1069/Documents/data/datasets/kummu_etal/GDP_PPP_1990_2015_5arcmin_v2.nc outputs/kummu_5arcmin.parquet" ] }, { "cell_type": "code", "execution_count": 71, "metadata": {}, "outputs": [], "source": [ "# the fine resolution file \n", "# At the moment fails\n", "# ! python -m scalenav.rast_converter /Users/cenv1069/Documents/data/datasets/kummu_etal/GDP_PPP_30arcsec_v3.nc outputs/kummu_30arcsec.parquet" ] }, { "cell_type": "code", "execution_count": 100, "metadata": {}, "outputs": [], "source": [ "layer_ = (\n", " layer\n", "# filter region\n", " .filter(_.lon>limits[0]\n", " ,_.lonlimits[1]\n", " ,_.lat 5)\n", " .sample(.1)\n", " .execute(\n", " limit=30_000\n", " )\n", ")" ] }, { "cell_type": "code", "execution_count": 101, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
lonlatband_var
08.46463512.01727613.0
18.46604912.01719513.0
28.46615012.01719513.0
38.46696012.01719513.0
48.46210312.01711313.0
\n", "
" ], "text/plain": [ " lon lat band_var\n", "0 8.464635 12.017276 13.0\n", "1 8.466049 12.017195 13.0\n", "2 8.466150 12.017195 13.0\n", "3 8.466960 12.017195 13.0\n", "4 8.462103 12.017113 13.0" ] }, "execution_count": 101, "metadata": {}, "output_type": "execute_result" } ], "source": [ "layer_.head()" ] }, { "cell_type": "code", "execution_count": 102, "metadata": {}, "outputs": [], "source": [ "# # sample_data = \"data/sample_data_kummu.parquet\"\n", "\n", "# # sample_data = \"data/sample_data_ghsl_msz.parquet\"\n", "\n", "# if not os.path.exists(sample_data):\n", "# layer.to_parquet(sample_data)" ] }, { "cell_type": "code", "execution_count": 103, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(24308, 3)" ] }, "execution_count": 103, "metadata": {}, "output_type": "execute_result" } ], "source": [ "layer_.shape" ] }, { "cell_type": "code", "execution_count": 104, "metadata": {}, "outputs": [], "source": [ "# layer_.count()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Deck map" ] }, { "cell_type": "code", "execution_count": 105, "metadata": {}, "outputs": [], "source": [ "cat_palette = \"Classic_10\"\n", "num_palette = \"Burg\"\n", "\n", "gdp_pc_pal = pypal.load_cmap(cat_palette,reverse=True)" ] }, { "cell_type": "code", "execution_count": 106, "metadata": {}, "outputs": [], "source": [ "deck_data = layer_#.execute(limit=500_000)" ] }, { "cell_type": "code", "execution_count": 107, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
lonlatband_var
08.46463512.01727613.0
18.46604912.01719513.0
28.46615012.01719513.0
38.46696012.01719513.0
48.46210312.01711313.0
\n", "
" ], "text/plain": [ " lon lat band_var\n", "0 8.464635 12.017276 13.0\n", "1 8.466049 12.017195 13.0\n", "2 8.466150 12.017195 13.0\n", "3 8.466960 12.017195 13.0\n", "4 8.462103 12.017113 13.0" ] }, "execution_count": 107, "metadata": {}, "output_type": "execute_result" } ], "source": [ "deck_data.head()" ] }, { "cell_type": "code", "execution_count": 108, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Max input : 24.00, palette colors : 10\n" ] } ], "source": [ "deck_data[\"col\"] = cmap(deck_data.band_var,palette=gdp_pc_pal,log=False)\n", "deck_data[\"value\"] = deck_data[\"band_var\"].apply(lambda x: msz_categories[str(int(x))])\n" ] }, { "cell_type": "code", "execution_count": 109, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(24308, 5)\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
lonlatband_varcolvalue
08.46463512.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...
18.46604912.01719513.0[148, 103, 189, 255]Built spaces, residential building 6m < height...
28.46615012.01719513.0[148, 103, 189, 255]Built spaces, residential building 6m < height...
38.46696012.01719513.0[148, 103, 189, 255]Built spaces, residential building 6m < height...
48.46210312.01711313.0[148, 103, 189, 255]Built spaces, residential building 6m < height...
\n", "
" ], "text/plain": [ " lon lat band_var col \\\n", "0 8.464635 12.017276 13.0 [148, 103, 189, 255] \n", "1 8.466049 12.017195 13.0 [148, 103, 189, 255] \n", "2 8.466150 12.017195 13.0 [148, 103, 189, 255] \n", "3 8.466960 12.017195 13.0 [148, 103, 189, 255] \n", "4 8.462103 12.017113 13.0 [148, 103, 189, 255] \n", "\n", " value \n", "0 Built spaces, residential building 6m < height... \n", "1 Built spaces, residential building 6m < height... \n", "2 Built spaces, residential building 6m < height... \n", "3 Built spaces, residential building 6m < height... \n", "4 Built spaces, residential building 6m < height... " ] }, "execution_count": 109, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print(deck_data.shape)\n", "deck_data.head()" ] }, { "cell_type": "code", "execution_count": 110, "metadata": {}, "outputs": [], "source": [ "viewstate = pdk.data_utils.compute_view(deck_data[[\"lon\",\"lat\"]])" ] }, { "cell_type": "code", "execution_count": 111, "metadata": {}, "outputs": [], "source": [ "\n", "deck_points = pdk.Layer(type=\"ScatterplotLayer\"\n", " ,data=deck_data\n", " ,get_position = [\"lon\",\"lat\"]\n", " ,get_radius=10\n", " ,get_fill_color= \"col\"\n", " ,pickable=True\n", " ,opacity=0.2\n", " ,stroked=True\n", " ,filled=True\n", " ,radius_scale=2\n", " ,radius_min_pixels=3\n", " ,radius_max_pixels=5\n", " ,line_width_min_pixels=1\n", " ,get_line_color=[0, 0, 0, 0]\n", " )" ] }, { "cell_type": "code", "execution_count": 112, "metadata": {}, "outputs": [], "source": [ "deck = pdk.Deck(layers=[deck_points],\n", " initial_view_state=viewstate,\n", " tooltip= {\"text\": \" {band_var} : {value}\"},\n", " map_style=\"road\")" ] }, { "cell_type": "code", "execution_count": 113, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 113, "metadata": {}, "output_type": "execute_result" } ], "source": [ "deck.to_html(\"../../deck_maps/ghsl_msz_vis.html\",iframe_height=700)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gridding function" ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [], "source": [ "# square = snd.square_poly(lat=limits[1],lon=limits[0],distance=grid_param)\n", "\n", "# ref_cell = h3.latlng_to_cell(lat=limits[1],lng=limits[0],res=base_res)\n", "# ref_cell_ij = h3.cell_to_local_ij(origin=ref_cell,h=ref_cell)\n", "\n", "# cells = h3.geo_to_cells(square,res=base_res)\n", "\n", "# cells_ij = [h3.cell_to_local_ij(origin=ref_cell,h=cell) for cell in cells]\n", "\n", "# snd.rast_to_h3[str(grid_param)][\"nn\"] = [(cell_i-ref_cell_ij[0],cell_j-ref_cell_ij[1]) for (cell_i,cell_j) in cells_ij]\n", "\n", "# cells_redone = [h3.local_ij_to_cell(origin=ref_cell,i=cells_i+ref_cell_ij[0],j=cells_j+ref_cell_ij[1]) for (cells_i,cells_j) in snd.rast_to_h3[str(grid_param)][\"nn\"]]\n", "\n", "# gpd.GeoSeries([h3.cells_to_h3shape(cells=cells_redone)],crs=\"epsg:4326\").explore()\n" ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [], "source": [ "# limits = [-1.486,51.138,0.352,51.882]\n", "\n", "# wide europe window\n", "# -14.06\n", "# 48.34\n", "# 31.95\n", "# 59.27" ] }, { "cell_type": "code", "execution_count": 80, "metadata": {}, "outputs": [], "source": [ "layer_.rename(columns={\"lon\":\"x\",\"lat\":\"y\"},inplace=True)" ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
xyband_varcolvalue
08.46493812.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...
18.46655712.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...
28.46665812.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...
38.46706312.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...
48.46665712.01719513.0[148, 103, 189, 255]Built spaces, residential building 6m < height...
\n", "
" ], "text/plain": [ " x y band_var col \\\n", "0 8.464938 12.017276 13.0 [148, 103, 189, 255] \n", "1 8.466557 12.017276 13.0 [148, 103, 189, 255] \n", "2 8.466658 12.017276 13.0 [148, 103, 189, 255] \n", "3 8.467063 12.017276 13.0 [148, 103, 189, 255] \n", "4 8.466657 12.017195 13.0 [148, 103, 189, 255] \n", "\n", " value \n", "0 Built spaces, residential building 6m < height... \n", "1 Built spaces, residential building 6m < height... \n", "2 Built spaces, residential building 6m < height... \n", "3 Built spaces, residential building 6m < height... \n", "4 Built spaces, residential building 6m < height... " ] }, "execution_count": 81, "metadata": {}, "output_type": "execute_result" } ], "source": [ "layer_.head()" ] }, { "cell_type": "code", "execution_count": 82, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
xyband_varcolvalueh3_id
08.46493812.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41831bbaf
18.46655712.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41822cd5f
28.46665812.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41822c8d7
38.46706312.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41822dc6f
48.46665712.01719513.0[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41822c037
\n", "
" ], "text/plain": [ " x y band_var col \\\n", "0 8.464938 12.017276 13.0 [148, 103, 189, 255] \n", "1 8.466557 12.017276 13.0 [148, 103, 189, 255] \n", "2 8.466658 12.017276 13.0 [148, 103, 189, 255] \n", "3 8.467063 12.017276 13.0 [148, 103, 189, 255] \n", "4 8.466657 12.017195 13.0 [148, 103, 189, 255] \n", "\n", " value h3_id \n", "0 Built spaces, residential building 6m < height... 8e580a41831bbaf \n", "1 Built spaces, residential building 6m < height... 8e580a41822cd5f \n", "2 Built spaces, residential building 6m < height... 8e580a41822c8d7 \n", "3 Built spaces, residential building 6m < height... 8e580a41822dc6f \n", "4 Built spaces, residential building 6m < height... 8e580a41822c037 " ] }, "execution_count": 82, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# layer_df = layer.execute(limit=100_000)\n", "# base_res = 14\n", "\n", "layer_df = snd.df_project_on_grid(layer_,res=base_res) \n", "layer_df.head()" ] }, { "cell_type": "code", "execution_count": 83, "metadata": {}, "outputs": [], "source": [ "layer_df.rename(columns={\"band\":\"band_var\"},inplace=True)" ] }, { "cell_type": "code", "execution_count": 84, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(24217, 6)\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
xyband_varcolvalueh3_id
08.46493812.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41831bbaf
18.46655712.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41822cd5f
28.46665812.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41822c8d7
38.46706312.01727613.0[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41822dc6f
48.46665712.01719513.0[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41822c037
\n", "
" ], "text/plain": [ " x y band_var col \\\n", "0 8.464938 12.017276 13.0 [148, 103, 189, 255] \n", "1 8.466557 12.017276 13.0 [148, 103, 189, 255] \n", "2 8.466658 12.017276 13.0 [148, 103, 189, 255] \n", "3 8.467063 12.017276 13.0 [148, 103, 189, 255] \n", "4 8.466657 12.017195 13.0 [148, 103, 189, 255] \n", "\n", " value h3_id \n", "0 Built spaces, residential building 6m < height... 8e580a41831bbaf \n", "1 Built spaces, residential building 6m < height... 8e580a41822cd5f \n", "2 Built spaces, residential building 6m < height... 8e580a41822c8d7 \n", "3 Built spaces, residential building 6m < height... 8e580a41822dc6f \n", "4 Built spaces, residential building 6m < height... 8e580a41822c037 " ] }, "execution_count": 84, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print(layer_df.shape)\n", "layer_df.head()\n" ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [], "source": [ "# layer_df[\"band_var\"].astype(copy=)" ] }, { "cell_type": "code", "execution_count": 86, "metadata": {}, "outputs": [], "source": [ "cat = True\n", "\n", "if cat:\n", " layer_df[\"band_var\"]=layer_df[\"band_var\"].astype(\"int\").astype(\"str\")" ] }, { "cell_type": "code", "execution_count": 87, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(0, 0),\n", " (2, 1),\n", " (1, -1),\n", " (1, 1),\n", " (0, 2),\n", " (-1, 0),\n", " (0, -2),\n", " (-2, -1),\n", " (1, 2),\n", " (2, 0),\n", " (0, 1),\n", " (-1, -1),\n", " (1, -2),\n", " (3, 1),\n", " (1, 0),\n", " (2, 2),\n", " (-1, -2),\n", " (0, -1)]" ] }, "execution_count": 87, "metadata": {}, "output_type": "execute_result" } ], "source": [ "neighbs" ] }, { "cell_type": "code", "execution_count": 88, "metadata": {}, "outputs": [], "source": [ "def layer_gridded(layer : pd.DataFrame, neighbs : list[tuple[int]]): \n", " \"\"\"\" Take a layer of raster centroids expressed as H3 cells and do the 'gridding'. It means filling in all the gaps respecting the original raster grid size.\n", " Duplicates are removed, the value from the single cell is uniformly distributed across the set of cells associated to it. \n", " \n", " \"\"\"\n", " \n", " vars_columns = [x for x in layer.columns if re.search(pattern=\"_var$\",string=x)]\n", "\n", " layer[\"h3_gridded\"] = pd.Series(layer.apply(lambda row: snd.centre_cell_to_square(h3_cell = row[\"h3_id\"],neighbs=neighbs),axis=1).tolist())\n", " layer.dropna(subset=[\"h3_gridded\"],inplace=True)\n", " gridded_layer = layer.explode(\"h3_gridded\")\n", " grid_layer_b = gridded_layer.shape[0]\n", " gridded_layer = gridded_layer[~gridded_layer.duplicated(subset=[\"h3_gridded\"])]\n", " grid_layer_a = gridded_layer.shape[0]\n", " \n", " print(\"Duplicated ratio : \",1-np.round(grid_layer_a/grid_layer_b,3))\n", " \n", " gridded_layer.reset_index(inplace=True,drop=True)\n", " gridded_layer.set_index(\"h3_id\",inplace=True,drop=True)\n", " gridded_layer = gridded_layer.join(gridded_layer.value_counts(\"h3_id\"),how=\"left\",sort=False,validate=\"m:1\")\n", "\n", " # since only one layer is being read at the moment, just check if it's string or not. \n", " # in the future, separate into categorical and numerical variables\n", " if len(vars_columns)>1:\n", " if (str not in gridded_layer[vars_columns].dtypes.values.tolist()):\n", " gridded_layer[vars_columns] = np.round(gridded_layer[vars_columns].div(gridded_layer[\"count\"],axis=0),3)\n", " \n", " # if len(vars_columns)==1:\n", " # if gridded_layer[vars_columns].dtype!=str:\n", " # gridded_layer[vars_columns] = np.round(gridded_layer[vars_columns].div(gridded_layer[\"count\"],axis=0),3)\n", " \n", " gridded_layer.reset_index(inplace=True,drop=False)\n", " gridded_layer.dropna(axis=0,subset=[\"h3_gridded\"],inplace=True)\n", "\n", " return gridded_layer" ] }, { "cell_type": "code", "execution_count": 89, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Duplicated ratio : 0.013000000000000012\n" ] } ], "source": [ "layer_grid = layer_gridded(layer=layer_df,neighbs=neighbs)" ] }, { "cell_type": "code", "execution_count": 90, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(430149, 8)\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
h3_idxyband_varcolvalueh3_griddedcount
08e580a41831bbaf8.46493812.01727613[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41831bbaf18
18e580a41831bbaf8.46493812.01727613[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41831b84f18
28e580a41831bbaf8.46493812.01727613[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41831bb0718
38e580a41831bbaf8.46493812.01727613[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41831bba718
48e580a41831bbaf8.46493812.01727613[148, 103, 189, 255]Built spaces, residential building 6m < height...8e580a41831bb9718
\n", "
" ], "text/plain": [ " h3_id x y band_var col \\\n", "0 8e580a41831bbaf 8.464938 12.017276 13 [148, 103, 189, 255] \n", "1 8e580a41831bbaf 8.464938 12.017276 13 [148, 103, 189, 255] \n", "2 8e580a41831bbaf 8.464938 12.017276 13 [148, 103, 189, 255] \n", "3 8e580a41831bbaf 8.464938 12.017276 13 [148, 103, 189, 255] \n", "4 8e580a41831bbaf 8.464938 12.017276 13 [148, 103, 189, 255] \n", "\n", " value h3_gridded count \n", "0 Built spaces, residential building 6m < height... 8e580a41831bbaf 18 \n", "1 Built spaces, residential building 6m < height... 8e580a41831b84f 18 \n", "2 Built spaces, residential building 6m < height... 8e580a41831bb07 18 \n", "3 Built spaces, residential building 6m < height... 8e580a41831bba7 18 \n", "4 Built spaces, residential building 6m < height... 8e580a41831bb97 18 " ] }, "execution_count": 90, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print(layer_grid.shape)\n", "layer_grid.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mapping gridded" ] }, { "cell_type": "code", "execution_count": 96, "metadata": {}, "outputs": [], "source": [ "sample_size = np.min([layer_grid.shape[0],30_000])\n", "plotting_layer = (layer_grid\n", " # .iloc[0:sample_size]\n", " .sample(sample_size))" ] }, { "cell_type": "code", "execution_count": 97, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Max input : 24.00, palette colors : 10\n" ] } ], "source": [ "plotting_layer[\"color\"] = cmap(plotting_layer[\"band_var\"].astype(\"int\",copy=False)\n", " ,palette=gdp_pc_pal\n", " ,log=False)" ] }, { "cell_type": "code", "execution_count": 98, "metadata": {}, "outputs": [], "source": [ "#####\n", "\n", "layer_gridded = pdk.Layer(\n", " \"H3HexagonLayer\",\n", " plotting_layer,\n", " pickable=True,\n", " stroked=True,\n", " filled=True,\n", " opacity=.5,\n", " extruded=False,\n", " # elevation_scale = 1e-4,\n", " # elevation_range = [0,100],\n", " get_elevation=\"band_var\",\n", " get_hexagon= \"h3_gridded\",\n", " get_fill_color = \"color\",\n", " get_line_color= [0, 0, 0, 100],\n", " line_width_min_pixels=0,\n", " line_width_max_pixels=1,\n", ")\n" ] }, { "cell_type": "code", "execution_count": 99, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 99, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Create mapdeck object\n", "layer_gridded_deck = pdk.Deck(layers=[layer_gridded]\n", " ,initial_view_state=viewstate\n", " ,tooltip= {\"text\": \"{band_var} : {value}\"}\n", " ,map_style=\"road\"\n", " )\n", "\n", "layer_gridded_deck.to_html(\"../../deck_maps/msz_ghsl_gridded.html\",iframe_height=800,open_browser=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.12.6" } }, "nbformat": 4, "nbformat_minor": 2 }