Proxying layers

Among the many layers of data available, some have actually high granularity and will be used as a renormalising function to distribute values that would otherwise be uniform in space. The idea is to link certain economic activity related variables to some infrastructure and geographic layers, and assuming some form of homogeneity, distribute the economic activity proposrtionally to the the distribution, say, of infrastrucutre in a location. This notebook will lay some basis for this.

#Base
import itertools as iter
import os
import re
import math
import numpy as np
import random as rd
import glob
import dask

# data eng
import h3
import shapely as shp
import geopandas as gpd
import pandas as pd
import rioxarray as rx
import osmnx as ox
import pyarrow
import ibis as ib
from ibis import _
ib.options.interactive = True
import ibis.selectors as s
# import polars as pl

# module
import scalenav.scale_nav as sd
import scalenav.data as dt
from scalenav.plotting import cmap

### Visualisation related
from IPython.display import display
import ipywidgets
# import contextily as cx
import pydeck as pdk
from matplotlib import pyplot as plt
import seaborn as sns
import pypalettes as pypal

Relevant docs

ddb spatial extension web: https://duckdb.org/docs/extensions/spatial.html

ddb spatial extension github functions ref: https://github.com/duckdb/duckdb_spatial/blob/main/docs/functions.md

H3 + duckdb extension : https://github.com/isaacbrodsky/h3-duckdb

ibis + geospatial :

general : https://ibis-project.org/reference/expression-tables#attributes

blogpost : https://ibis-project.org/posts/ibis-duckdb-geospatial/#functions-supported-and-next-steps

functions ref : https://ibis-project.org/reference/expression-geospatial

# # plotting support function

# def cmap(input, palette):
#     m = np.max(input.to_list())
#     l = palette.N
#     print("Max input : {}, palette colors : {}".format(m,l))
#     return [[int(255*j) for j in palette(int(x/m*l))] for x in input] #

grid_param = 100

ghsl_file = f"/Users/cenv1069/Documents/data/datasets/JRC/S_{grid_param}m/**"
# "/Users/cenv1069/Documents/agriculture/ghsl/non_res/**"
# "/Users/cenv1069/Documents/data/datasets/JRC/S_10m/**"

tiffs = [x for x in glob.glob(ghsl_file,recursive=True) if re.search(pattern=".tif$",string=x)]
print(len(tiffs))
tiffs

4

['/Users/cenv1069/Documents/data/datasets/JRC/S_100m/GHS_BUILT_S_NRES_E2015_GLOBE_R2023A_54009_100_V1_0_R8_C18/GHS_BUILT_S_NRES_E2015_GLOBE_R2023A_54009_100_V1_0_R8_C18.tif',
 '/Users/cenv1069/Documents/data/datasets/JRC/S_100m/GHS_BUILT_S_NRES_E2015_GLOBE_R2023A_54009_100_V1_0_R8_C20/GHS_BUILT_S_NRES_E2015_GLOBE_R2023A_54009_100_V1_0_R8_C20.tif',
 '/Users/cenv1069/Documents/data/datasets/JRC/S_100m/GHS_BUILT_S_NRES_E2015_GLOBE_R2023A_54009_100_V1_0_R8_C19/GHS_BUILT_S_NRES_E2015_GLOBE_R2023A_54009_100_V1_0_R8_C19.tif',
 '/Users/cenv1069/Documents/data/datasets/JRC/S_100m/GHS_BUILT_S_NRES_E2015_GLOBE_R2023A_54009_100_V1_0_R9_C19/GHS_BUILT_S_NRES_E2015_GLOBE_R2023A_54009_100_V1_0_R9_C19.tif']

out_folder = "outputs"

out_filename = f'S_{grid_param}m_raw_grid.parquet'

out_file = f"../../{out_folder}/{out_filename}"

if os.path.exists(out_file):
    print("Reading existing file")
    raw_grid = pd.read_parquet(out_file)
else :
    print("Converting rasters to centroids.")
    raw_grid = dt.rast_to_centroid(out_file,in_paths=tiffs)
    raw_grid.to_parquet(out_file)

Reading existing file

print(raw_grid.shape)
raw_grid.head()

(67835, 3)
band lon lat
142378 113 -803150.0 1998550.0
152379 133 -803050.0 1998450.0
322425 112 -798450.0 1996750.0
451074 51 -933550.0 1995450.0
461073 80 -933650.0 1995350.0 
raw_grid = gpd.GeoDataFrame(raw_grid,geometry=gpd.GeoSeries(raw_grid.apply(lambda row: shp.Point(row.loc["lon"],row.loc["lat"]),axis=1),crs="ESRI:54009")).to_crs(epsg=4326)
raw_grid[["x","y"]] = raw_grid.get_coordinates()

raw_grid.drop(columns=["lon","lat"],inplace=True)
raw_grid.head()
band geometry x y
142378 113 POINT (-8.21791 16.24726) -8.217905 16.247257
152379 133 POINT (-8.21686 16.24644) -8.216861 16.246436
322425 112 POINT (-8.16943 16.23247) -8.169435 16.232471
451074 51 POINT (-9.55141 16.22179) -9.551406 16.221792
461073 80 POINT (-9.5524 16.22097) -9.552404 16.220971 
# raw_grid.crs
# import pyarrow.parquet
# pyarrow.parquet.read_metadata(out_file)

raw_grid.band.hist()

<Axes: >
../_images/759ff3c3c3b6b32491b87260954b2b6b0139de935c0dfcd8e0f46b058cde103d.png 
raw_grid["log_band"] = np.log1p(raw_grid.band)

# raw_grid.sample(n=20_000).explore()

ghsl_df = raw_grid[["band","x","y"]].copy()

ghsl_df
band x y
142378 113 -8.217905 16.247257
152379 133 -8.216861 16.246436
322425 112 -8.169435 16.232471
451074 51 -9.551406 16.221792
461073 80 -9.552404 16.220971
... ... ... ...
98557016 1560 6.591825 0.116867
98567015 468 6.590827 0.116058
98567016 1471 6.591825 0.116058
99116953 12 6.528960 0.071576
99416960 40 6.535942 0.047313

67835 rows × 3 columns

Projecting the old way

# ghsl_gdf = raw_grid.copy() #gpd.read_file("/Users/cenv1069/Documents/agriculture/data/ghsl_clean_data.geojson")

# ghsl_gdf.reset_index(inplace=True,drop=True)
# ghsl_gdf.loc[0:10000].explore("band")
# dt.rast_to_h3_map(x = ghsl_gdf.x.mean(),y=ghsl_gdf.y.mean())
# # h3_resolution = 12
# grid_param=100
# h3_resolution = dt.rast_to_h3["100"]["h3_res"]
# ghsl_grid = dt.df_project_on_grid(ghsl_gdf,res=h3_resolution)
# ghsl_grid["h3_gridded"] = pd.Series(ghsl_grid.apply(lambda row: dt.centre_cell_to_square(row["h3_id"],grid_param=grid_param),axis=1).tolist())
# ghsl_grid.drop(columns=["geometry"],inplace=True)
# ghsl_grid.dropna(subset=["h3_gridded"],inplace=True)
# ghsl_gridded = ghsl_grid.explode("h3_gridded")
# print(ghsl_gridded.shape)
# ghsl_gridded.head()
# ghsl_gridded = ghsl_gridded[~ghsl_gridded.duplicated(subset=["h3_gridded"])]
# print(ghsl_gridded.shape)
# ghsl_gridded.reset_index(inplace=True,drop=True)
# ghsl_gridded.head()
# ghsl_gridded.set_index("h3_id",inplace=True,drop=True)
# ghsl_gridded = ghsl_gridded.join(ghsl_gridded.value_counts("h3_id"),how="left",sort=False,validate="m:1")
# ghsl_gridded["band"] /= ghsl_gridded["count"]
# ghsl_gridded.reset_index(inplace=True,drop=False)
# ghsl_gridded.head()
# ghsl_gridded.dropna(subset=["h3_id","h3_gridded"],inplace=True)
# ghsl_gridded.band.hist()
# import pypalettes as pypal

# bmlunge = pypal.load_cmap("Bmlunge")
# bmlunge.N
# # ghsl_gridded["log_band"] = np.log1p(ghsl_gridded.band)
# ghsl_gridded["cols"] = cmap(ghsl_gridded.log_band,bmlunge)
# # centroid = [ghsl_gridded.x.mean(),ghsl_gridded.y.mean()]
# # centroid = gpd.GeoSeries(shp.Point(centroid),crs=ghsl_gridded.crs).to_crs("EPSG:4326")
# centroid = gpd.GeoSeries(shp.Point(0,0),crs=4326)
# lat_center = float(centroid.geometry.y.iloc[0])
# lon_center = float(centroid.geometry.x.iloc[0])
# # view_state = pdk.ViewState(latitude=lat_center, longitude=lon_center, zoom=4, bearing=0, pitch=0)
# view_state = pdk.data_utils.compute_view(points=ghsl_df[["x","y"]])
# # subset_data = grid_4326[(grid_4326.band>0) and (grid_4326.band < max_val)]
# print(ghsl_gridded.shape)
# ghsl_gridded.head()
# deck_data_size = np.min([500_000,ghsl_gridded.shape[0]])

# h3_layer_ghsl = pdk.Layer(
#     "H3HexagonLayer",
#     ghsl_gridded.loc[0:deck_data_size], # .sample(n=200_000)
#     pickable=True,
#     stroked=True,
#     filled=True,
#     opacity=.5,
#     extruded=False,
#     get_hexagon= "h3_gridded",
#     get_fill_color = [200,200,100,255], # "cols",
#     get_line_color= [0, 0, 0, 100],
#     line_width_min_pixels=0,
#     line_width_max_pixels=1,
# )

# # Create mapdeck object
# deck_map_ghsl = pdk.Deck(layers=[h3_layer_ghsl]
#              ,initial_view_state=view_state
#              ,tooltip= {"text": "Value :  {band}"}
#              ,height=800
#              ,map_style="road"
#              )

# deck_map_ghsl.to_html("../../deck_maps/deck_ghsl_nres_grid.html")

H3 + duckdb scale up

conn = ib.connect("duckdb://")

res = conn.raw_sql("""INSTALL h3 FROM community;
                   LOAD h3;""").df()

res
Success 
# loading the spatial extension.
res = conn.raw_sql("install spatial; load spatial").df()
res
Success 
ghsl_df.head()
band x y
142378 113 -8.217905 16.247257
152379 133 -8.216861 16.246436
322425 112 -8.169435 16.232471
451074 51 -9.551406 16.221792
461073 80 -9.552404 16.220971

H3 and grid parameters

dose_wdi_geo_file = "/Users/cenv1069/Documents/data/datasets/local_data/dose-wdi/dose_wdi_geo.parquet"

dwg = conn.read_parquet(dose_wdi_geo_file,table_name="dwg")

dwg.head()

┏━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━━┓
┃ gid_0  ┃ country ┃ gid_1   ┃ grp_usd_2015 ┃ services_usd_2015 ┃ manufacturing_usd_2015 ┃ agriculture_usd_2015 ┃ centr                                                           ┃ geometry                                                                                                                                                                                                                                    ┃ color             ┃ radius    ┃ x         ┃ y         ┃
┡━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━━┩
│ string │ string  │ string  │ float64      │ float64           │ float64                │ float64              │ binary                                                          │ binary                                                                                                                                                                                                                                      │ array<int64>      │ float64   │ float64   │ float64   │
├────────┼─────────┼─────────┼──────────────┼───────────────────┼────────────────────────┼──────────────────────┼─────────────────────────────────────────────────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┼───────────────────┼───────────┼───────────┼───────────┤
│ ALB    │ Albania │ ALB.1_1 │ 4.341915e+08 │      2.009631e+08 │           9.449606e+07 │         8.588401e+07 │ b'\x01\x01\x00\x00\x00@\xa0\x88\x02=\x174@\x9d\xa3\xc5\xb3rPD@' │ b'\x01\x03\x00\x00\x00\x01\x00\x00\x00\xbe\x02\x00\x00\x07h\x11@ a4@\xd3\xa9\x02 \xb47D@\xcd\x15\x10 d`4@\xc4G\xfe?\xd27D@\xa0*\xf0_j_4@\xd8\x98\xf8\x1f\xf07D@$\xb3\xfe\xbf\x7f^4@.:\xf6?\x028D@\x88\xce\xf8'+11165                        │ [238, 97, ... +2] │ 19.888996 │ 20.090775 │ 40.628500 │
│ ALB    │ Albania │ ALB.2_1 │ 3.719427e+08 │      1.721516e+08 │           8.094841e+07 │         7.357106e+07 │ b'\x01\x01\x00\x00\x00\xe6M\xcan\xe1?4@\x87\xf1\xed\xcdg\xcbD@' │ b'\x01\x03\x00\x00\x00\x01\x00\x00\x00\xbb\x02\x00\x00\xa0\x16\x07\xc0\rl4@\xf5\x0e\xf9_\x15\xa4D@\xdfC\x15@)l4@u\xe3\xfe\x9f\x1a\xa4D@\x1f\x99\xf5\xff\xfdk4@:\xa9\xf5\x9f1\xa4D@\xb5O\t\x80\xd8k4@\xab\xc4\x02\x80_\xa4D@\x9fY\x12'+11117 │ [238, 97, ... +2] │ 19.734250 │ 20.249534 │ 41.589105 │
│ ALB    │ Albania │ ALB.3_1 │ 1.113524e+09 │      5.153885e+08 │           2.423439e+08 │         2.202575e+08 │ b'\x01\x01\x00\x00\x00\xb5?\xc3\x1e`\xa13@\x14\x15\xcd2B\xb9D@' │ b'\x01\x03\x00\x00\x00\x01\x00\x00\x00\x8f\x0c\x00\x00\xb3@\xfd_\\\xa33@9\xad\x04 \xa8\xabD@\x85\x1a\xf0\xbf\xb6\xa23@\x13U\t`\x14\xabD@~r\x14 \n\xa33@+\xf0\xfc?\xda\xaaD@\xa8O\xf6\xbf\xf1\xa33@\xd1\x9a\xf6\xff7\xaaD@\xc0@\x10'+51373   │ [238, 97, ... +2] │ 20.830796 │ 19.630373 │ 41.447333 │
│ ALB    │ Albania │ ALB.4_1 │ 7.954835e+08 │      3.681851e+08 │           1.731265e+08 │         1.573483e+08 │ b'\x01\x01\x00\x00\x00\x17\xb66{\xbb/4@m\xda[\xa0\x1f\x85D@'    │ b'\x01\x03\x00\x00\x00\x01\x00\x00\x00\xc5\x02\x00\x00\xe3\xe0\xf3?\xd0o4@\xe3\xa1\xf7\xdf\xadlD@ti\x15`\xcco4@\xb0\xb3\xf8\x7f\x9blD@A\x93\x0e \x8dp4@\x1f\x03\x05 flD@\xd8E\x13 \xf1p4@\xd5\x0f\x03\xc0\xa5kD@\xa6,\x01'+11277            │ [238, 97, ... +2] │ 20.494461 │ 20.186454 │ 41.040028 │
│ ALB    │ Albania │ ALB.5_1 │ 1.345159e+09 │      6.225993e+08 │           2.927561e+08 │         2.660753e+08 │ b"\x01\x01\x00\x00\x00\xee\x16g\xe7\xda\x9e3@Y\x90\xfc/'cD@"    │ b'\x01\x03\x00\x00\x00\x01\x00\x00\x00,\n\x00\x00x\x06\xf4_s\xbe3@\xfc\x0c\xfb\x7f\xcd9D@+\xb1\x00\xe0\xb7\xbe3@\x8c\xc1\xfd_y:D@\xf2.\x0f\x80U\xbe3@\x159\x02\xc0\xca:D@\xdf\xed\xee\xbf\x9d\xbd3@\xf6\xb7\x04\xe0\x1f;D@,B\x14'+41597     │ [238, 97, ... +2] │ 21.019778 │ 19.620528 │ 40.774633 │
└────────┴─────────┴─────────┴──────────────┴───────────────────┴────────────────────────┴──────────────────────┴─────────────────────────────────────────────────────────────────┴─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴───────────────────┴───────────┴───────────┴───────────┘

# conn.drop_view("ghsl_h3_",force=True)
# conn.drop_view("ghsl_h3",force=True)

print("Assuming raster grid resolution : ",grid_param)

h3_res = dt.rast_to_h3[str(grid_param)]["h3_res"]

print("Using base H3 resolution : ", h3_res)

res = conn.raw_sql(f"""create or replace view ghsl_h3 as
               (select    
                    band,
                    h3_h3_to_string(h3_latlng_to_cell(y,x,{h3_res})) as h3_id,
                    x,
                    y
                   from ghsl_df);""")

res.df()

Assuming raster grid resolution :  100
Using base H3 resolution :  12
Count 
ghsl_h3 = conn.table("ghsl_h3")

print(ghsl_h3.count())
ghsl_h3.head()



┌───────┐
│ 67835 │
└───────┘





┏━━━━━━━━┳━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━━┓
┃ band   ┃ h3_id           ┃ x         ┃ y         ┃
┡━━━━━━━━╇━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━━┩
│ uint16 │ string          │ float64   │ float64   │
├────────┼─────────────────┼───────────┼───────────┤
│    113 │ 8c555aad8160bff │ -8.217905 │ 16.247257 │
│    133 │ 8c555aad8b9e7ff │ -8.216861 │ 16.246436 │
│    112 │ 8c555aacbd413ff │ -8.169435 │ 16.232471 │
│     51 │ 8c5429101a725ff │ -9.551406 │ 16.221792 │
│     80 │ 8c5429101ae55ff │ -9.552404 │ 16.220971 │
└────────┴─────────────────┴───────────┴───────────┘











ghsl_h3.h3_id.nunique()











┌───────┐
│ 67835 │
└───────┘













conn.list_tables()









['dwg', 'ghsl_h3']













# ghsl_h3 = conn.create_view("ghsl_h3",obj=ghsl_h3.group_by("h3_id").agg(band=_.band.sum(),x=_.x.first(),y=_.y.first()))
ghsl_h3 = ghsl_h3.group_by("h3_id").agg(band=_.band.sum(),x=_.x.first(),y=_.y.first()).alias("ghsl_h3")













conn.sql("select count(*) as N from ghsl_h3;")









┏━━━━━━━┓
┃ N     ┃
┡━━━━━━━┩
│ int64 │
├───────┤
│ 67835 │
└───────┘











ghsl_df













  
    

      
      band
      x
      y
    

  
  
    

      142378
      113
      -8.217905
      16.247257
    

    

      152379
      133
      -8.216861
      16.246436
    

    

      322425
      112
      -8.169435
      16.232471
    

    

      451074
      51
      -9.551406
      16.221792
    

    

      461073
      80
      -9.552404
      16.220971
    

    

      ...
      ...
      ...
      ...
    

    

      98557016
      1560
      6.591825
      0.116867
    

    

      98567015
      468
      6.590827
      0.116058
    

    

      98567016
      1471
      6.591825
      0.116058
    

    

      99116953
      12
      6.528960
      0.071576
    

    

      99416960
      40
      6.535942
      0.047313
    

  



67835 rows × 3 columns










Coordinates to plot only points






h3_layer_ghsl_h3 = pdk.Layer(
    "H3HexagonLayer",
    ghsl_h3.schema(),#.execute(limit=100_000),
    pickable=True,
    stroked=True,
    filled=True,
    opacity=1,
    extruded=False,
    get_hexagon= "h3_id",
    get_fill_color = [200,200,100,255],
    get_line_color=[0, 0, 0, 100],
    line_width_min_pixels=0,
    line_width_max_pixels=1,
)


points_layer = pdk.Layer(
    "ScatterplotLayer",
    ghsl_df,
    pickable=False,
    opacity=0.8,
    stroked=False,
    filled=True,
    radius_scale=1,
    radius_min_pixels=1,
    radius_max_pixels=20,
    line_width_min_pixels=1,
    get_position = ["x","y"],
    get_radius=50,
    get_fill_color=[255, 140, 0, 255],
    get_line_color=[0, 0, 0, 0],
)













view_state_global = pdk.data_utils.compute_view(ghsl_df[["x","y"]])

# Create mapdeck object
deck_map_ghsl_h3 = pdk.Deck(layers=[
    # h3_layer_ghsl_h3,
    points_layer
    ]
             ,initial_view_state=view_state_global
             ,tooltip= {"text": "Value :  {band}"}
             ,map_style="road"
             )

deck_map_ghsl_h3.to_html("../../deck_maps/deck_ghsl_h3_grid.html",iframe_height=800,open_browser=False)










        
    








Selecting an area






dwg.filter(_.country.re_search("Nigeria")).select("gid_0","country")#.sql("Select ST_GeomFromWKB(geometry) as geom from dwg;")









┏━━━━━━━━┳━━━━━━━━━┓
┃ gid_0  ┃ country ┃
┡━━━━━━━━╇━━━━━━━━━┩
│ string │ string  │
├────────┼─────────┤
│ NGA    │ Nigeria │
│ NGA    │ Nigeria │
│ NGA    │ Nigeria │
│ NGA    │ Nigeria │
│ NGA    │ Nigeria │
│ NGA    │ Nigeria │
│ NGA    │ Nigeria │
│ NGA    │ Nigeria │
│ NGA    │ Nigeria │
│ NGA    │ Nigeria │
│ …      │ …       │
└────────┴─────────┘











selected_country = "NGA"















Rescaling approach to econ data






conn.drop_view("country",force=True)













country = conn.sql(f"""
SELECT sel_country.* EXCLUDE (color,centr,geom,radius,x,y),
         ghsl.*
         FROM ghsl_h3 AS ghsl
         JOIN (SELECT * EXCLUDE geometry, ST_GeomFromWKB(geometry) as geom FROM dwg where gid_0='{selected_country}') AS sel_country
         ON ST_Intersects(ST_Point(ghsl.x,ghsl.y),sel_country.geom);
         
""")













print(country.count())
country.head(3)











┌───────┐
│ 18401 │
└───────┘





┏━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓
┃ gid_0  ┃ country ┃ gid_1    ┃ grp_usd_2015 ┃ services_usd_2015 ┃ manufacturing_usd_2015 ┃ agriculture_usd_2015 ┃ band   ┃ h3_id           ┃ x        ┃ y         ┃
┡━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩
│ string │ string  │ string   │ float64      │ float64           │ float64                │ float64              │ uint16 │ string          │ float64  │ float64   │
├────────┼─────────┼──────────┼──────────────┼───────────────────┼────────────────────────┼──────────────────────┼────────┼─────────────────┼──────────┼───────────┤
│ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │    209 │ 8c580a0a26087ff │ 8.643287 │ 12.237987 │
│ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │     83 │ 8c580a0f19907ff │ 8.780907 │ 12.237171 │
│ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │     50 │ 8c580a0f19989ff │ 8.781919 │ 12.237171 │
└────────┴─────────┴──────────┴──────────────┴───────────────────┴────────────────────────┴──────────────────────┴────────┴─────────────────┴──────────┴───────────┘











gid_1_band = country.group_by("gid_1").agg(total_band=_.band.sum())
gid_1_band









┏━━━━━━━━━━┳━━━━━━━━━━━━┓
┃ gid_1    ┃ total_band ┃
┡━━━━━━━━━━╇━━━━━━━━━━━━┩
│ string   │ int64      │
├──────────┼────────────┤
│ NGA.23_1 │    1546077 │
│ NGA.28_1 │    4120591 │
│ NGA.12_1 │     407735 │
│ NGA.11_1 │     247390 │
│ NGA.5_1  │     401362 │
│ NGA.19_1 │    3657441 │
│ NGA.13_1 │     176415 │
│ NGA.10_1 │    2777156 │
│ NGA.33_1 │    4243386 │
│ NGA.16_1 │     199230 │
│ …        │          … │
└──────────┴────────────┘











country_gdf = (country
           .join(gid_1_band,country.gid_1==gid_1_band.gid_1,how="left")
           .mutate(band_frac=_.band/_.total_band,build_gdp=_.services_usd_2015+_.manufacturing_usd_2015)
           .mutate(h3_id_gdp=_.band_frac*_.build_gdp)).to_pandas()

















print(country_gdf.shape)
country_gdf.head()









(18401, 16)









  
    

      
      gid_0
      country
      gid_1
      grp_usd_2015
      services_usd_2015
      manufacturing_usd_2015
      agriculture_usd_2015
      band
      h3_id
      x
      y
      gid_1_right
      total_band
      band_frac
      build_gdp
      h3_id_gdp
    

  
  
    

      0
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8c580a0a26087ff
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
    

    

      1
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      83
      8c580a0f19907ff
      8.780907
      12.237171
      NGA.20_1
      3625664
      0.000023
      9.189241e+09
      210363.400722
    

    

      2
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      50
      8c580a0f19989ff
      8.781919
      12.237171
      NGA.20_1
      3625664
      0.000014
      9.189241e+09
      126724.940194
    

    

      3
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      371
      8c581d228cc19ff
      8.074492
      12.236355
      NGA.20_1
      3625664
      0.000102
      9.189241e+09
      940299.056241
    

    

      4
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      55
      8c581d2680d2dff
      8.237413
      12.235539
      NGA.20_1
      3625664
      0.000015
      9.189241e+09
      139397.434214
    

  









Check that things add up






country_gdf[["gid_1","band_frac"]].groupby("gid_1").agg("sum")













  
    

      
      band_frac
    

    

      gid_1
      
    

  
  
    

      NGA.10_1
      1.0
    

    

      NGA.11_1
      1.0
    

    

      NGA.12_1
      1.0
    

    

      NGA.13_1
      1.0
    

    

      NGA.15_1
      1.0
    

    

      NGA.16_1
      1.0
    

    

      NGA.18_1
      1.0
    

    

      NGA.19_1
      1.0
    

    

      NGA.20_1
      1.0
    

    

      NGA.23_1
      1.0
    

    

      NGA.27_1
      1.0
    

    

      NGA.28_1
      1.0
    

    

      NGA.29_1
      1.0
    

    

      NGA.30_1
      1.0
    

    

      NGA.31_1
      1.0
    

    

      NGA.33_1
      1.0
    

    

      NGA.37_1
      1.0
    

    

      NGA.3_1
      1.0
    

    

      NGA.4_1
      1.0
    

    

      NGA.5_1
      1.0
    

    

      NGA.6_1
      1.0
    

    

      NGA.9_1
      1.0
    

  











country_gdf[["h3_id_gdp"]].sum()









h3_id_gdp    2.000056e+11
dtype: float64













country_gdf[["h3_id_gdp"]].sum()/country_gdf.shape[0]









h3_id_gdp    1.086928e+07
dtype: float64













# country_gdf[["h3_id_gdp"]].sum()
# country_gdf.groupby("gid_1").first()













Mapping






import scalenav.plotting as snplt













col_pal = "viridis"

col_pal = pypal.load_cmap(col_pal)













cols =  snplt.cmap((country_gdf.h3_id_gdp),palette=col_pal) #pd.Series([[255*j for j in x] for x in bmlunge(compact_geo_downscaled["log_value"])])
country_gdf["color"] = cols
country_gdf.head(3)
# print(cols)
# print(unique_vals)









Max input : 329053197.20, palette colors : 10









  
    

      
      gid_0
      country
      gid_1
      grp_usd_2015
      services_usd_2015
      manufacturing_usd_2015
      agriculture_usd_2015
      band
      h3_id
      x
      y
      gid_1_right
      total_band
      band_frac
      build_gdp
      h3_id_gdp
      color
    

  
  
    

      0
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8c580a0a26087ff
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
    

    

      1
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      83
      8c580a0f19907ff
      8.780907
      12.237171
      NGA.20_1
      3625664
      0.000023
      9.189241e+09
      210363.400722
      [72, 33, 115, 255]
    

    

      2
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      50
      8c580a0f19989ff
      8.781919
      12.237171
      NGA.20_1
      3625664
      0.000014
      9.189241e+09
      126724.940194
      [72, 33, 115, 255]
    

  











view_state = pdk.data_utils.compute_view(points=country_gdf[["x","y"]])













# # this map is not really usefull at high resolutions

# h3_layer_down = pdk.Layer(
#     "H3HexagonLayer",
#     country_gdf[["gid_0","gid_1","h3_id","total_band","h3_id_gdp","color"]],
#     pickable=True,
#     stroked=True,
#     filled=True,
#     opacity=.8,
#     extruded=False,
#     get_hexagon= "h3_id",
#     get_fill_color = "color",
#     get_line_color= [0, 0, 0, 100],
#     line_width_min_pixels=0,
#     line_width_max_pixels=1,
# )


# layer_down = pdk.Layer(
#     "ScatterplotLayer",
#     country_gdf[["h3_id_gdp","x","y"]],
#     pickable=False,
#     opacity=0.4,
#     stroked=False,
#     filled=True,
#     radius_scale=6,
#     radius_min_pixels=1,
#     radius_max_pixels=10,
#     line_width_min_pixels=1,
#     get_position = ["x","y"],
#     get_radius=50,
#     get_fill_color=[255, 140, 0, 255],
#     get_line_color=[0, 0, 0, 0],
# )













# # Create mapdeck object
# deck_map_down = pdk.Deck(layers=[
#     h3_layer_down,
#     layer_down
#     ]
#              ,initial_view_state=view_state
#             #  ,tooltip= {"text": "Value :  {h3_id_gdp}"}
#              ,height=800
#              ,map_style="road"
#              )

# deck_map_down.to_html("deck_maps/deck_ghsl_h3_down_grid.html")













Gridding






rast_to_h3 = dt.rast_to_h3_map(x=country_gdf["x"].mean(),y = country_gdf["y"].mean())

# h3_res = dt.rast_to_h3["100"]["h3_res"]

# square = dt.square_poly(lat=country_gdf.loc[0,"y"],lon=country_gdf.loc[0,"x"],distance=100)

# ref_cell = h3.latlng_to_cell(lat=country_gdf.loc[0,"y"],lng=country_gdf.loc[0,"x"],res=h3_res )
# ref_cell_ij = h3.cell_to_local_ij(origin=ref_cell,h=ref_cell)

# cells = h3.geo_to_cells(square,res=h3_res )

# cells_ij = [h3.cell_to_local_ij(origin=ref_cell,h=cell) for cell in cells]

# neighbs = rast_to_h3["100"]["nn"]
# # [(cell_i-ref_cell_ij[0],cell_j-ref_cell_ij[1]) for (cell_i,cell_j) in cells_ij]

# ####
# neighb_cells = [h3.local_ij_to_cell(origin=ref_cell,i=neighb[0]+ref_cell_ij[0],j=neighb[1]+ref_cell_ij[1]) for neighb in neighbs]

# reconstructed = dt.centre_cell_to_square(ref_cell,grid_param=grid_param,neighbs=neighbs)









Using angles for meter grid.
Using angles for meter grid.
Using angles for meter grid.
Using angles for meter grid.
Using angles for meter grid.













# gpd.GeoSeries([
#     square,
#                h3.cells_to_h3shape(neighb_cells)
#                ,h3.cells_to_h3shape(reconstructed)
#                ],crs="EPSG:4326").explore()













# ideally, this is precomputed, but for some reason it has not worked...
neighbs = rast_to_h3[str(grid_param)]["nn"]

country_gdf["h3_gridded"] = pd.Series(country_gdf.apply(lambda row: dt.centre_cell_to_square(h3_cell = row["h3_id"],neighbs=neighbs),axis=1).tolist())













country_gdf.dropna(subset=["h3_gridded"],inplace=True)













Gridded data set






ghsl_gridded = country_gdf.explode("h3_gridded")













print(ghsl_gridded.shape)
ghsl_gridded.head(3)









(662436, 18)









  
    

      
      gid_0
      country
      gid_1
      grp_usd_2015
      services_usd_2015
      manufacturing_usd_2015
      agriculture_usd_2015
      band
      h3_id
      x
      y
      gid_1_right
      total_band
      band_frac
      build_gdp
      h3_id_gdp
      color
      h3_gridded
    

  
  
    

      0
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8c580a0a26087ff
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a26081ff
    

    

      0
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8c580a0a26087ff
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a26089ff
    

    

      0
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8c580a0a26087ff
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a260adff
    

  











ghsl_gridded = ghsl_gridded[~ghsl_gridded.duplicated(subset=["h3_gridded"])]













print(ghsl_gridded.shape)
ghsl_gridded.reset_index(inplace=True,drop=True)









(632998, 18)













ghsl_gridded.head(3)













  
    

      
      gid_0
      country
      gid_1
      grp_usd_2015
      services_usd_2015
      manufacturing_usd_2015
      agriculture_usd_2015
      band
      h3_id
      x
      y
      gid_1_right
      total_band
      band_frac
      build_gdp
      h3_id_gdp
      color
      h3_gridded
    

  
  
    

      0
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8c580a0a26087ff
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a26081ff
    

    

      1
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8c580a0a26087ff
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a26089ff
    

    

      2
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8c580a0a26087ff
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a260adff
    

  











ghsl_gridded.set_index("h3_id",inplace=True,drop=True)













ghsl_gridded = ghsl_gridded.join(ghsl_gridded.value_counts("h3_id"),how="left",sort=False,validate="m:1")













ghsl_gridded["h3_id_gdp_gridded"] = np.round(ghsl_gridded["h3_id_gdp"]/ghsl_gridded["count"])













ghsl_gridded.reset_index(inplace=True,drop=False)













ghsl_gridded.dropna(axis=0,subset=["h3_gridded"],inplace=True)













print(ghsl_gridded.shape)
ghsl_gridded.head(3)









(632998, 20)









  
    

      
      h3_id
      gid_0
      country
      gid_1
      grp_usd_2015
      services_usd_2015
      manufacturing_usd_2015
      agriculture_usd_2015
      band
      x
      y
      gid_1_right
      total_band
      band_frac
      build_gdp
      h3_id_gdp
      color
      h3_gridded
      count
      h3_id_gdp_gridded
    

  
  
    

      0
      8c580a0a26087ff
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a26081ff
      36
      14714.0
    

    

      1
      8c580a0a26087ff
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a26089ff
      36
      14714.0
    

    

      2
      8c580a0a26087ff
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a260adff
      36
      14714.0
    

  











# view_state = pdk.data_utils.compute_view(points=ghsl_gridded.sample(10000)["h3_gridded"].apply(lambda cell: h3.cell_to_latlng(cell)[::-1]))

# lagos: latitude = 6.450151,longitude = 3.531193 !!!! DOES NOT WORK BECAUSE IT IS EXCLUDED FROM DOSE!!!
# Abidjan :  latitude = 5.367, longitude = -4.009
# Kano : latitude = 11.985, longitude = 8.533

zoom_lat=11.985
zoom_lon=8.533
zoom_width=70_000

view_state_zoomed = pdk.ViewState(latitude = zoom_lat,longitude = zoom_lon,zoom=12,pitch=30,bearing=0) 













# extracting data for the plotting region : 

plotting_region = dt.square_poly(lat=zoom_lat,lon=zoom_lon,distance=zoom_width)
# gpd.GeoSeries(plotting_region,crs="EPSG:4326").explore()









Using angles for arc grid.













limits = [np.round(x,3) for x in plotting_region.bounds]
limits









[8.212, 11.664, 8.854, 12.306]













# ghsl_gridded

conn.sql("""select *, 
    array_extract(h3_cell_to_latlng(h3_gridded),2) as x_h3, 
    array_extract(h3_cell_to_latlng(h3_gridded),1) as y_h3
    from ghsl_gridded""")









┏━━━━━━━━━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━┳━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━┳━━━━━━━┳━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━┳━━━━━━━━━━━┓
┃ h3_id           ┃ gid_0  ┃ country ┃ gid_1    ┃ grp_usd_2015 ┃ services_usd_2015 ┃ manufacturing_usd_2015 ┃ agriculture_usd_2015 ┃ band   ┃ x        ┃ y         ┃ gid_1_right ┃ total_band ┃ band_frac ┃ build_gdp    ┃ h3_id_gdp     ┃ color            ┃ h3_gridded      ┃ count ┃ h3_id_gdp_gridded ┃ x_h3     ┃ y_h3      ┃
┡━━━━━━━━━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━╇━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━╇━━━━━━━╇━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━╇━━━━━━━━━━━┩
│ string          │ string │ string  │ string   │ float64      │ float64           │ float64                │ float64              │ uint16 │ float64  │ float64   │ string      │ int64      │ float64   │ float64      │ float64       │ array<int32>     │ string          │ int64 │ float64           │ float64  │ float64   │
├─────────────────┼────────┼─────────┼──────────┼──────────────┼───────────────────┼────────────────────────┼──────────────────────┼────────┼──────────┼───────────┼─────────────┼────────────┼───────────┼──────────────┼───────────────┼──────────────────┼─────────────────┼───────┼───────────────────┼──────────┼───────────┤
│ 8c580a0a26087ff │ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │    209 │ 8.643287 │ 12.237987 │ NGA.20_1    │    3625664 │  0.000058 │ 9.189241e+09 │ 529710.250012 │ [72, 33, ... +2] │ 8c580a0a26081ff │    36 │           14714.0 │ 8.643200 │ 12.238081 │
│ 8c580a0a26087ff │ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │    209 │ 8.643287 │ 12.237987 │ NGA.20_1    │    3625664 │  0.000058 │ 9.189241e+09 │ 529710.250012 │ [72, 33, ... +2] │ 8c580a0a26089ff │    36 │           14714.0 │ 8.643154 │ 12.238244 │
│ 8c580a0a26087ff │ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │    209 │ 8.643287 │ 12.237987 │ NGA.20_1    │    3625664 │  0.000058 │ 9.189241e+09 │ 529710.250012 │ [72, 33, ... +2] │ 8c580a0a260adff │    36 │           14714.0 │ 8.642837 │ 12.237747 │
│ 8c580a0a26087ff │ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │    209 │ 8.643287 │ 12.237987 │ NGA.20_1    │    3625664 │  0.000058 │ 9.189241e+09 │ 529710.250012 │ [72, 33, ... +2] │ 8c580a0a26199ff │    36 │           14714.0 │ 8.642762 │ 12.237472 │
│ 8c580a0a26087ff │ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │    209 │ 8.643287 │ 12.237987 │ NGA.20_1    │    3625664 │  0.000058 │ 9.189241e+09 │ 529710.250012 │ [72, 33, ... +2] │ 8c580a0a26091ff │    36 │           14714.0 │ 8.643658 │ 12.238089 │
│ 8c580a0a26087ff │ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │    209 │ 8.643287 │ 12.237987 │ NGA.20_1    │    3625664 │  0.000058 │ 9.189241e+09 │ 529710.250012 │ [72, 33, ... +2] │ 8c580a0a260e3ff │    36 │           14714.0 │ 8.642911 │ 12.238021 │
│ 8c580a0a26087ff │ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │    209 │ 8.643287 │ 12.237987 │ NGA.20_1    │    3625664 │  0.000058 │ 9.189241e+09 │ 529710.250012 │ [72, 33, ... +2] │ 8c580a0a260e7ff │    36 │           14714.0 │ 8.642790 │ 12.237910 │
│ 8c580a0a26087ff │ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │    209 │ 8.643287 │ 12.237987 │ NGA.20_1    │    3625664 │  0.000058 │ 9.189241e+09 │ 529710.250012 │ [72, 33, ... +2] │ 8c580a0a26099ff │    36 │           14714.0 │ 8.643611 │ 12.238252 │
│ 8c580a0a26087ff │ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │    209 │ 8.643287 │ 12.237987 │ NGA.20_1    │    3625664 │  0.000058 │ 9.189241e+09 │ 529710.250012 │ [72, 33, ... +2] │ 8c580a0a260bdff │    36 │           14714.0 │ 8.643294 │ 12.237754 │
│ 8c580a0a26087ff │ NGA    │ Nigeria │ NGA.20_1 │ 1.151516e+10 │      8.028410e+09 │           1.160831e+09 │         2.325922e+09 │    209 │ 8.643287 │ 12.237987 │ NGA.20_1    │    3625664 │  0.000058 │ 9.189241e+09 │ 529710.250012 │ [72, 33, ... +2] │ 8c580a0a260a9ff │    36 │           14714.0 │ 8.642958 │ 12.237858 │
│ …               │ …      │ …       │ …        │            … │                 … │                      … │                    … │      … │        … │         … │ …           │          … │         … │            … │             … │ …                │ …               │     … │                 … │        … │         … │
└─────────────────┴────────┴─────────┴──────────┴──────────────┴───────────────────┴────────────────────────┴──────────────────────┴────────┴──────────┴───────────┴─────────────┴────────────┴───────────┴──────────────┴───────────────┴──────────────────┴─────────────────┴───────┴───────────────────┴──────────┴───────────┘











ghsl_gridded_zoom = conn.sql(
f"""
Select * from 
    (select *
    EXCLUDE (band,x,y,gid_1_right),
    band::DECIMAL(6,2) as band,
    array_extract(h3_cell_to_latlng(h3_gridded),2) as x_h3, 
    array_extract(h3_cell_to_latlng(h3_gridded),1) as y_h3
    from ghsl_gridded)
where (x_h3 <= {limits[2]}) and (x_h3 >= {limits[0]}) and (y_h3 <= {limits[3]}) and (y_h3 >= {limits[1]});
""").execute()


print(ghsl_gridded_zoom.shape)
ghsl_gridded_zoom.head(3)









(84208, 19)









  
    

      
      h3_id
      gid_0
      country
      gid_1
      grp_usd_2015
      services_usd_2015
      manufacturing_usd_2015
      agriculture_usd_2015
      total_band
      band_frac
      build_gdp
      h3_id_gdp
      color
      h3_gridded
      count
      h3_id_gdp_gridded
      band
      x_h3
      y_h3
    

  
  
    

      0
      8c580a0a26087ff
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a26081ff
      36
      14714.0
      209.00
      8.643200
      12.238081
    

    

      1
      8c580a0a26087ff
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a26089ff
      36
      14714.0
      209.00
      8.643154
      12.238244
    

    

      2
      8c580a0a26087ff
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a260adff
      36
      14714.0
      209.00
      8.642837
      12.237747
    

  











band_pal = pypal.load_cmap("Bmsurface")
ghsl_gridded_zoom["band_color"] = cmap(ghsl_gridded_zoom.band_frac,palette=band_pal)
cmap([0.002457481],palette=band_pal)









Max input : 0.00, palette colors : 6
Max input : 0.00, palette colors : 6





[[56, 65, 87, 255]]













ghsl_gridded_zoom.columns









Index(['h3_id', 'gid_0', 'country', 'gid_1', 'grp_usd_2015',
       'services_usd_2015', 'manufacturing_usd_2015', 'agriculture_usd_2015',
       'total_band', 'band_frac', 'build_gdp', 'h3_id_gdp', 'color',
       'h3_gridded', 'count', 'h3_id_gdp_gridded', 'band', 'x_h3', 'y_h3',
       'band_color'],
      dtype='object')













h3_layer_ghsl_h3_gridded = pdk.Layer(
    "H3HexagonLayer",
    # ghsl_gridded.loc[0:200_000],
    ghsl_gridded_zoom[["band_frac","h3_gridded","color","band_color"]],
    pickable=True,
    stroked=True,
    filled=True,
    opacity=.5,
    extruded=False,
    get_hexagon= "h3_gridded",
    get_fill_color = "band_color",
    get_line_color= [0, 0, 0, 100],
    line_width_min_pixels=0,
    line_width_max_pixels=1,
)


# layer_down = pdk.Layer(
#     "ScatterplotLayer",
#     country_gdf[["h3_id_gdp","x","y"]],
#     pickable=False,
#     opacity=0.8,
#     stroked=False,
#     filled=True,
#     radius_scale=10,
#     radius_min_pixels=1,
#     radius_max_pixels=20,
#     line_width_min_pixels=0,
#     line_width_max_pixels=0,
#     get_position = ["x","y"],
#     get_radius=25,
#     get_fill_color=[255, 140, 0, 255],
#     get_line_color=[0, 0, 0, 0],
# )













# Create mapdeck object

deck_map_ghsl_h3_gridded = pdk.Deck(layers=[h3_layer_ghsl_h3_gridded, 
                                            # layer_down
                                            ]
             ,initial_view_state=view_state_zoomed
             ,tooltip= {"text": "Value :  {band_frac}"}
             ,map_style="road"
             )

deck_map_ghsl_h3_gridded.to_html("../../deck_maps/deck_ghsl_nres_h3_gridded.html",open_browser=True,iframe_height=800)










        
    










Changing Scales






ghsl_gridded.rename(columns={"h3_id_gdp_gridded" : "h3_id_gdp_gridded_var",
                             "h3_id" : "h3_gridded_parent"
                             },inplace=True)













ghsl_gridded.rename(columns={"h3_gridded" : "h3_id"
                             },inplace=True)













ghsl_gridded.head(3)













  
    

      
      h3_gridded_parent
      gid_0
      country
      gid_1
      grp_usd_2015
      services_usd_2015
      manufacturing_usd_2015
      agriculture_usd_2015
      band
      x
      y
      gid_1_right
      total_band
      band_frac
      build_gdp
      h3_id_gdp
      color
      h3_id
      count
      h3_id_gdp_gridded_var
    

  
  
    

      0
      8c580a0a26087ff
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a26081ff
      36
      14714.0
    

    

      1
      8c580a0a26087ff
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a26089ff
      36
      14714.0
    

    

      2
      8c580a0a26087ff
      NGA
      Nigeria
      NGA.20_1
      1.151516e+10
      8.028410e+09
      1.160831e+09
      2.325922e+09
      209
      8.643287
      12.237987
      NGA.20_1
      3625664
      0.000058
      9.189241e+09
      529710.250012
      [72, 33, 115, 255]
      8c580a0a260adff
      36
      14714.0
    

  











h3_gridded_new_scale = sd.change_scale(ghsl_gridded,level=-4)

# h3_gridded_new_scale = sd.change_scale(h3_gridded_new_scale,level=1)













print(h3_gridded_new_scale.shape)
h3_gridded_new_scale.head()









(4691, 3)









  
    

      
      h3_id
      child_cells
      h3_id_gdp_gridded_var
    

  
  
    

      0
      8858022411fffff
      [89580224107ffff]
      9140961.0
    

    

      1
      8858022413fffff
      [89580224133ffff, 8958022413bffff]
      76653316.0
    

    

      2
      8858022c8bfffff
      [8958022c8a3ffff, 8958022c8b7ffff]
      5570626.0
    

    

      3
      8858094823fffff
      [89580948223ffff, 89580948227ffff, 8958094822b...
      201764238.0
    

    

      4
      88580a0053fffff
      [89580a00523ffff, 89580a0053bffff]
      2825964.0
    

  











h3_gridded_new_scale["color"] = cmap(np.log1p(h3_gridded_new_scale.h3_id_gdp_gridded_var),palette=col_pal)









Max input : 21.78, palette colors : 10













#####

h3_layer_ghsl_h3_gridded_new_scale = pdk.Layer(
    "H3HexagonLayer",
    h3_gridded_new_scale,
    pickable=True,
    stroked=True,
    filled=True,
    opacity=.5,
    extruded=False,
    get_hexagon= "h3_id",
    get_fill_color = "color",
    get_line_color= [0, 0, 0, 100],
    line_width_min_pixels=0,
    line_width_max_pixels=1,
)













# Create mapdeck object
deck_map_ghsl_h3_gridded_rescaled = pdk.Deck(layers=[h3_layer_ghsl_h3_gridded_new_scale]
             ,initial_view_state=view_state
             ,tooltip= {"text": "Value :  {h3_id_gdp_gridded_var}"}
             ,map_style="road"
             )

deck_map_ghsl_h3_gridded_rescaled.to_html("../../deck_maps/deck_ghsl_nres_h3_gridded_new_scale.html",iframe_height=800,open_browser=True)