{ "cells": [ { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "#Base\n", "import itertools as iter\n", "import os\n", "import re\n", "import math\n", "import numpy as np\n", "import random as rd\n", "\n", "\n", "# data eng\n", "import h3\n", "import shapely as shp\n", "import geopandas as gpd\n", "import pandas as pd\n", "import rioxarray as rx\n", "import osmnx as ox\n", "import h3\n", "# import ibis as ib\n", "# import polars as pl\n", "\n", "# module\n", "import scalenav.scale_nav as sd\n", "import scalenav.data as dt\n", "\n", "### Visualisation related\n", "from IPython.display import display\n", "import ipywidgets\n", "import pydeck as pdk\n", "from matplotlib import pyplot as plt\n", "import seaborn as sns\n", "import pypalettes as pypal\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Introduction\n", "This document elaborates on the techniques to downscale geospatial data sets across different resolutions. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "nice blog post: https://strimas.com/post/hexagonal-grids/" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "point(-0.1275,51.5131)\n", "\n", "poly(-0.12496,51.51047,-0.12029,51.51305)\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def make_poly(x_min,y_min,x_max,y_max):\n", " return shp.Polygon([(x_min,y_min),(x_max,y_min),(x_max,y_max,),(x_min,y_max)])" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# x = 51.5119\n", "# y = -0.1227\n", "\n", "x_1 = -0.12496\n", "x_2 = -0.12029\n", "y_1 = 51.51047\n", "y_2 = 51.51305\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "square = make_poly(x_1,y_1,x_2,y_2)\n", "square\n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "y,x = square.centroid.x,square.centroid.y\n" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "# res_s2 = 14\n", "\n", "# s2_cell = s2.geo_to_s2(lat=y,lon=x,res=res_s2)\n", "# s2_cell_coords = s2.s2_to_geo_boundary(s2_cell)\n", "# s2_cell_geo = Polygon(s2_cell_coords)\n" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "51.51176 -0.12262499999999998\n" ] } ], "source": [ "print(x,y)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res_h3 = 11\n", "\n", "h3_cell = h3.latlng_to_cell(lng=y,lat=x,res=res_h3)\n", "h3_geo = h3.cells_to_h3shape([h3_cell])\n", "h3_geo" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "43.84" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "h3_d = np.round(np.sqrt(h3.cell_area(h3_cell,unit=\"m^2\")),2)\n", "h3_d\n", "# np.sqrt(h3.cell_area(h=h3_cell,unit=\"m^2\"))" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['8b194ad3245efff',\n", " '8b194ad324edfff',\n", " '8b194ad324e8fff',\n", " '8b194ad324eafff',\n", " '8b194ad324c1fff']" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# relative neighbors\n", "\n", "i_range = range(-4,4)\n", "j_range = range(-4,4)\n", "\n", "coords = [(i,j) for i in i_range for j in j_range]\n", "\n", "ref_cell = h3.cell_to_local_ij(origin=h3_cell,h=h3_cell)\n", "\n", "# coords\n", "neighbs = [h3.local_ij_to_cell(origin=h3_cell,i=ref_cell[0] + id[0],j=ref_cell[1] + id[1]) for id in coords]\n", "neighbs[0:5]" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-5681" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ref_cell[0]" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "neighbs_geo = h3.cells_to_h3shape(neighbs)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "# relative coordinates of the neighbors if approximating the shape with the h3 function \n", "\n", "neighbs = [h3.cell_to_local_ij(origin=h3_cell,h=cell) for cell in h3.geo_to_cells(square,res=11)] \n", "\n", "neighb_id= [(neighb[0]-ref_cell[0],neighb[1]-ref_cell[1]) for neighb in neighbs]\n" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "49\n" ] }, { "data": { "text/plain": [ "[(-2, -4), (3, 4), (-2, -2), (2, 1), (2, 3)]" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print(len(neighb_id))\n", "neighb_id[0:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Approximating other shapes with hexes" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "\n", "rast_to_h3 = {\n", " \"100\": 12,\n", " \"1000\" : 11,\n", " \"10000\" : 9\n", "}\n" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res_h3 = 12\n", "\n", "x = 51.51176 \n", "y = -0.1227\n", "\n", "h3_cell = h3.latlng_to_cell(lng=y,lat=x,res=res_h3)\n", "h3_geo = h3.cells_to_h3shape([h3_cell])\n", "h3_geo\n", "\n", "gpd.GeoSeries([h3_geo],crs=\"EPSG:4326\").explore() # should show a hex over Covent garden in London\n", "\n", "# Getting the ij neighbors h3 ids.\n", "\n", "i_range = range(-4,5)\n", "j_range = range(-4,5)\n", "\n", "coords = [(i,j) for i in i_range for j in j_range]\n", "# coords\n", "ref_cell = h3.cell_to_local_ij(origin=h3_cell,h=h3_cell)\n", "\n", "neighbs = [h3.local_ij_to_cell(origin=h3_cell,i=ref_cell[0]+id[0],j=ref_cell[1]+id[1]) for id in neighb_id]\n", "neighbs\n", "neighbs_geo = h3.cells_to_h3shape(neighbs)\n", "\n", "# plots in the Atlantic Ocean\n", "gpd.GeoSeries([h3_geo,neighbs_geo,square],crs=\"EPSG:4326\").explore()\n" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[(-3, -4),\n", " (-3, -3),\n", " (-3, -2),\n", " (-3, -1),\n", " (-3, 0),\n", " (-3, 1),\n", " (-3, 2),\n", " (-2, -4),\n", " (-2, -3),\n", " (-2, -2),\n", " (-2, -1),\n", " (-2, 0),\n", " (-2, 1),\n", " (-2, 2),\n", " (-1, -3),\n", " (-1, -2),\n", " (-1, -1),\n", " (-1, 0),\n", " (-1, 1),\n", " (-1, 2),\n", " (-1, 3),\n", " (0, -3),\n", " (0, -2),\n", " (0, -1),\n", " (0, 0),\n", " (0, 1),\n", " (0, 2),\n", " (0, 3),\n", " (1, -3),\n", " (1, -2),\n", " (1, -1),\n", " (1, 0),\n", " (1, 1),\n", " (1, 2),\n", " (1, 3),\n", " (2, -2),\n", " (2, -1),\n", " (2, 0),\n", " (2, 1),\n", " (2, 2),\n", " (2, 3),\n", " (2, 4),\n", " (3, -2),\n", " (3, -1),\n", " (3, 0),\n", " (3, 1),\n", " (3, 2),\n", " (3, 3),\n", " (3, 4)]" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[x for x in coords if x in neighb_id]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Naive approach" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "\n", "def plot_df(df, column=None, ax=None):\n", " 'Plot based on the `geometry` column of a GeoPandas dataframe'\n", " df = df.copy()\n", " df = df.to_crs(epsg=3857) # web mercator\n", "\n", " if ax is None:\n", " fig, ax = plt.subplots(figsize=(8,8))\n", " ax.get_xaxis().set_visible(False)\n", " ax.get_yaxis().set_visible(False)\n", " \n", " df.plot(\n", " ax=ax,\n", " alpha=0.5, edgecolor='k',\n", " column=column, categorical=True,\n", " legend=True, legend_kwds={'loc': 'upper left'}, \n", " )\n", " cx.add_basemap(ax, crs=df.crs, source=cx.providers.CartoDB.Positron)\n", "\n", "def plot_shape(shape, ax=None):\n", " df = gpd.GeoDataFrame({'geometry': [shape]}, crs='EPSG:4326')\n", " plot_df(df, ax=ax)\n", "\n", "def plot_cells(cells, ax=None):\n", " shape = h3.cells_to_h3shape(cells)\n", " plot_shape(shape, ax=ax)\n", "\n", "def plot_shape_and_cells(shape, res=9):\n", " fig, axs = plt.subplots(1,2, figsize=(10,5), sharex=True, sharey=True)\n", " plot_shape(shape, ax=axs[0])\n", " plot_cells(h3.h3shape_to_cells(shape, res), ax=axs[1])\n", " fig.tight_layout()" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "child_num = 7\n", "res_ = 3" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "POLYGON ((10.371 0.725, 10.371 18.208, -15.161 18.208, -15.161 0.725, 10.371 0.725))\n" ] } ], "source": [ "test_grid = [-15.161,0.725,10.371,18.208]\n", "\n", "bbox_polygon = shp.geometry.box(*test_grid)\n", "print(bbox_polygon.wkt)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "h3_shape = h3.geo_to_h3shape(bbox_polygon)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "h3_shape_cells = h3.h3shape_to_cells(h3_shape,res=res_)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "h3_df = pd.DataFrame(pd.Series(h3_shape_cells),columns=[\"h3_id\"],dtype='str')\n", "\n", "# with polars\n", "# h3_pl = pl.DataFrame(h3_df)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "530" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(h3_shape_cells)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "h3_df[\"geom\"] = gpd.GeoSeries(h3_df.h3_id.apply(lambda x: h3.cells_to_h3shape([x])))\n", "h3_df = gpd.GeoDataFrame(h3_df,geometry=\"geom\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data manipulations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "adding some randomly sampled data" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [], "source": [ "alpha = 1\n", "h3_df[\"values_e\"] = [rd.expovariate(lambd=2*alpha) for x in h3_df.h3_id]\n", "\n", "h3_df[\"values_i\"] = [rd.expovariate(lambd=alpha) for x in h3_df.h3_id]\n", "\n", "h3_df = sd.add_geom(h3_df)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "h3_df.explore(column=\"values_e\")" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [], "source": [ "h3_df[\"child\"] = h3_df.h3_id.apply(h3.cell_to_children)\n" ] }, { "cell_type": "code", "execution_count": 40, "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", "
h3_idvalues_ivalues_e
08454041ffffffff0.602140.180218
08454043ffffffff0.602140.180218
08454045ffffffff0.602140.180218
08454047ffffffff0.602140.180218
08454049ffffffff0.602140.180218
\n", "
" ], "text/plain": [ " h3_id values_i values_e\n", "0 8454041ffffffff 0.60214 0.180218\n", "0 8454043ffffffff 0.60214 0.180218\n", "0 8454045ffffffff 0.60214 0.180218\n", "0 8454047ffffffff 0.60214 0.180218\n", "0 8454049ffffffff 0.60214 0.180218" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "h3_df_expl = h3_df[[\"child\",\"values_i\",\"values_e\"]].explode(\"child\").rename(columns={\"child\":\"h3_id\"})\n", "h3_df_expl.head()" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": [ "h3_df_expl[\"values_e\"] = h3_df_expl[[\"h3_id\",\"values_e\"]].filter(regex=\"_e\").mul({\"values_e\":1/child_num})" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [], "source": [ "h3_df_expl = sd.add_geom(h3_df_expl)" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "h3_df_expl.explore(column = \"values_e\")" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "osm_admin = ox.geocode_to_gdf(query=\"Ile de France\")\n", "# osm_poly = h3.latlng_to_cell(osm_admin,res=res_)\n", "osm_admin.explore()" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "osm_df = pd.DataFrame()" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 861fb6b4fffffff\n", "1 861fb4607ffffff\n", "2 861fb0ae7ffffff\n", "3 861fb5497ffffff\n", "4 861fb6a6fffffff\n", " ... \n", "366 861fb6b07ffffff\n", "367 861fb0197ffffff\n", "368 861fb0a9fffffff\n", "369 861fb6a27ffffff\n", "370 861fb44dfffffff\n", "Name: geometry, Length: 371, dtype: object" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "osm_cells = osm_admin.geometry.apply(lambda x: h3.geo_to_cells(x,res=6)).explode().reset_index(drop=True)\n", "osm_cells\n" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [], "source": [ "osm_df[\"compact\"] = pd.Series(h3.compact_cells(osm_cells),index=None)\n" ] }, { "cell_type": "code", "execution_count": 48, "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", "
compact
0861fb0ae7ffffff
1861fb636fffffff
2861fb0307ffffff
3861fb448fffffff
4861fb0067ffffff
\n", "
" ], "text/plain": [ " compact\n", "0 861fb0ae7ffffff\n", "1 861fb636fffffff\n", "2 861fb0307ffffff\n", "3 861fb448fffffff\n", "4 861fb0067ffffff" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "osm_df.head()" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [], "source": [ "osm_df[\"geom\"] = osm_df.compact.apply(lambda x: h3.cells_to_h3shape([x]))\n", "\n", "osm_df = gpd.GeoDataFrame(osm_df,geometry=\"geom\",crs=\"EPSG:4326\")" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "osm_df.explore()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Projecting crops " ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [], "source": [ "mapspam_dir = \"/Users/cenv1069/Documents/agriculture/mapspam/spam2017v2r1_ssa_val_prod/\"\n", "\n", "mapspan_files = os.listdir(mapspam_dir)" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [], "source": [ "mapspam = rx.open_rasterio(mapspam_dir+mapspan_files[0],)" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (band: 1, y: 2160, x: 4320)> Size: 37MB\n",
       "[9331200 values with dtype=float32]\n",
       "Coordinates:\n",
       "  * band         (band) int64 8B 1\n",
       "  * x            (x) float64 35kB -180.0 -179.9 -179.8 ... 179.8 179.9 180.0\n",
       "  * y            (y) float64 17kB 89.96 89.88 89.79 ... -89.79 -89.87 -89.96\n",
       "    spatial_ref  int64 8B 0\n",
       "Attributes:\n",
       "    AREA_OR_POINT:  Area\n",
       "    _FillValue:     -1.0\n",
       "    scale_factor:   1.0\n",
       "    add_offset:     0.0
" ], "text/plain": [ " Size: 37MB\n", "[9331200 values with dtype=float32]\n", "Coordinates:\n", " * band (band) int64 8B 1\n", " * x (x) float64 35kB -180.0 -179.9 -179.8 ... 179.8 179.9 180.0\n", " * y (y) float64 17kB 89.96 89.88 89.79 ... -89.79 -89.87 -89.96\n", " spatial_ref int64 8B 0\n", "Attributes:\n", " AREA_OR_POINT: Area\n", " _FillValue: -1.0\n", " scale_factor: 1.0\n", " add_offset: 0.0" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mapspam" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [], "source": [ "mapspam_df_ = mapspam[0,:,:].to_pandas().melt(ignore_index=False).reset_index().rename(columns={\"value\":\"value_var\"})" ] }, { "cell_type": "code", "execution_count": 55, "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", "
yxvalue_var
089.958333-179.958334-1.0
189.875000-179.958334-1.0
289.791667-179.958334-1.0
389.708334-179.958334-1.0
489.625001-179.958334-1.0
\n", "
" ], "text/plain": [ " y x value_var\n", "0 89.958333 -179.958334 -1.0\n", "1 89.875000 -179.958334 -1.0\n", "2 89.791667 -179.958334 -1.0\n", "3 89.708334 -179.958334 -1.0\n", "4 89.625001 -179.958334 -1.0" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mapspam_df_.head()" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [], "source": [ "mapspam_df_raw = mapspam_df_[mapspam_df_.value_var>0].reset_index(drop=True)\n", "# mapspam_df = mapspam_df_[mapspam_df_.value_var>0].reset_index(drop=True)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGvCAYAAABSC3+tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAueUlEQVR4nO3df1RU953/8deEHyMQmBUp4JxAlqSE1YKpwXwRTapWBV3Fetyt3SWdNV1rzNFoWOVka92exa2BrL9iVzYedT1Rgy795ljT1Bg65KQhJfiTLCeiHps0fv3RgJiIoOgOE5jvHz3ekxFjZhQc+OT5OIdzMp/7nvv53PeEk1c+M5ex+Xw+nwAAAAx0T6gXAAAA0FcIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAY4WHegGh1N3drU8++USxsbGy2WyhXg4AAAiAz+fT5cuX5XQ6dc89t96z+VoHnU8++UQpKSmhXgYAALgNZ8+e1X333XfLmq910ImNjZX050bFxcX16rm9Xq/cbrfy8vIUERHRq+c2CX0KDH0KDH0KDH0KDH0KTCj61N7erpSUFOu/47fytQ4619+uiouL65OgEx0drbi4OH5BboE+BYY+BYY+BYY+BYY+BSaUfQrkYyd8GBkAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWOGhXoDpMkt+K0/XV3+NfH/x/16YFuolAADQa9jRAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADBWUEGnpKRENpvN7yc5Odk67vP5VFJSIqfTqaioKI0fP17Hjh3zO4fH49GiRYuUkJCgmJgYzZgxQ+fOnfOraW1tlcvlksPhkMPhkMvl0qVLl/xqzpw5o4KCAsXExCghIUGLFy9WZ2dnkJcPAABMFvSOzre+9S01NTVZP0ePHrWOrVq1SuvWrVN5ebkOHz6s5ORkTZ48WZcvX7ZqioqKtGfPHlVWVqq2tlZXrlzR9OnT1dXVZdUUFhaqoaFBVVVVqqqqUkNDg1wul3W8q6tL06ZNU0dHh2pra1VZWandu3dr6dKlt9sHAABgoPCgnxAe7reLc53P59P69eu1fPlyzZo1S5K0fft2JSUladeuXZo/f77a2tq0detWvfLKK5o0aZIkqaKiQikpKXrrrbeUn5+vEydOqKqqSgcOHFBOTo4kacuWLcrNzdXJkyeVkZEht9ut48eP6+zZs3I6nZKktWvX6sknn9Tzzz+vuLi4224IAAAwR9BB58MPP5TT6ZTdbldOTo5KS0v1wAMP6NSpU2publZeXp5Va7fbNW7cONXV1Wn+/Pmqr6+X1+v1q3E6ncrMzFRdXZ3y8/O1f/9+ORwOK+RI0ujRo+VwOFRXV6eMjAzt379fmZmZVsiRpPz8fHk8HtXX12vChAk3XbvH45HH47Eet7e3S5K8Xq+8Xm+wrbil6+ez3+Pr1fP2td7uQ6Dz3e15Bxr6FBj6FBj6FBj6FJhQ9CmYuYIKOjk5OdqxY4ceeughnT9/XitXrtSYMWN07NgxNTc3S5KSkpL8npOUlKTTp09LkpqbmxUZGanBgwf3qLn+/ObmZiUmJvaYOzEx0a/mxnkGDx6syMhIq+ZmysrKtGLFih7jbrdb0dHRX3X5t+Xno7r75Lx9Zd++fSGZt7q6OiTzDjT0KTD0KTD0KTD0KTB3s09Xr14NuDaooDN16lTrn7OyspSbm6sHH3xQ27dv1+jRoyVJNpvN7zk+n6/H2I1urLlZ/e3U3GjZsmVasmSJ9bi9vV0pKSnKy8vr9be7vF6vqqur9bMj98jTfevr708aS/Lv6nzX+zR58mRFRETc1bkHEvoUGPoUGPoUGPoUmFD06fo7MoEI+q2rL4qJiVFWVpY+/PBDzZw5U9Kfd1uGDh1q1bS0tFi7L8nJyers7FRra6vfrk5LS4vGjBlj1Zw/f77HXBcuXPA7z8GDB/2Ot7a2yuv19tjp+SK73S673d5jPCIios9eHE+3TZ6ugRN0QvXL3JevgUnoU2DoU2DoU2DoU2DuZp+CmeeO/o6Ox+PRiRMnNHToUKWlpSk5Odlv66qzs1M1NTVWiMnOzlZERIRfTVNTkxobG62a3NxctbW16dChQ1bNwYMH1dbW5lfT2NiopqYmq8btdstutys7O/tOLgkAABgkqB2d4uJiFRQUKDU1VS0tLVq5cqXa29s1Z84c2Ww2FRUVqbS0VOnp6UpPT1dpaamio6NVWFgoSXI4HJo7d66WLl2qIUOGKD4+XsXFxcrKyrLuwho2bJimTJmiefPmadOmTZKkp556StOnT1dGRoYkKS8vT8OHD5fL5dLq1at18eJFFRcXa968edxxBQAALEEFnXPnzunv//7v9emnn+ob3/iGRo8erQMHDuj++++XJD333HO6du2aFixYoNbWVuXk5Mjtdis2NtY6x4svvqjw8HDNnj1b165d08SJE7Vt2zaFhYVZNTt37tTixYutu7NmzJih8vJy63hYWJjeeOMNLViwQGPHjlVUVJQKCwu1Zs2aO2oGAAAwS1BBp7Ky8pbHbTabSkpKVFJS8qU1gwYN0oYNG7Rhw4YvrYmPj1dFRcUt50pNTdXevXtvWQMAAL7e+K4rAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYdxR0ysrKZLPZVFRUZI35fD6VlJTI6XQqKipK48eP17Fjx/ye5/F4tGjRIiUkJCgmJkYzZszQuXPn/GpaW1vlcrnkcDjkcDjkcrl06dIlv5ozZ86ooKBAMTExSkhI0OLFi9XZ2XknlwQAAAxy20Hn8OHD2rx5s0aMGOE3vmrVKq1bt07l5eU6fPiwkpOTNXnyZF2+fNmqKSoq0p49e1RZWana2lpduXJF06dPV1dXl1VTWFiohoYGVVVVqaqqSg0NDXK5XNbxrq4uTZs2TR0dHaqtrVVlZaV2796tpUuX3u4lAQAAw9xW0Lly5YqeeOIJbdmyRYMHD7bGfT6f1q9fr+XLl2vWrFnKzMzU9u3bdfXqVe3atUuS1NbWpq1bt2rt2rWaNGmSRo4cqYqKCh09elRvvfWWJOnEiROqqqrSf/3Xfyk3N1e5ubnasmWL9u7dq5MnT0qS3G63jh8/roqKCo0cOVKTJk3S2rVrtWXLFrW3t99pXwAAgAHCb+dJCxcu1LRp0zRp0iStXLnSGj916pSam5uVl5dnjdntdo0bN051dXWaP3++6uvr5fV6/WqcTqcyMzNVV1en/Px87d+/Xw6HQzk5OVbN6NGj5XA4VFdXp4yMDO3fv1+ZmZlyOp1WTX5+vjwej+rr6zVhwoQe6/Z4PPJ4PNbj64HI6/XK6/XeTiu+1PXz2e/x9ep5+1pv9yHQ+e72vAMNfQoMfQoMfQoMfQpMKPoUzFxBB53Kykq9//77Onz4cI9jzc3NkqSkpCS/8aSkJJ0+fdqqiYyM9NsJul5z/fnNzc1KTEzscf7ExES/mhvnGTx4sCIjI62aG5WVlWnFihU9xt1ut6Kjo2/6nDv181HdfXLevrJv376QzFtdXR2SeQca+hQY+hQY+hQY+hSYu9mnq1evBlwbVNA5e/asnn32Wbndbg0aNOhL62w2m99jn8/XY+xGN9bcrP52ar5o2bJlWrJkifW4vb1dKSkpysvLU1xc3C3XFyyv16vq6mr97Mg98nTf+tr7k8aS/Ls63/U+TZ48WREREXd17oGEPgWGPgWGPgWGPgUmFH0K5iMqQQWd+vp6tbS0KDs72xrr6urSu+++q/LycuvzM83NzRo6dKhV09LSYu2+JCcnq7OzU62trX67Oi0tLRozZoxVc/78+R7zX7hwwe88Bw8e9Dve2toqr9fbY6fnOrvdLrvd3mM8IiKiz14cT7dNnq6BE3RC9cvcl6+BSehTYOhTYOhTYOhTYO5mn4KZJ6gPI0+cOFFHjx5VQ0OD9TNq1Cg98cQTamho0AMPPKDk5GS/7avOzk7V1NRYISY7O1sRERF+NU1NTWpsbLRqcnNz1dbWpkOHDlk1Bw8eVFtbm19NY2OjmpqarBq32y273e4XxAAAwNdXUDs6sbGxyszM9BuLiYnRkCFDrPGioiKVlpYqPT1d6enpKi0tVXR0tAoLCyVJDodDc+fO1dKlSzVkyBDFx8eruLhYWVlZmjRpkiRp2LBhmjJliubNm6dNmzZJkp566ilNnz5dGRkZkqS8vDwNHz5cLpdLq1ev1sWLF1VcXKx58+b1+ttQAABgYLqtu65u5bnnntO1a9e0YMECtba2KicnR263W7GxsVbNiy++qPDwcM2ePVvXrl3TxIkTtW3bNoWFhVk1O3fu1OLFi627s2bMmKHy8nLreFhYmN544w0tWLBAY8eOVVRUlAoLC7VmzZreviQAADBA3XHQeeedd/we22w2lZSUqKSk5EufM2jQIG3YsEEbNmz40pr4+HhVVFTccu7U1FTt3bs3mOUCAICvEb7rCgAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgrKCCzsaNGzVixAjFxcUpLi5Oubm5evPNN63jPp9PJSUlcjqdioqK0vjx43Xs2DG/c3g8Hi1atEgJCQmKiYnRjBkzdO7cOb+a1tZWuVwuORwOORwOuVwuXbp0ya/mzJkzKigoUExMjBISErR48WJ1dnYGefkAAMBkQQWd++67Ty+88IKOHDmiI0eO6Lvf/a6+973vWWFm1apVWrduncrLy3X48GElJydr8uTJunz5snWOoqIi7dmzR5WVlaqtrdWVK1c0ffp0dXV1WTWFhYVqaGhQVVWVqqqq1NDQIJfLZR3v6urStGnT1NHRodraWlVWVmr37t1aunTpnfYDAAAYJDyY4oKCAr/Hzz//vDZu3KgDBw5o+PDhWr9+vZYvX65Zs2ZJkrZv366kpCTt2rVL8+fPV1tbm7Zu3apXXnlFkyZNkiRVVFQoJSVFb731lvLz83XixAlVVVXpwIEDysnJkSRt2bJFubm5OnnypDIyMuR2u3X8+HGdPXtWTqdTkrR27Vo9+eSTev755xUXF3fHjQEAAANfUEHni7q6uvTqq6+qo6NDubm5OnXqlJqbm5WXl2fV2O12jRs3TnV1dZo/f77q6+vl9Xr9apxOpzIzM1VXV6f8/Hzt379fDofDCjmSNHr0aDkcDtXV1SkjI0P79+9XZmamFXIkKT8/Xx6PR/X19ZowYcJN1+zxeOTxeKzH7e3tkiSv1yuv13u7rbip6+ez3+Pr1fP2td7uQ6Dz3e15Bxr6FBj6FBj6FBj6FJhQ9CmYuYIOOkePHlVubq7+93//V/fee6/27Nmj4cOHq66uTpKUlJTkV5+UlKTTp09LkpqbmxUZGanBgwf3qGlubrZqEhMTe8ybmJjoV3PjPIMHD1ZkZKRVczNlZWVasWJFj3G3263o6OivuvTb8vNR3X1y3r6yb9++kMxbXV0dknkHGvoUGPoUGPoUGPoUmLvZp6tXrwZcG3TQycjIUENDgy5duqTdu3drzpw5qqmpsY7bbDa/ep/P12PsRjfW3Kz+dmputGzZMi1ZssR63N7erpSUFOXl5fX6211er1fV1dX62ZF75Om+9fX3J40l+Xd1vut9mjx5siIiIu7q3AMJfQoMfQoMfQoMfQpMKPp0/R2ZQAQddCIjI/XNb35TkjRq1CgdPnxYv/jFL/TP//zPkv682zJ06FCrvqWlxdp9SU5OVmdnp1pbW/12dVpaWjRmzBir5vz58z3mvXDhgt95Dh486He8tbVVXq+3x07PF9ntdtnt9h7jERERffbieLpt8nQNnKATql/mvnwNTEKfAkOfAkOfAkOfAnM3+xTMPHf8d3R8Pp88Ho/S0tKUnJzst3XV2dmpmpoaK8RkZ2crIiLCr6apqUmNjY1WTW5urtra2nTo0CGr5uDBg2pra/OraWxsVFNTk1Xjdrtlt9uVnZ19p5cEAAAMEdSOzk9/+lNNnTpVKSkpunz5siorK/XOO++oqqpKNptNRUVFKi0tVXp6utLT01VaWqro6GgVFhZKkhwOh+bOnaulS5dqyJAhio+PV3FxsbKysqy7sIYNG6YpU6Zo3rx52rRpkyTpqaee0vTp05WRkSFJysvL0/Dhw+VyubR69WpdvHhRxcXFmjdvHndcAQAAS1BB5/z583K5XGpqapLD4dCIESNUVVWlyZMnS5Kee+45Xbt2TQsWLFBra6tycnLkdrsVGxtrnePFF19UeHi4Zs+erWvXrmnixInatm2bwsLCrJqdO3dq8eLF1t1ZM2bMUHl5uXU8LCxMb7zxhhYsWKCxY8cqKipKhYWFWrNmzR01AwAAmCWooLN169ZbHrfZbCopKVFJScmX1gwaNEgbNmzQhg0bvrQmPj5eFRUVt5wrNTVVe/fuvWUNAAD4euO7rgAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYwUVdMrKyvToo48qNjZWiYmJmjlzpk6ePOlX4/P5VFJSIqfTqaioKI0fP17Hjh3zq/F4PFq0aJESEhIUExOjGTNm6Ny5c341ra2tcrlccjgccjgccrlcunTpkl/NmTNnVFBQoJiYGCUkJGjx4sXq7OwM5pIAAIDBggo6NTU1WrhwoQ4cOKDq6mp9/vnnysvLU0dHh1WzatUqrVu3TuXl5Tp8+LCSk5M1efJkXb582aopKirSnj17VFlZqdraWl25ckXTp09XV1eXVVNYWKiGhgZVVVWpqqpKDQ0Ncrlc1vGuri5NmzZNHR0dqq2tVWVlpXbv3q2lS5feST8AAIBBwoMprqqq8nv88ssvKzExUfX19frOd74jn8+n9evXa/ny5Zo1a5Ykafv27UpKStKuXbs0f/58tbW1aevWrXrllVc0adIkSVJFRYVSUlL01ltvKT8/XydOnFBVVZUOHDignJwcSdKWLVuUm5urkydPKiMjQ263W8ePH9fZs2fldDolSWvXrtWTTz6p559/XnFxcXfcHAAAMLAFFXRu1NbWJkmKj4+XJJ06dUrNzc3Ky8uzaux2u8aNG6e6ujrNnz9f9fX18nq9fjVOp1OZmZmqq6tTfn6+9u/fL4fDYYUcSRo9erQcDofq6uqUkZGh/fv3KzMz0wo5kpSfny+Px6P6+npNmDChx3o9Ho88Ho/1uL29XZLk9Xrl9XrvpBU9XD+f/R5fr563r/V2HwKd727PO9DQp8DQp8DQp8DQp8CEok/BzHXbQcfn82nJkiV67LHHlJmZKUlqbm6WJCUlJfnVJiUl6fTp01ZNZGSkBg8e3KPm+vObm5uVmJjYY87ExES/mhvnGTx4sCIjI62aG5WVlWnFihU9xt1ut6Kjo7/ymm/Hz0d198l5+8q+fftCMm91dXVI5h1o6FNg6FNg6FNg6FNg7mafrl69GnDtbQedZ555Rh988IFqa2t7HLPZbH6PfT5fj7Eb3Vhzs/rbqfmiZcuWacmSJdbj9vZ2paSkKC8vr9ff6vJ6vaqurtbPjtwjT/etr70/aSzJv6vzXe/T5MmTFRERcVfnHkjoU2DoU2DoU2DoU2BC0afr78gE4raCzqJFi/T666/r3Xff1X333WeNJycnS/rzbsvQoUOt8ZaWFmv3JTk5WZ2dnWptbfXb1WlpadGYMWOsmvPnz/eY98KFC37nOXjwoN/x1tZWeb3eHjs919ntdtnt9h7jERERffbieLpt8nQNnKATql/mvnwNTEKfAkOfAkOfAkOfAnM3+xTMPEHddeXz+fTMM8/oV7/6ld5++22lpaX5HU9LS1NycrLf9lVnZ6dqamqsEJOdna2IiAi/mqamJjU2Nlo1ubm5amtr06FDh6yagwcPqq2tza+msbFRTU1NVo3b7Zbdbld2dnYwlwUAAAwV1I7OwoULtWvXLv36179WbGys9VkYh8OhqKgo2Ww2FRUVqbS0VOnp6UpPT1dpaamio6NVWFho1c6dO1dLly7VkCFDFB8fr+LiYmVlZVl3YQ0bNkxTpkzRvHnztGnTJknSU089penTpysjI0OSlJeXp+HDh8vlcmn16tW6ePGiiouLNW/ePO64AgAAkoIMOhs3bpQkjR8/3m/85Zdf1pNPPilJeu6553Tt2jUtWLBAra2tysnJkdvtVmxsrFX/4osvKjw8XLNnz9a1a9c0ceJEbdu2TWFhYVbNzp07tXjxYuvurBkzZqi8vNw6HhYWpjfeeEMLFizQ2LFjFRUVpcLCQq1ZsyaoBgAAAHMFFXR8vq++Vdpms6mkpEQlJSVfWjNo0CBt2LBBGzZs+NKa+Ph4VVRU3HKu1NRU7d279yvXBAAAvp74risAAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFhBB513331XBQUFcjqdstlseu211/yO+3w+lZSUyOl0KioqSuPHj9exY8f8ajwejxYtWqSEhATFxMRoxowZOnfunF9Na2urXC6XHA6HHA6HXC6XLl265Fdz5swZFRQUKCYmRgkJCVq8eLE6OzuDvSQAAGCooINOR0eHHn74YZWXl9/0+KpVq7Ru3TqVl5fr8OHDSk5O1uTJk3X58mWrpqioSHv27FFlZaVqa2t15coVTZ8+XV1dXVZNYWGhGhoaVFVVpaqqKjU0NMjlclnHu7q6NG3aNHV0dKi2tlaVlZXavXu3li5dGuwlAQAAQ4UH+4SpU6dq6tSpNz3m8/m0fv16LV++XLNmzZIkbd++XUlJSdq1a5fmz5+vtrY2bd26Va+88oomTZokSaqoqFBKSoreeust5efn68SJE6qqqtKBAweUk5MjSdqyZYtyc3N18uRJZWRkyO126/jx4zp79qycTqckae3atXryySf1/PPPKy4u7rYaAgAAzBF00LmVU6dOqbm5WXl5edaY3W7XuHHjVFdXp/nz56u+vl5er9evxul0KjMzU3V1dcrPz9f+/fvlcDiskCNJo0ePlsPhUF1dnTIyMrR//35lZmZaIUeS8vPz5fF4VF9frwkTJvRYn8fjkcfjsR63t7dLkrxer7xeb2+2wjqf/R5fr563r/V2HwKd727PO9DQp8DQp8DQp8DQp8CEok/BzNWrQae5uVmSlJSU5DeelJSk06dPWzWRkZEaPHhwj5rrz29ublZiYmKP8ycmJvrV3DjP4MGDFRkZadXcqKysTCtWrOgx7na7FR0dHcglBu3no7r75Lx9Zd++fSGZt7q6OiTzDjT0KTD0KTD0KTD0KTB3s09Xr14NuLZXg851NpvN77HP5+sxdqMba25Wfzs1X7Rs2TItWbLEetze3q6UlBTl5eX1+ltdXq9X1dXV+tmRe+TpvvW19yeNJfl3db7rfZo8ebIiIiLu6twDCX0KDH0KDH0KDH0KTCj6dP0dmUD0atBJTk6W9OfdlqFDh1rjLS0t1u5LcnKyOjs71dra6rer09LSojFjxlg158+f73H+Cxcu+J3n4MGDfsdbW1vl9Xp77PRcZ7fbZbfbe4xHRET02Yvj6bbJ0zVwgk6ofpn78jUwCX0KDH0KDH0KDH0KzN3sUzDz9Orf0UlLS1NycrLf9lVnZ6dqamqsEJOdna2IiAi/mqamJjU2Nlo1ubm5amtr06FDh6yagwcPqq2tza+msbFRTU1NVo3b7Zbdbld2dnZvXhYAABiggt7RuXLlij766CPr8alTp9TQ0KD4+HilpqaqqKhIpaWlSk9PV3p6ukpLSxUdHa3CwkJJksPh0Ny5c7V06VINGTJE8fHxKi4uVlZWlnUX1rBhwzRlyhTNmzdPmzZtkiQ99dRTmj59ujIyMiRJeXl5Gj58uFwul1avXq2LFy+quLhY8+bN444rAAAg6TaCzpEjR/zuaLr+mZc5c+Zo27Zteu6553Tt2jUtWLBAra2tysnJkdvtVmxsrPWcF198UeHh4Zo9e7auXbumiRMnatu2bQoLC7Nqdu7cqcWLF1t3Z82YMcPvb/eEhYXpjTfe0IIFCzR27FhFRUWpsLBQa9asCb4LAADASEEHnfHjx8vn+/Jbpm02m0pKSlRSUvKlNYMGDdKGDRu0YcOGL62Jj49XRUXFLdeSmpqqvXv3fuWaAQDA1xPfdQUAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMNaADzovvfSS0tLSNGjQIGVnZ+v3v/99qJcEAAD6ifBQL+BO/PKXv1RRUZFeeukljR07Vps2bdLUqVN1/Phxpaamhnp5A9Jf/uSNuzqfPcynVf9Hyiz5rTxdtts6x/97YVovrwoAYIoBvaOzbt06zZ07Vz/+8Y81bNgwrV+/XikpKdq4cWOolwYAAPqBAbuj09nZqfr6ev3kJz/xG8/Ly1NdXd1Nn+PxeOTxeKzHbW1tkqSLFy/K6/X26vq8Xq+uXr2qcO896uq+vZ2Kr4Pwbp+uXu2+oz59s/j/9vKq+t7BZRODqr/+79Nnn32miIiIPlrVwEefAkOfAkOfAhOKPl2+fFmS5PP5vrJ2wAadTz/9VF1dXUpKSvIbT0pKUnNz802fU1ZWphUrVvQYT0tL65M1IjCFoV5ACCSsDfUKAGDgu3z5shwOxy1rBmzQuc5m898F8Pl8PcauW7ZsmZYsWWI97u7u1sWLFzVkyJAvfc7tam9vV0pKis6ePau4uLhePbdJ6FNg6FNg6FNg6FNg6FNgQtEnn8+ny5cvy+l0fmXtgA06CQkJCgsL67F709LS0mOX5zq73S673e439hd/8Rd9tURJUlxcHL8gAaBPgaFPgaFPgaFPgaFPgbnbffqqnZzrBuyHkSMjI5Wdna3q6mq/8erqao0ZMyZEqwIAAP3JgN3RkaQlS5bI5XJp1KhRys3N1ebNm3XmzBk9/fTToV4aAADoBwZ00PnBD36gzz77TP/2b/+mpqYmZWZmat++fbr//vtDvTTZ7Xb967/+a4+3yuCPPgWGPgWGPgWGPgWGPgWmv/fJ5gvk3iwAAIABaMB+RgcAAOCrEHQAAICxCDoAAMBYBB0AAGAsgk4feOmll5SWlqZBgwYpOztbv//970O9pH7n3XffVUFBgZxOp2w2m1577bVQL6nfKSsr06OPPqrY2FglJiZq5syZOnnyZKiX1e9s3LhRI0aMsP5YWW5urt58881QL6vfKysrk81mU1FRUaiX0q+UlJTIZrP5/SQnJ4d6Wf3Sn/70J/3whz/UkCFDFB0drW9/+9uqr68P9bJ6IOj0sl/+8pcqKirS8uXL9T//8z96/PHHNXXqVJ05cybUS+tXOjo69PDDD6u8vDzUS+m3ampqtHDhQh04cEDV1dX6/PPPlZeXp46OjlAvrV+577779MILL+jIkSM6cuSIvvvd7+p73/uejh07Fuql9VuHDx/W5s2bNWLEiFAvpV/61re+paamJuvn6NGjoV5Sv9Pa2qqxY8cqIiJCb775po4fP661a9f2+bcN3A5uL+9lOTk5euSRR7Rx40ZrbNiwYZo5c6bKyspCuLL+y2azac+ePZo5c2aol9KvXbhwQYmJiaqpqdF3vvOdUC+nX4uPj9fq1as1d+7cUC+l37ly5YoeeeQRvfTSS1q5cqW+/e1va/369aFeVr9RUlKi1157TQ0NDaFeSr/2k5/8RO+9996AeMeCHZ1e1NnZqfr6euXl5fmN5+Xlqa6uLkSrgina2tok/fk/4ri5rq4uVVZWqqOjQ7m5uaFeTr+0cOFCTZs2TZMmTQr1UvqtDz/8UE6nU2lpafq7v/s7ffzxx6FeUr/z+uuva9SoUfr+97+vxMREjRw5Ulu2bAn1sm6KoNOLPv30U3V1dfX4UtGkpKQeXz4KBMPn82nJkiV67LHHlJmZGerl9DtHjx7VvffeK7vdrqefflp79uzR8OHDQ72sfqeyslLvv/8+u8u3kJOTox07dui3v/2ttmzZoubmZo0ZM0afffZZqJfWr3z88cfauHGj0tPT9dvf/lZPP/20Fi9erB07doR6aT0M6K+A6K9sNpvfY5/P12MMCMYzzzyjDz74QLW1taFeSr+UkZGhhoYGXbp0Sbt379acOXNUU1ND2PmCs2fP6tlnn5Xb7dagQYNCvZx+a+rUqdY/Z2VlKTc3Vw8++KC2b9+uJUuWhHBl/Ut3d7dGjRql0tJSSdLIkSN17Ngxbdy4Uf/wD/8Q4tX5Y0enFyUkJCgsLKzH7k1LS0uPXR4gUIsWLdLrr7+u3/3ud7rvvvtCvZx+KTIyUt/85jc1atQolZWV6eGHH9YvfvGLUC+rX6mvr1dLS4uys7MVHh6u8PBw1dTU6D/+4z8UHh6urq6uUC+xX4qJiVFWVpY+/PDDUC+lXxk6dGiP/5EYNmxYv7zxhqDTiyIjI5Wdna3q6mq/8erqao0ZMyZEq8JA5fP59Mwzz+hXv/qV3n77baWlpYV6SQOGz+eTx+MJ9TL6lYkTJ+ro0aNqaGiwfkaNGqUnnnhCDQ0NCgsLC/US+yWPx6MTJ05o6NChoV5KvzJ27Ngef+7iD3/4Q7/4Uu0b8dZVL1uyZIlcLpdGjRql3Nxcbd68WWfOnNHTTz8d6qX1K1euXNFHH31kPT516pQaGhoUHx+v1NTUEK6s/1i4cKF27dqlX//614qNjbV2Ch0Oh6KiokK8uv7jpz/9qaZOnaqUlBRdvnxZlZWVeuedd1RVVRXqpfUrsbGxPT7fFRMToyFDhvC5ry8oLi5WQUGBUlNT1dLSopUrV6q9vV1z5swJ9dL6lX/6p3/SmDFjVFpaqtmzZ+vQoUPavHmzNm/eHOql9eRDr/vP//xP3/333++LjIz0PfLII76amppQL6nf+d3vfueT1ONnzpw5oV5av3Gz/kjyvfzyy6FeWr/yj//4j9bv2ze+8Q3fxIkTfW63O9TLGhDGjRvne/bZZ0O9jH7lBz/4gW/o0KG+iIgIn9Pp9M2aNct37NixUC+rX/rNb37jy8zM9Nntdt9f/dVf+TZv3hzqJd0Uf0cHAAAYi8/oAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAAC96t1331VBQYGcTqdsNptee+21oM/h8/m0Zs0aPfTQQ7Lb7UpJSbG+RDQYfAUEAADoVR0dHXr44Yf1ox/9SH/zN39zW+d49tln5Xa7tWbNGmVlZamtrU2ffvpp0OfhLyMDAIA+Y7PZtGfPHs2cOdMa6+zs1L/8y79o586dunTpkjIzM/Xv//7vGj9+vCTpxIkTGjFihBobG5WRkXFH8/PWFQAAuKt+9KMf6b333lNlZaU++OADff/739eUKVP04YcfSpJ+85vf6IEHHtDevXuVlpamv/zLv9SPf/xjXbx4Mei5CDoAAOCu+eMf/6j//u//1quvvqrHH39cDz74oIqLi/XYY4/p5ZdfliR9/PHHOn36tF599VXt2LFD27ZtU319vf72b/826Pn4jA4AALhr3n//ffl8Pj300EN+4x6PR0OGDJEkdXd3y+PxaMeOHVbd1q1blZ2drZMnTwb1dhZBBwAA3DXd3d0KCwtTfX29wsLC/I7de++9kqShQ4cqPDzcLwwNGzZMknTmzBmCDgAA6J9Gjhyprq4utbS06PHHH79pzdixY/X555/rj3/8ox588EFJ0h/+8AdJ0v333x/UfNx1BQAAetWVK1f00UcfSfpzsFm3bp0mTJig+Ph4paam6oc//KHee+89rV27ViNHjtSnn36qt99+W1lZWfrrv/5rdXd369FHH9W9996r9evXq7u7WwsXLlRcXJzcbndQayHoAACAXvXOO+9owoQJPcbnzJmjbdu2yev1auXKldqxY4f+9Kc/aciQIcrNzdWKFSuUlZUlSfrkk0+0aNEiud1uxcTEaOrUqVq7dq3i4+ODWgtBBwAAGIvbywEAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAw1v8HKSpuqWLsR0wAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mapspam_df_raw.value_var.hist()" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [], "source": [ "res_crop = 6" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [], "source": [ "# def get_square(h3_cell: str, d: float = 10000):\n", " \n", "# h3_d = np.round(np.sqrt(h3.cell_area(h3_cell,unit=\"m^2\")),2)\n", "# frag = int(np.ceil(d/h3_d))\n", "\n", "# h3_ref = h3.cell_to_local_ij(origin=h3_cell,h=h3_cell)\n", "\n", "# i_range = range(-frag,frag+1)\n", "# j_range = range(-frag,frag+1)\n", "\n", "# coords = [(i,j) for i in i_range for j in j_range]\n", "# # coords\n", "# ref_cell = h3.cell_to_local_ij(origin=h3_cell,h=h3_cell)\n", "\n", "# return [h3.local_ij_to_cell(origin=h3_cell,i=ref_cell[0]+id[0],j=ref_cell[1]+id[1]) for id in coords]" ] }, { "cell_type": "code", "execution_count": 60, "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", "
yxvalue_var
013.375306-16.7923193647.500000
113.291973-16.79231911340.200195
214.291969-16.70898763.299999
314.208636-16.708987189.699997
414.125303-16.708987253.000000
\n", "
" ], "text/plain": [ " y x value_var\n", "0 13.375306 -16.792319 3647.500000\n", "1 13.291973 -16.792319 11340.200195\n", "2 14.291969 -16.708987 63.299999\n", "3 14.208636 -16.708987 189.699997\n", "4 14.125303 -16.708987 253.000000" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mapspam_df_raw.head()" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [], "source": [ "# get_square(\"8854a93225fffff\")" ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [], "source": [ "mapspam_df_raw[\"geom\"] = mapspam_df_raw.filter(regex=\"h3_id_\").apply(h3.cells_to_h3shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "mapspam_gdf = gpd.GeoDataFrame(mapspam_df\n", " ,geometry=\"geom\"\n", " ,crs=\"EPSG:4326\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "mapspam_gdf.iloc[0:100].explore(column=\"value\")\n", "mapspam_gdf.head()" ] }, { "cell_type": "code", "execution_count": 91, "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", "
yxvalue_vargeomh3_id
013.375306-16.7923193647.500000NaN8654a9337ffffff
113.291973-16.79231911340.200195NaN8654a931fffffff
214.291969-16.70898763.299999NaN8654ad387ffffff
314.208636-16.708987189.699997NaN8654ad2a7ffffff
414.125303-16.708987253.000000NaN8654ad2afffffff
\n", "
" ], "text/plain": [ " y x value_var geom h3_id\n", "0 13.375306 -16.792319 3647.500000 NaN 8654a9337ffffff\n", "1 13.291973 -16.792319 11340.200195 NaN 8654a931fffffff\n", "2 14.291969 -16.708987 63.299999 NaN 8654ad387ffffff\n", "3 14.208636 -16.708987 189.699997 NaN 8654ad2a7ffffff\n", "4 14.125303 -16.708987 253.000000 NaN 8654ad2afffffff" ] }, "execution_count": 91, "metadata": {}, "output_type": "execute_result" } ], "source": [ "try : mapspam_df_raw.drop(labels=[\"h3_id\"],axis={0},inplace=True)\n", "except : Warning(\"Nothing to remove, skipping\")\n", "mapspam_df_raw[\"h3_id\"] = mapspam_df_raw.apply(lambda row: h3.latlng_to_cell(row[\"y\"],row[\"x\"],res=res_crop),axis=1)\n", "mapspam_df_raw.head()" ] }, { "cell_type": "code", "execution_count": 90, "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", "
h3_idvalue_var
0865280477ffffff590.099976
1865280807ffffff1376.800049
286528080fffffff524.500000
38652808c7ffffff1311.199951
48652808d7ffffff1901.199951
\n", "
" ], "text/plain": [ " h3_id value_var\n", "0 865280477ffffff 590.099976\n", "1 865280807ffffff 1376.800049\n", "2 86528080fffffff 524.500000\n", "3 8652808c7ffffff 1311.199951\n", "4 8652808d7ffffff 1901.199951" ] }, "execution_count": 90, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mapspam_df = mapspam_df_raw[[\"h3_id\",\"value_var\"]].groupby(\"h3_id\").agg(\"sum\").reset_index()\n", "mapspam_df.head()" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [], "source": [ "compact_geo = sd.add_geom(mapspam_df)" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(53680, 3)" ] }, "execution_count": 66, "metadata": {}, "output_type": "execute_result" } ], "source": [ "compact_geo.shape" ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [], "source": [ "compact_geo[\"log_value_var\"] = np.round(np.log1p(compact_geo[\"value_var\"]))" ] }, { "cell_type": "code", "execution_count": 68, "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", "
h3_idvalue_vargeomlog_value_var
0865280477ffffff590.099976POLYGON ((38.27947 12.86167, 38.29737 12.83463...6.0
1865280807ffffff1376.800049POLYGON ((38.5914 13.4915, 38.57351 13.51876, ...7.0
286528080fffffff524.500000POLYGON ((38.5914 13.4915, 38.57484 13.46133, ...6.0
38652808c7ffffff1311.199951POLYGON ((38.57751 13.34657, 38.59536 13.31938...7.0
48652808d7ffffff1901.199951POLYGON ((38.57751 13.34657, 38.5431 13.34362,...8.0
\n", "
" ], "text/plain": [ " h3_id value_var \\\n", "0 865280477ffffff 590.099976 \n", "1 865280807ffffff 1376.800049 \n", "2 86528080fffffff 524.500000 \n", "3 8652808c7ffffff 1311.199951 \n", "4 8652808d7ffffff 1901.199951 \n", "\n", " geom log_value_var \n", "0 POLYGON ((38.27947 12.86167, 38.29737 12.83463... 6.0 \n", "1 POLYGON ((38.5914 13.4915, 38.57351 13.51876, ... 7.0 \n", "2 POLYGON ((38.5914 13.4915, 38.57484 13.46133, ... 6.0 \n", "3 POLYGON ((38.57751 13.34657, 38.59536 13.31938... 7.0 \n", "4 POLYGON ((38.57751 13.34657, 38.5431 13.34362,... 8.0 " ] }, "execution_count": 68, "metadata": {}, "output_type": "execute_result" } ], "source": [ "compact_geo.head()" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 69, "metadata": {}, "output_type": "execute_result" } ], "source": [ "compact_geo.loc[0:1000].explore(column=\"log_value_var\")" ] }, { "cell_type": "code", "execution_count": 92, "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", "
h3_idvalue_vargeomlog_value_var
0865280477ffffff590.099976POLYGON ((38.27947 12.86167, 38.29737 12.83463...6.0
1865280807ffffff1376.800049POLYGON ((38.5914 13.4915, 38.57351 13.51876, ...7.0
286528080fffffff524.500000POLYGON ((38.5914 13.4915, 38.57484 13.46133, ...6.0
38652808c7ffffff1311.199951POLYGON ((38.57751 13.34657, 38.59536 13.31938...7.0
48652808d7ffffff1901.199951POLYGON ((38.57751 13.34657, 38.5431 13.34362,...8.0
\n", "
" ], "text/plain": [ " h3_id value_var \\\n", "0 865280477ffffff 590.099976 \n", "1 865280807ffffff 1376.800049 \n", "2 86528080fffffff 524.500000 \n", "3 8652808c7ffffff 1311.199951 \n", "4 8652808d7ffffff 1901.199951 \n", "\n", " geom log_value_var \n", "0 POLYGON ((38.27947 12.86167, 38.29737 12.83463... 6.0 \n", "1 POLYGON ((38.5914 13.4915, 38.57351 13.51876, ... 7.0 \n", "2 POLYGON ((38.5914 13.4915, 38.57484 13.46133, ... 6.0 \n", "3 POLYGON ((38.57751 13.34657, 38.59536 13.31938... 7.0 \n", "4 POLYGON ((38.57751 13.34657, 38.5431 13.34362,... 8.0 " ] }, "execution_count": 92, "metadata": {}, "output_type": "execute_result" } ], "source": [ "compact_geo.head()" ] }, { "cell_type": "code", "execution_count": 93, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Skipping geom\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", "
h3_idchild_cellsvalue_varlog_value_var
0835280fffffffff[8452805ffffffff, 8452809ffffffff, 845280dffff...4.969470e+0410.81
1835282fffffffff[8452821ffffffff, 8452823ffffffff, 8452825ffff...1.189984e+0613.99
2835283fffffffff[8452831ffffffff, 8452835ffffffff, 8452837ffff...3.338324e+0512.72
3835284fffffffff[8452841ffffffff, 8452843ffffffff, 8452845ffff...1.562956e+0614.26
4835285fffffffff[8452853ffffffff, 845285dffffffff]1.291540e+049.47
\n", "
" ], "text/plain": [ " h3_id child_cells \\\n", "0 835280fffffffff [8452805ffffffff, 8452809ffffffff, 845280dffff... \n", "1 835282fffffffff [8452821ffffffff, 8452823ffffffff, 8452825ffff... \n", "2 835283fffffffff [8452831ffffffff, 8452835ffffffff, 8452837ffff... \n", "3 835284fffffffff [8452841ffffffff, 8452843ffffffff, 8452845ffff... \n", "4 835285fffffffff [8452853ffffffff, 845285dffffffff] \n", "\n", " value_var log_value_var \n", "0 4.969470e+04 10.81 \n", "1 1.189984e+06 13.99 \n", "2 3.338324e+05 12.72 \n", "3 1.562956e+06 14.26 \n", "4 1.291540e+04 9.47 " ] }, "execution_count": 93, "metadata": {}, "output_type": "execute_result" } ], "source": [ "compact_geo_downscaled = sd.change_scale(grid=compact_geo, level=-3)\n", "compact_geo_downscaled[\"log_value_var\"] = np.round(np.log1p(compact_geo_downscaled[\"value_var\"]),decimals=2)\n", "compact_geo_downscaled.head()" ] }, { "cell_type": "code", "execution_count": 94, "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", "
h3_idchild_cellsvalue_varlog_value_vargeom
0835280fffffffff[8452805ffffffff, 8452809ffffffff, 845280dffff...4.969470e+0410.81POLYGON ((37.97019 12.83405, 38.46072 12.47793...
1835282fffffffff[8452821ffffffff, 8452823ffffffff, 8452825ffff...1.189984e+0613.99POLYGON ((37.27428 11.9689, 37.76603 11.61796,...
2835283fffffffff[8452831ffffffff, 8452835ffffffff, 8452837ffff...3.338324e+0512.72POLYGON ((38.35619 11.87175, 38.83967 11.52114...
3835284fffffffff[8452841ffffffff, 8452843ffffffff, 8452845ffff...1.562956e+0614.26POLYGON ((38.07384 13.44699, 38.67299 13.69993...
4835285fffffffff[8452853ffffffff, 845285dffffffff]1.291540e+049.47POLYGON ((39.75747 13.58785, 39.8701 14.19935,...
\n", "
" ], "text/plain": [ " h3_id child_cells \\\n", "0 835280fffffffff [8452805ffffffff, 8452809ffffffff, 845280dffff... \n", "1 835282fffffffff [8452821ffffffff, 8452823ffffffff, 8452825ffff... \n", "2 835283fffffffff [8452831ffffffff, 8452835ffffffff, 8452837ffff... \n", "3 835284fffffffff [8452841ffffffff, 8452843ffffffff, 8452845ffff... \n", "4 835285fffffffff [8452853ffffffff, 845285dffffffff] \n", "\n", " value_var log_value_var \\\n", "0 4.969470e+04 10.81 \n", "1 1.189984e+06 13.99 \n", "2 3.338324e+05 12.72 \n", "3 1.562956e+06 14.26 \n", "4 1.291540e+04 9.47 \n", "\n", " geom \n", "0 POLYGON ((37.97019 12.83405, 38.46072 12.47793... \n", "1 POLYGON ((37.27428 11.9689, 37.76603 11.61796,... \n", "2 POLYGON ((38.35619 11.87175, 38.83967 11.52114... \n", "3 POLYGON ((38.07384 13.44699, 38.67299 13.69993... \n", "4 POLYGON ((39.75747 13.58785, 39.8701 14.19935,... " ] }, "execution_count": 94, "metadata": {}, "output_type": "execute_result" } ], "source": [ "max_geom = np.min([50_000,compact_geo_downscaled.shape[0]])\n", "\n", "compact_geo_downscaled_sub = sd.add_geom(compact_geo_downscaled.loc[0:max_geom].copy())\n", "\n", "compact_geo_downscaled_sub.reset_index(drop=True,inplace=True)\n", "\n", "compact_geo_downscaled_sub.head()" ] }, { "cell_type": "code", "execution_count": 73, "metadata": {}, "outputs": [], "source": [ "# compact_geo.to_file(\"outputs/compact_geo.geojson\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Color palette" ] }, { "cell_type": "code", "execution_count": 95, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgAAAABACAYAAABsv8+/AAAAFnRFWHRUaXRsZQB2aXJpZGlzIGNvbG9ybWFwrE0mCwAAABx0RVh0RGVzY3JpcHRpb24AdmlyaWRpcyBjb2xvcm1hcAtjl3IAAAAwdEVYdEF1dGhvcgBNYXRwbG90bGliIHYzLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZ2GZxVMAAAAydEVYdFNvZnR3YXJlAE1hdHBsb3RsaWIgdjMuOS4yLCBodHRwczovL21hdHBsb3RsaWIub3JnTz9adAAAAiJJREFUeJzt1kGSmzAURdEv2FqWkP0vJfQgMhQCGceV2Ttn4pL0EVQPum771X5vVVXVWv39XfrPeV193V5zS98f1sf5/fPjey73zu6/3Hv/uz2cz57f9vP68rxO9+/zre7nhvvG+et6vH92bw3PDfcsD+eX59+/53n96f3362/f87/vf5yr93Of72/fPV9P89tX3zGeH3OT8/07Zs+/32+TuXZZD8/VODf8W5uuH/b7vctlfuv7NazH8/t7ZnP7bz2cD3NL+/Ph3Hl+/efz83vWun/vuL++nquH9eu9w/uu6/vvOO49f/8xf77vOj+8b7Y/fMfse9ca/y7nv+d62a++X+f1vt+G/b7u+/u6TxzzS//tc2053QMABBEAABBIAABAIAEAAIEEAAAEEgAAEEgAAEAgAQAAgQQAAAQSAAAQSAAAQCABAACBBAAABBIAABBIAABAIAEAAIEEAAAEEgAAEEgAAEAgAQAAgQQAAAQSAAAQSAAAQCABAACBBAAABBIAABBIAABAIAEAAIEEAAAEEgAAEEgAAEAgAQAAgQQAAAQSAAAQSAAAQCABAACBBAAABBIAABBIAABAIAEAAIEEAAAEEgAAEEgAAEAgAQAAgQQAAAQSAAAQSAAAQCABAACBBAAABBIAABBIAABAIAEAAIEEAAAEEgAAEEgAAEAgAQAAgQQAAAQSAAAQSAAAQCABAACBBAAABBIAABBIAABAoB9ucImHxcKZtAAAAABJRU5ErkJggg==", "text/html": [ "
viridis
\"viridis
under
bad
over
" ], "text/plain": [ "" ] }, "execution_count": 95, "metadata": {}, "output_type": "execute_result" } ], "source": [ "viridis = sns.color_palette(\"viridis\", as_cmap=True)\n", "viridis" ] }, { "cell_type": "code", "execution_count": 96, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "6" ] }, "execution_count": 96, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bmlunge = pypal.load_cmap(\"Bmlunge\")\n", "bmlunge.N" ] }, { "cell_type": "code", "execution_count": 97, "metadata": {}, "outputs": [], "source": [ "def cmap(input, palette):\n", " m = np.max(input.to_list())\n", " l = palette.N\n", " print(\"Max input : {}, palette colors : {}\".format(m,l))\n", " return [[int(255*j) for j in palette(int(x/m*l))] for x in input] #" ] }, { "cell_type": "code", "execution_count": 98, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 10.81\n", "1 13.99\n", "2 12.72\n", "3 14.26\n", "4 9.47\n", "dtype: float32" ] }, "execution_count": 98, "metadata": {}, "output_type": "execute_result" } ], "source": [ "unique_vals = pd.Series(compact_geo_downscaled[\"log_value_var\"].unique())\n", "unique_vals.head()" ] }, { "cell_type": "code", "execution_count": 99, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Max input : 18.690000534057617, palette colors : 6\n" ] } ], "source": [ "cols = cmap(compact_geo_downscaled_sub[\"log_value_var\"],palette=bmlunge) #pd.Series([[255*j for j in x] for x in bmlunge(compact_geo_downscaled[\"log_value\"])])\n", "compact_geo_downscaled_sub[\"color\"] = cols\n" ] }, { "cell_type": "code", "execution_count": 100, "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", "
h3_idchild_cellsvalue_varlog_value_vargeomcolor
0835280fffffffff[8452805ffffffff, 8452809ffffffff, 845280dffff...4.969470e+0410.81POLYGON ((37.97019 12.83405, 38.46072 12.47793...[52, 170, 182, 255]
1835282fffffffff[8452821ffffffff, 8452823ffffffff, 8452825ffff...1.189984e+0613.99POLYGON ((37.27428 11.9689, 37.76603 11.61796,...[61, 118, 136, 255]
2835283fffffffff[8452831ffffffff, 8452835ffffffff, 8452837ffff...3.338324e+0512.72POLYGON ((38.35619 11.87175, 38.83967 11.52114...[61, 118, 136, 255]
3835284fffffffff[8452841ffffffff, 8452843ffffffff, 8452845ffff...1.562956e+0614.26POLYGON ((38.07384 13.44699, 38.67299 13.69993...[61, 118, 136, 255]
4835285fffffffff[8452853ffffffff, 845285dffffffff]1.291540e+049.47POLYGON ((39.75747 13.58785, 39.8701 14.19935,...[52, 170, 182, 255]
\n", "
" ], "text/plain": [ " h3_id child_cells \\\n", "0 835280fffffffff [8452805ffffffff, 8452809ffffffff, 845280dffff... \n", "1 835282fffffffff [8452821ffffffff, 8452823ffffffff, 8452825ffff... \n", "2 835283fffffffff [8452831ffffffff, 8452835ffffffff, 8452837ffff... \n", "3 835284fffffffff [8452841ffffffff, 8452843ffffffff, 8452845ffff... \n", "4 835285fffffffff [8452853ffffffff, 845285dffffffff] \n", "\n", " value_var log_value_var \\\n", "0 4.969470e+04 10.81 \n", "1 1.189984e+06 13.99 \n", "2 3.338324e+05 12.72 \n", "3 1.562956e+06 14.26 \n", "4 1.291540e+04 9.47 \n", "\n", " geom color \n", "0 POLYGON ((37.97019 12.83405, 38.46072 12.47793... [52, 170, 182, 255] \n", "1 POLYGON ((37.27428 11.9689, 37.76603 11.61796,... [61, 118, 136, 255] \n", "2 POLYGON ((38.35619 11.87175, 38.83967 11.52114... [61, 118, 136, 255] \n", "3 POLYGON ((38.07384 13.44699, 38.67299 13.69993... [61, 118, 136, 255] \n", "4 POLYGON ((39.75747 13.58785, 39.8701 14.19935,... [52, 170, 182, 255] " ] }, "execution_count": 100, "metadata": {}, "output_type": "execute_result" } ], "source": [ "compact_geo_downscaled_sub.head()" ] }, { "cell_type": "code", "execution_count": 101, "metadata": {}, "outputs": [], "source": [ "# viewport = pdk.data_utils.compute_view(points=compact_geo_downscaled[['x', 'y']], view_proportion=0.9)\n", "viewport = pdk.ViewState(longitude=0,latitude=0,zoom=3)\n", "# auto_zoom_map = pdk.Deck(layers=[], initial_view_state=viewport)\n", "# auto_zoom_map.show()" ] }, { "cell_type": "code", "execution_count": 88, "metadata": {}, "outputs": [], "source": [ "\n", "# scaled_layer = pdk.Layer(\n", "# \"GeoJsonLayer\",\n", "# compact_geo_downscaled_sub,\n", "# pickable=True,\n", "# extruded=False,\n", "# filled=True,\n", "# stroked=True,\n", "# opacity=.6,\n", "# get_geometry = \"geom\",\n", "# get_fill_color = \"color\", #\"[255*log_value/20,100,120]\",\n", "# get_line_color= [255, 255, 255, 0],\n", "# line_width_min_pixels=1,\n", "# )\n", "\n", "h3_layer = pdk.Layer(\n", " \"H3HexagonLayer\",\n", " compact_geo_downscaled_sub,\n", " pickable=True,\n", " stroked=True,\n", " filled=True,\n", " opacity=.5,\n", " extruded=False,\n", " get_hexagon=\"h3_id\",\n", " get_fill_color= \"color\",\n", " get_line_color=[0, 0, 0, 0],\n", " line_width_min_pixels=1,\n", ")\n", "\n", "# raster_centr_layer = pdk.Layer(\n", "# \"ScatterplotLayer\",\n", "# mapspam_df_raw,\n", "# pickable=True,\n", "# opacity=0.3,\n", "# stroked=True,\n", "# filled=True,\n", "# radius_scale=1,\n", "# radius_min_pixels=4,\n", "# radius_max_pixels=5,\n", "# line_width_min_pixels=1,\n", "# get_position=[\"x\",\"y\"],\n", "# get_radius=5,\n", "# get_fill_color=[255, 140, 0],\n", "# get_line_color=[0, 0, 0, 0],\n", "# )\n", "\n", "r = pdk.Deck(layers=[h3_layer]\n", " ,initial_view_state=viewport,tooltip=True)\n", "\n", "# h3_layer\n", "\n", "# Create an HTML header to display the year\n", "# display_el = ipywidgets.HTML('

Cropland in Sub-Saharan africa

')\n", "# display(display_el)\n", "# Show the current visualization\n", "# r.show()\n" ] }, { "cell_type": "code", "execution_count": 102, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 102, "metadata": {}, "output_type": "execute_result" } ], "source": [ "r.to_html(\"../../deck_maps/cropland_downscaled.html\",iframe_height=800)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Algorightm ameliorations" ] }, { "cell_type": "code", "execution_count": 103, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(53680, 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", "
yxvalue_vargeomh3_id
013.375306-16.7923193647.500000NaN8654a9337ffffff
113.291973-16.79231911340.200195NaN8654a931fffffff
214.291969-16.70898763.299999NaN8654ad387ffffff
\n", "
" ], "text/plain": [ " y x value_var geom h3_id\n", "0 13.375306 -16.792319 3647.500000 NaN 8654a9337ffffff\n", "1 13.291973 -16.792319 11340.200195 NaN 8654a931fffffff\n", "2 14.291969 -16.708987 63.299999 NaN 8654ad387ffffff" ] }, "execution_count": 103, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# mapspam_harv_area = glob.glob(\"mapspam/spam2017v2r1_ssa_harv_area/*\")\n", "# from re import search\n", "# mapspam_harv_area_sel = [x for x in mapspam_harv_area if search(string=x,pattern=\"TEAS\")]\n", "# mapspam_harv_area_sel\n", "# mapsmam_out_path = \"mapspam_teas_ha.geojson\"\n", "\n", "teas_grid = mapspam_df_raw.copy()\n", "# sn.rast_to_centroid(out_path=mapsmam_out_path,tiff_paths=mapspam_harv_area_sel[1:2]).reset_index(inplace=False,drop=True)\n", "\n", "print(teas_grid.shape)\n", "\n", "teas_grid.head(3)\n" ] }, { "cell_type": "code", "execution_count": 104, "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", "
yxvalue_vargeomh3_idgeometry
013.375306-16.7923193647.500000NaN8654a9337ffffffPOINT (-16.79232 13.37531)
113.291973-16.79231911340.200195NaN8654a931fffffffPOINT (-16.79232 13.29197)
214.291969-16.70898763.299999NaN8654ad387ffffffPOINT (-16.70899 14.29197)
314.208636-16.708987189.699997NaN8654ad2a7ffffffPOINT (-16.70899 14.20864)
414.125303-16.708987253.000000NaN8654ad2afffffffPOINT (-16.70899 14.1253)
\n", "
" ], "text/plain": [ " y x value_var geom h3_id \\\n", "0 13.375306 -16.792319 3647.500000 NaN 8654a9337ffffff \n", "1 13.291973 -16.792319 11340.200195 NaN 8654a931fffffff \n", "2 14.291969 -16.708987 63.299999 NaN 8654ad387ffffff \n", "3 14.208636 -16.708987 189.699997 NaN 8654ad2a7ffffff \n", "4 14.125303 -16.708987 253.000000 NaN 8654ad2afffffff \n", "\n", " geometry \n", "0 POINT (-16.79232 13.37531) \n", "1 POINT (-16.79232 13.29197) \n", "2 POINT (-16.70899 14.29197) \n", "3 POINT (-16.70899 14.20864) \n", "4 POINT (-16.70899 14.1253) " ] }, "execution_count": 104, "metadata": {}, "output_type": "execute_result" } ], "source": [ "teas_grid = gpd.GeoDataFrame(teas_grid, geometry=teas_grid.apply(lambda row: shp.Point(row[\"x\"],row[\"y\"]),axis=1),crs=4326)\n", "teas_grid.head()" ] }, { "cell_type": "code", "execution_count": 105, "metadata": {}, "outputs": [], "source": [ "\n", "centroid = [teas_grid.geometry.loc[0:100].x.median(),teas_grid.geometry.loc[0:100].y.median()]\n", "\n", "centroid = gpd.GeoSeries(shp.Point(centroid),crs=teas_grid.crs).to_crs(\"EPSG:4326\")\n", "\n", "lat_center = float(centroid.geometry.y.iloc[0])\n", "lon_center = float(centroid.geometry.x.iloc[0])" ] }, { "cell_type": "code", "execution_count": 106, "metadata": {}, "outputs": [], "source": [ "view_state = pdk.ViewState(latitude=lat_center, longitude=lon_center, zoom=6, bearing=0, pitch=0)\n", "teas_grid_4326 = teas_grid.to_crs(\"EPSG:4326\")\n", "teas_grid_4326[[\"x\",\"y\"]] = teas_grid_4326.get_coordinates()\n", "teas_grid_4326.reset_index(inplace=True,drop=True)\n" ] }, { "cell_type": "code", "execution_count": 107, "metadata": {}, "outputs": [], "source": [ "\n", "layer_total = pdk.Layer(\n", " \"ScatterplotLayer\",\n", " teas_grid_4326.loc[0:1000],\n", " pickable=True,\n", " opacity=0.5,\n", " stroked=True,\n", " filled=True,\n", " radius_scale=1,\n", " radius_min_pixels=3,\n", " radius_max_pixels=5,\n", " line_width_max_pixels=0,\n", " get_position= [\"x\",\"y\"],\n", " get_radius=2,\n", " get_fill_color= [40, 2, 61],\n", " get_line_color=[0, 0, 0]\n", ")\n", "\n", "# Create mapdeck object\n", "deck_map = pdk.Deck(layers=[layer_total\n", " ]\n", " ,initial_view_state=view_state\n", " ,tooltip= {\"text\": \"Value : {band}\"}\n", " ,height=800\n", " ,map_style=\"road\"\n", " )\n" ] }, { "cell_type": "code", "execution_count": 109, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 109, "metadata": {}, "output_type": "execute_result" } ], "source": [ "deck_map.to_html(\"../../deck_maps/deck_mapspam_teas_ha_grid.html\")\n" ] }, { "cell_type": "code", "execution_count": 110, "metadata": {}, "outputs": [], "source": [ "h3_resolution = 8" ] }, { "cell_type": "code", "execution_count": 111, "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", "
yxvalue_vargeomh3_idgeometry
013.375306-16.7923193647.500000NaN8854a93225fffffPOINT (-16.79232 13.37531)
113.291973-16.79231911340.200195NaN8854a931d3fffffPOINT (-16.79232 13.29197)
214.291969-16.70898763.299999NaN8854ad3819fffffPOINT (-16.70899 14.29197)
314.208636-16.708987189.699997NaN8854ad2a51fffffPOINT (-16.70899 14.20864)
414.125303-16.708987253.000000NaN8854ad2ab7fffffPOINT (-16.70899 14.1253)
\n", "
" ], "text/plain": [ " y x value_var geom h3_id \\\n", "0 13.375306 -16.792319 3647.500000 NaN 8854a93225fffff \n", "1 13.291973 -16.792319 11340.200195 NaN 8854a931d3fffff \n", "2 14.291969 -16.708987 63.299999 NaN 8854ad3819fffff \n", "3 14.208636 -16.708987 189.699997 NaN 8854ad2a51fffff \n", "4 14.125303 -16.708987 253.000000 NaN 8854ad2ab7fffff \n", "\n", " geometry \n", "0 POINT (-16.79232 13.37531) \n", "1 POINT (-16.79232 13.29197) \n", "2 POINT (-16.70899 14.29197) \n", "3 POINT (-16.70899 14.20864) \n", "4 POINT (-16.70899 14.1253) " ] }, "execution_count": 111, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "teas_grid_4326 = dt.df_project_on_grid(teas_grid_4326,res=h3_resolution)\n", "teas_grid_4326.head()\n" ] }, { "cell_type": "code", "execution_count": 133, "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" ] }, { "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", "
yxvalue_vargeomh3_idgeometryh3_gridded
013.375306-16.7923193647.500000NaN8854a93225fffffPOINT (-16.79232 13.37531)[8854a93201fffff, 8854a9320bfffff, 8854a93215f...
113.291973-16.79231911340.200195NaN8854a931d3fffffPOINT (-16.79232 13.29197)[8854a9308dfffff, 8854a93089fffff, 8854a93083f...
214.291969-16.70898763.299999NaN8854ad3819fffffPOINT (-16.70899 14.29197)[8854ad3a85fffff, 8854ad3a81fffff, 8854ad3ab9f...
314.208636-16.708987189.699997NaN8854ad2a51fffffPOINT (-16.70899 14.20864)[8854ad216bfffff, 8854ad2147fffff, 8854ad2101f...
414.125303-16.708987253.000000NaN8854ad2ab7fffffPOINT (-16.70899 14.1253)[8854ad2a93fffff, 8854ad2f4dfffff, 8854ad2f47f...
\n", "
" ], "text/plain": [ " y x value_var geom h3_id \\\n", "0 13.375306 -16.792319 3647.500000 NaN 8854a93225fffff \n", "1 13.291973 -16.792319 11340.200195 NaN 8854a931d3fffff \n", "2 14.291969 -16.708987 63.299999 NaN 8854ad3819fffff \n", "3 14.208636 -16.708987 189.699997 NaN 8854ad2a51fffff \n", "4 14.125303 -16.708987 253.000000 NaN 8854ad2ab7fffff \n", "\n", " geometry \\\n", "0 POINT (-16.79232 13.37531) \n", "1 POINT (-16.79232 13.29197) \n", "2 POINT (-16.70899 14.29197) \n", "3 POINT (-16.70899 14.20864) \n", "4 POINT (-16.70899 14.1253) \n", "\n", " h3_gridded \n", "0 [8854a93201fffff, 8854a9320bfffff, 8854a93215f... \n", "1 [8854a9308dfffff, 8854a93089fffff, 8854a93083f... \n", "2 [8854ad3a85fffff, 8854ad3a81fffff, 8854ad3ab9f... \n", "3 [8854ad216bfffff, 8854ad2147fffff, 8854ad2101f... \n", "4 [8854ad2a93fffff, 8854ad2f4dfffff, 8854ad2f47f... " ] }, "execution_count": 133, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scalenav.data import rast_to_h3_map\n", "\n", "rast_to_h3 = rast_to_h3_map(x=teas_grid_4326.x[0],y=teas_grid_4326.y[0],ref=\"m\")\n", "\n", "teas_grid_4326[\"h3_gridded\"] = pd.Series(teas_grid_4326[\"h3_id\"].apply(lambda cell: dt.centre_cell_to_square(h3_cell=cell,neighbs=rast_to_h3[\"10000\"][\"nn\"])).tolist())\n", "teas_grid_4326.head()\n" ] }, { "cell_type": "code", "execution_count": 134, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(7461520, 7)\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", "
yxvalue_vargeomh3_idgeometryh3_gridded
013.375306-16.7923193647.5NaN8854a93225fffffPOINT (-16.79232 13.37531)8854a93201fffff
013.375306-16.7923193647.5NaN8854a93225fffffPOINT (-16.79232 13.37531)8854a9320bfffff
013.375306-16.7923193647.5NaN8854a93225fffffPOINT (-16.79232 13.37531)8854a93215fffff
\n", "
" ], "text/plain": [ " y x value_var geom h3_id \\\n", "0 13.375306 -16.792319 3647.5 NaN 8854a93225fffff \n", "0 13.375306 -16.792319 3647.5 NaN 8854a93225fffff \n", "0 13.375306 -16.792319 3647.5 NaN 8854a93225fffff \n", "\n", " geometry h3_gridded \n", "0 POINT (-16.79232 13.37531) 8854a93201fffff \n", "0 POINT (-16.79232 13.37531) 8854a9320bfffff \n", "0 POINT (-16.79232 13.37531) 8854a93215fffff " ] }, "execution_count": 134, "metadata": {}, "output_type": "execute_result" } ], "source": [ "teas_h3_gridded = teas_grid_4326.explode(\"h3_gridded\")\n", "# teas_h3_gridded[\"value_var\"] = teas_h3_gridded[\"value_var\"]/teas_h3_gridded[\"coeff\"]\n", "print(teas_h3_gridded.shape)\n", "teas_h3_gridded.head(3)" ] }, { "cell_type": "code", "execution_count": 135, "metadata": {}, "outputs": [], "source": [ "centroid = [teas_grid_4326.iloc[0:100].x.mean(),teas_grid_4326.iloc[0:100].y.mean()]\n", "centroid = gpd.GeoSeries(shp.Point(centroid),crs=teas_grid_4326.crs).to_crs(\"EPSG:4326\")\n", "lat_center = float(centroid.geometry.y.iloc[0])\n", "lon_center = float(centroid.geometry.x.iloc[0])" ] }, { "cell_type": "code", "execution_count": 136, "metadata": {}, "outputs": [], "source": [ "view_state = pdk.ViewState(latitude=lat_center, longitude=lon_center, zoom=6, bearing=0, pitch=0)\n", "# subset_data = grid_4326[(grid_4326.band>0) and (grid_4326.band < max_val)]" ] }, { "cell_type": "code", "execution_count": 137, "metadata": {}, "outputs": [], "source": [ "layer_total = pdk.Layer(\n", " \"ScatterplotLayer\",\n", " teas_grid_4326.loc[0:100],\n", " pickable=True,\n", " opacity=1,\n", " stroked=True,\n", " filled=True,\n", " radius_scale=5,\n", " radius_min_pixels=3,\n", " radius_max_pixels=5,\n", " line_width_max_pixels=0,\n", " get_position= [\"x\",\"y\"],\n", " get_radius=1,\n", " get_fill_color= [200, 0, 150, 255],\n", " get_line_color=[0, 0, 0]\n", ")\n", "\n", "\n", "h3_layer = pdk.Layer(\n", " \"H3HexagonLayer\",\n", " teas_h3_gridded.loc[0:100],\n", " pickable=True,\n", " stroked=True,\n", " filled=True,\n", " opacity=.1,\n", " extruded=False,\n", " get_hexagon=\"h3_gridded\",\n", " get_fill_color= [0,0,0,255],\n", " get_line_color=[0, 0, 0, 0],\n", " line_width_min_pixels=1,\n", ")\n", "\n" ] }, { "cell_type": "code", "execution_count": 138, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 138, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "# Create mapdeck object\n", "deck_map = pdk.Deck(layers=[layer_total\n", " ,h3_layer]\n", " ,initial_view_state=view_state\n", " ,tooltip= {\"text\": \"Value : {value_var}\"}\n", " ,height=800\n", " ,map_style=\"road\"\n", " )\n", "\n", "deck_map.to_html(\"../../deck_maps/deck_mapspam_teas_ha_grid.html\")" ] } ], "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 }