Uzbekistan Boundary Comparison

import sys, os, importlib, json, multiprocessing
import folium, shapely, rasterio

import pandas as pd
import geopandas as gpd
import numpy as np

from shapely.geometry import Polygon, Point, mapping
from shapely.ops import unary_union
from urllib.request import urlopen

sys.path.insert(0, "../../gostrocks/src/")
import GOSTRocks.ntlMisc as ntl
import GOSTRocks.rasterMisc as rMisc

sys.path.append("../src")
from GOSTboundaries import boundary_helper

%load_ext autoreload
%autoreload 2

geobounds_url = 'https://www.geoboundaries.org/api/current/gbOpen/{iso3}/ADM{lvl}/'
sel_iso3 = 'UZB'
h3_lvl = 6
official_wb_bounds = f"/home/wb411133/projects/BOUNDARIES/Data/WB_med_res/{sel_iso3}_medium_res.geojson"
high_res_bounds = '/home/wb411133/projects/BOUNDARIES/Data/HighRes/WB_GAD_ADM_UZB/shp/WB_GAD_ADM2_UZB.shp'
output_folder = "/home/wb411133/projects/BOUNDARIES/"
html_maps_folder = "/home/wb411133/Code/Boundary_Comparison/docs/_static"
# Open official World Bank boundaries
selWB = gpd.read_file(official_wb_bounds)

inH = gpd.read_file(high_res_bounds)
selWB = selWB.to_crs(4326)
selH = inH.loc[inH['ISO_A3'] == sel_iso3]

# Something dumb is broken somewhere in the geometries
selWB['geometry'] = selWB['geometry'].apply(lambda x: x.buffer(0))
inH['geometry'] = inH['geometry'].apply(lambda x: x.buffer(0))

comparer = boundary_helper.country_boundary(sel_iso3, selWB, "OBJECTID", geoBounds=selH)
comparer.run_all(run_h3_summary=True, run_comparison=False, run_zonal=True)

Matching bounds2 to bounds 1: 100%|██████████| 199/199 [00:05<00:00, 35.37it/s]
Matching h3 to bounds 1: 100%|██████████| 10596/10596 [03:38<00:00, 48.39it/s]

Index(['geometry', 'shape_id', 'med_id', 'med_per'], dtype='object')

Matching h3 to bounds 2: 100%|██████████| 10596/10596 [03:52<00:00, 45.59it/s]

Index(['geometry', 'shape_id', 'med_id', 'med_per', 'geo_match_id',
       'geo_match_per'],
      dtype='object')

comparer.write_output(os.path.join(output_folder, sel_iso3), write_slivers=False)

Comparing Boundaries

There are substantive differences between the medium-resolution boundaries and the high-resolution boundaries; these areas are most prevalent in coastal islands and small, urban administrative units.

# Tashkent
comparer.map_boundary_comparison(Point(69.23295,41.27852), 11, buffer_dist=0.1)

Make this Notebook Trusted to load map: File -> Trust Notebook

# Tra Ban Islands
comparer.map_boundary_comparison(Point(66.24777,39.73277), 10, buffer_dist=0.3)

Make this Notebook Trusted to load map: File -> Trust Notebook

This leads to variation in certain zonal statistics in these areas. In this experiment we ran the following comparisons:

1. ESA Globcover - determine the majority landcover class in each district

2. Nighttime Lights - calculate sum of lights for most recent month (2023-02)

3. H3 grid comparison - joined h3 grid cells to admin divisions and identified which h3 grid cells join the different administrative divisions based on resolution.

ESA Globcover

This ~300m2 resolution landcover dataset classifies landcover into 23 categories.

For each admin division we calculated the majority class and compared between the medium and high resolution datasets; of the 681 admin 2 divisions in Vietnam, 12 have a different majority landcover class.

cur_plot = comparer.static_map_lc()
cur_plot

Nighttime Lights

We calculated sum of lights for each dataset, and then compared them, as % change from medium-resolution boundary to high-resolution.

ntl_summary = comparer.ntl_summary()
ntl_summary[0]

H3 grid comparison

Using h3 grid level 6 (~36km2), there are 10596 individual h3 grid cells in Vietnam. Of those, 12 are identified as having different parental admin boundaries based on resolutions

mapp = comparer.static_map_h3()
mapp