Skip to content

Commit

Permalink
Bugfix/58/h3 parent approach (#69)
Browse files Browse the repository at this point in the history 
* Add failing test

* Implement recursive polyfill with higher hex resolution

* Implement support for touches (intersects)

* Remove finding exterior hexagons for touches

* Re-add MultiPolygon support to generate_h3_ids

* Convert typing from list to set when applicable

* Re-add deleted ingest cli

* Rename recursive parent to get parent

* Remove _get_parent and add Literal to spatial_join_method typing

* Add Kenya administrative boundary example

Andres notebook raised the initial bug for the small geometries. The notebook reproduced the issue and now demonstrates the fix.

---------

Co-authored-by: Anthony Lukach <anthonylukach@gmail.com>
  • Loading branch information
@zacdezgeo @alukach
zacdezgeo and alukach authored Sep 30, 2024
1 parent 966084a commit e2dfb0f
Show file tree
Hide file tree
Showing 4 changed files with 477 additions and 45 deletions.
319 changes: 319 additions & 0 deletions notebooks/space2stats_api_adm_example.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,319 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [],
"source": [
"import requests\n",
"import geopandas as gpd\n",
"import pandas as pd\n",
"import json\n",
"from shapely.geometry import shape\n",
"import numpy as np\n",
"from lonboard import Map, ScatterplotLayer"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [],
"source": [
"BASE_URL = \"http://127.0.0.1:8000\"\n",
"FIELDS_ENDPOINT = f\"{BASE_URL}/fields\"\n",
"SUMMARY_ENDPOINT = f\"{BASE_URL}/summary\""
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [],
"source": [
"def fetch_admin_boundaries(iso3: str, adm: str) -> gpd.GeoDataFrame:\n",
" \"\"\"Fetch administrative boundaries from GeoBoundaries API.\"\"\"\n",
" url = f'https://www.geoboundaries.org/api/current/gbOpen/{iso3}/{adm}/'\n",
" res = requests.get(url).json()\n",
" return gpd.read_file(res['gjDownloadURL'])\n",
"\n",
"def fetch_summary_data(feature: Dict) -> pd.DataFrame:\n",
" \"\"\"Fetch summary data for each administrative feature.\"\"\"\n",
" request_payload = {\n",
" \"aoi\": feature,\n",
" \"spatial_join_method\": \"touches\",\n",
" \"fields\": [\"sum_pop_2020\"], \n",
" \"geometry\": \"point\"\n",
" }\n",
" response = requests.post(SUMMARY_ENDPOINT, json=request_payload)\n",
" if response.status_code != 200:\n",
" raise Exception(f\"Failed to get summary: {response.text}\")\n",
" \n",
" summary_data = response.json()\n",
" if not summary_data:\n",
" print(f\"Failed to get summary for {feature['id']}\")\n",
" return pd.DataFrame() # Return an empty DataFrame if no data\n",
" \n",
" df = pd.DataFrame(summary_data)\n",
" df['adm_id'] = int(feature['id'])\n",
" df['adm_name'] = feature['properties']['shapeName']\n",
" df['geometry'] = df['geometry'].apply(lambda geom: shape(geom))\n",
" return df"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Turkana 0\n",
"Marsabit 1\n",
"Mandera 2\n",
"Wajir 3\n",
"West Pokot 4\n",
"Samburu 5\n",
"Isiolo 6\n",
"Baringo 7\n",
"Elgeyo-Marakwet 8\n",
"Trans Nzoia 9\n",
"Bungoma 10\n",
"Garissa 11\n",
"Uasin Gishu 12\n",
"Kakamega 13\n",
"Laikipia 14\n",
"Busia 15\n",
"Meru 16\n",
"Nandi 17\n",
"Siaya 18\n",
"Nakuru 19\n",
"Vihiga 20\n",
"Nyandarua 21\n",
"Tharaka 22\n",
"Kericho 23\n",
"Kisumu 24\n",
"Nyeri 25\n",
"Tana River 26\n",
"Kitui 27\n",
"Kirinyaga 28\n",
"Embu 29\n",
"Homa Bay 30\n",
"Bomet 31\n",
"Nyamira 32\n",
"Narok 33\n",
"Kisii 34\n",
"Murang'a 35\n",
"Migori 36\n",
"Kiambu 37\n",
"Machakos 38\n",
"Kajiado 39\n",
"Nairobi 40\n",
"Makueni 41\n",
"Lamu 42\n",
"Kilifi 43\n",
"Taita Taveta 44\n",
"Kwale 45\n",
"Mombasa 46\n"
]
},
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>hex_id</th>\n",
" <th>geometry</th>\n",
" <th>sum_pop_2020</th>\n",
" <th>adm_id</th>\n",
" <th>adm_name</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>866a4a407ffffff</td>\n",
" <td>POINT (35.68873 4.94892)</td>\n",
" <td>554.655021</td>\n",
" <td>0</td>\n",
" <td>Turkana</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>866a4a417ffffff</td>\n",
" <td>POINT (35.59512 1.81273)</td>\n",
" <td>207.977046</td>\n",
" <td>0</td>\n",
" <td>Turkana</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>866a4a41fffffff</td>\n",
" <td>POINT (35.17288 5.20383)</td>\n",
" <td>637.218817</td>\n",
" <td>0</td>\n",
" <td>Turkana</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>866a4a427ffffff</td>\n",
" <td>POINT (35.61627 2.62098)</td>\n",
" <td>948.876730</td>\n",
" <td>0</td>\n",
" <td>Turkana</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>866a4a437ffffff</td>\n",
" <td>POINT (35.89027 3.48441)</td>\n",
" <td>112.805780</td>\n",
" <td>0</td>\n",
" <td>Turkana</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" hex_id geometry sum_pop_2020 adm_id adm_name\n",
"0 866a4a407ffffff POINT (35.68873 4.94892) 554.655021 0 Turkana\n",
"1 866a4a417ffffff POINT (35.59512 1.81273) 207.977046 0 Turkana\n",
"2 866a4a41fffffff POINT (35.17288 5.20383) 637.218817 0 Turkana\n",
"3 866a4a427ffffff POINT (35.61627 2.62098) 948.876730 0 Turkana\n",
"4 866a4a437ffffff POINT (35.89027 3.48441) 112.805780 0 Turkana"
]
},
"execution_count": 55,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ISO3 = 'KEN'\n",
"ADM = 'ADM1'\n",
"adm_boundaries = fetch_admin_boundaries(ISO3, ADM)\n",
"geojson_str = adm_boundaries.to_json()\n",
"adm_geojson = json.loads(geojson_str)\n",
"adm_features = adm_geojson['features']\n",
"\n",
"gdfs = []\n",
"for i, feature in enumerate(adm_features):\n",
" df = fetch_summary_data(feature)\n",
" if not df.empty:\n",
" gdfs.append(gpd.GeoDataFrame(df, geometry='geometry', crs='EPSG:4326'))\n",
" print(feature['properties']['shapeName'], i)\n",
"\n",
"# Concatenate all GeoDataFrames into a single GeoDataFrame\n",
"gdf = pd.concat(gdfs, ignore_index=True)\n",
"\n",
"# Display the GeoDataFrame structure\n",
"gdf.head()"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [],
"source": [
"# Define custom breaks and corresponding RGBA colors for visualization\n",
"breaks = [0, 1, 1000, 10000, 50000, 100000, 200000, gdf[\"sum_pop_2020\"].max()]\n",
"colors = np.array([\n",
" [211, 211, 211, 255], # Light gray for 0\n",
" [255, 255, 0, 255], # Yellow for 1-1000\n",
" [255, 165, 0, 255], # Orange for 1000-10000\n",
" [255, 0, 0, 255], # Red for 10000-50000\n",
" [128, 0, 128, 255], # Purple for 50000-100000\n",
" [0, 0, 255, 255], # Blue for 100000-200000\n",
" [0, 0, 139, 255], # Dark blue for 200000+\n",
"])\n",
"\n",
"def assign_color(value: float) -> list:\n",
" \"\"\"Assign colors based on population value.\"\"\"\n",
" for i in range(len(breaks) - 1):\n",
" if breaks[i] <= value < breaks[i + 1]:\n",
" return colors[i].tolist() # Convert to list\n",
" return colors[-1].tolist() # In case value exceeds all breaks"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "63b834db2cde4dea96705c9b656cea45",
"version_major": 2,
"version_minor": 1
},
"text/plain": [
"Map(layers=[ScatterplotLayer(get_fill_color=<pyarrow.lib.FixedSizeListArray object at 0x13b3b78e0>\n",
"[\n",
" [\n",
" 2…"
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Map sum_pop_2020 values to colors using the custom function\n",
"gdf['color'] = gdf[\"sum_pop_2020\"].apply(assign_color)\n",
"\n",
"# Flatten the color list into a 2D array for ScatterplotLayer\n",
"color_list = np.array(gdf['color'].tolist(), dtype=np.uint8)\n",
"\n",
"# Create the scatterplot layer with the assigned colors\n",
"layer = ScatterplotLayer.from_geopandas(gdf, get_radius=2000, get_fill_color=color_list)\n",
"\n",
"# Create the map with the scatterplot layer\n",
"m = Map(layer)\n",
"m"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "spacestats",
"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.11.9"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Loading

0 comments on commit e2dfb0f

Please sign in to comment.