Skip to content

Commit 040cda6

Browse files
bpstewarbpstewar
committed
FATHOM updates
1 parent a245cd9 commit 040cda6

14 files changed

+3444
-3235
lines changed

‎notebooks/BUILT_compare_buildings_wsf.ipynb‎

Lines changed: 298 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,298 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": 30,
6+
"metadata": {},
7+
"outputs": [
8+
{
9+
"name": "stdout",
10+
"output_type": "stream",
11+
"text": [
12+
"The autoreload extension is already loaded. To reload it, use:\n",
13+
" %reload_ext autoreload\n"
14+
]
15+
}
16+
],
17+
"source": [
18+
"import sys, os\n",
19+
"import rasterio\n",
20+
"\n",
21+
"import pandas as pd\n",
22+
"import geopandas as gpd\n",
23+
"import numpy as np\n",
24+
"\n",
25+
"from shapely.wkt import loads\n",
26+
"\n",
27+
"sys.path.insert(0, \"../src\")\n",
28+
"\n",
29+
"import GOSTrocks.rasterMisc as rMisc\n",
30+
"import GOSTrocks.dataMisc as dMisc\n",
31+
"\n",
32+
"%load_ext autoreload\n",
33+
"%autoreload 2"
34+
]
35+
},
36+
{
37+
"cell_type": "code",
38+
"execution_count": 5,
39+
"metadata": {},
40+
"outputs": [
41+
{
42+
"name": "stderr",
43+
"output_type": "stream",
44+
"text": [
45+
"C:\\Users\\WB411133\\AppData\\Local\\Temp\\ipykernel_25180\\3039392467.py:14: FutureWarning: The geopandas.dataset module is deprecated and will be removed in GeoPandas 1.0. You can get the original 'naturalearth_lowres' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.\n",
46+
" world_filepath = gpd.datasets.get_path('naturalearth_lowres')\n"
47+
]
48+
}
49+
],
50+
"source": [
51+
"# Local/input files\n",
52+
"iso3 = 'KHM'\n",
53+
"out_folder = \"c:/WBG/Work/KHM_Energy/data\"\n",
54+
"wsf_file = os.path.join(out_folder, \"WSF\", \"wsf.tif\")\n",
55+
"ghsl_file = os.path.join(out_folder, \"GHSL\", \"ghsl.tif\")\n",
56+
"overture_buildings = os.path.join(out_folder, \"overture\", \"overture_download_2024_03_29.csv\")\n",
57+
"overture_raster = os.path.join(out_folder, \"overture\", \"overture_download_2024_03_29.tif\")\n",
58+
"overture_raster_points = os.path.join(out_folder, \"overture\", \"overture_download_2024_03_29_points.tif\")\n",
59+
"for file in [wsf_file, ghsl_file]:\n",
60+
" if not os.path.exists(os.path.dirname(file)):\n",
61+
" os.makedirs(os.path.dirname(file))\n",
62+
"\n",
63+
"# get country extent from geopandas\n",
64+
"world_filepath = gpd.datasets.get_path('naturalearth_lowres')\n",
65+
"world = gpd.read_file(world_filepath)\n",
66+
"country = world[world.iso_a3 == iso3]"
67+
]
68+
},
69+
{
70+
"cell_type": "code",
71+
"execution_count": 6,
72+
"metadata": {},
73+
"outputs": [
74+
{
75+
"name": "stdout",
76+
"output_type": "stream",
77+
"text": [
78+
"FUBAR\n"
79+
]
80+
}
81+
],
82+
"source": [
83+
"\"\"\" Not working with World Bank Firewall\n",
84+
"# Download the WSF data\n",
85+
"if not os.path.exists(wsf_file):\n",
86+
" print(\"Downloading WSF data\")\n",
87+
" wsf_data, wsf_profile = dMisc.download_WSF(country, out_file = wsf_file)\n",
88+
"\n",
89+
"wsf_r = rasterio.open(wsf_file)\n",
90+
"\"\"\"\n",
91+
"print(\"FUBAR\")"
92+
]
93+
},
94+
{
95+
"cell_type": "code",
96+
"execution_count": 7,
97+
"metadata": {},
98+
"outputs": [],
99+
"source": [
100+
"#Clip GHSL using local files\n",
101+
"local_version = r\"J:\\Data\\GLOBAL\\GHSL\\Built\\GHS_BUILT_S_E2020_GLOBE_R2023A_54009_100_V1_0.tif\"\n",
102+
"if not os.path.exists(ghsl_file):\n",
103+
" ghsl_raster = rasterio.open(local_version)\n",
104+
" data, profile = rMisc.clipRaster(ghsl_raster, country)\n",
105+
" with rasterio.open(ghsl_file, 'w', **profile) as dst:\n",
106+
" dst.write(data)\n",
107+
"ghsl_r = rasterio.open(ghsl_file)"
108+
]
109+
},
110+
{
111+
"cell_type": "code",
112+
"execution_count": 8,
113+
"metadata": {},
114+
"outputs": [
115+
{
116+
"data": {
117+
"text/html": [
118+
"<div>\n",
119+
"<style scoped>\n",
120+
" .dataframe tbody tr th:only-of-type {\n",
121+
" vertical-align: middle;\n",
122+
" }\n",
123+
"\n",
124+
" .dataframe tbody tr th {\n",
125+
" vertical-align: top;\n",
126+
" }\n",
127+
"\n",
128+
" .dataframe thead th {\n",
129+
" text-align: right;\n",
130+
" }\n",
131+
"</style>\n",
132+
"<table border=\"1\" class=\"dataframe\">\n",
133+
" <thead>\n",
134+
" <tr style=\"text-align: right;\">\n",
135+
" <th></th>\n",
136+
" <th>class</th>\n",
137+
" <th>height</th>\n",
138+
" <th>wkt</th>\n",
139+
" <th>geometry</th>\n",
140+
" </tr>\n",
141+
" </thead>\n",
142+
" <tbody>\n",
143+
" <tr>\n",
144+
" <th>0</th>\n",
145+
" <td>NaN</td>\n",
146+
" <td>NaN</td>\n",
147+
" <td>POLYGON ((103.3991092 13.6076154, 103.3991148 ...</td>\n",
148+
" <td>POLYGON ((103.39911 13.60762, 103.39911 13.607...</td>\n",
149+
" </tr>\n",
150+
" <tr>\n",
151+
" <th>1</th>\n",
152+
" <td>NaN</td>\n",
153+
" <td>NaN</td>\n",
154+
" <td>POLYGON ((103.819971 13.2832912, 103.8199974 1...</td>\n",
155+
" <td>POLYGON ((103.81997 13.28329, 103.82000 13.283...</td>\n",
156+
" </tr>\n",
157+
" <tr>\n",
158+
" <th>2</th>\n",
159+
" <td>NaN</td>\n",
160+
" <td>NaN</td>\n",
161+
" <td>POLYGON ((103.8152555 13.2890315, 103.8151626 ...</td>\n",
162+
" <td>POLYGON ((103.81526 13.28903, 103.81516 13.289...</td>\n",
163+
" </tr>\n",
164+
" <tr>\n",
165+
" <th>3</th>\n",
166+
" <td>NaN</td>\n",
167+
" <td>NaN</td>\n",
168+
" <td>POLYGON ((105.5873344 12.3655821, 105.5873972 ...</td>\n",
169+
" <td>POLYGON ((105.58733 12.36558, 105.58740 12.365...</td>\n",
170+
" </tr>\n",
171+
" <tr>\n",
172+
" <th>4</th>\n",
173+
" <td>NaN</td>\n",
174+
" <td>NaN</td>\n",
175+
" <td>POLYGON ((103.0758412 13.2660819, 103.0758304 ...</td>\n",
176+
" <td>POLYGON ((103.07584 13.26608, 103.07583 13.266...</td>\n",
177+
" </tr>\n",
178+
" </tbody>\n",
179+
"</table>\n",
180+
"</div>"
181+
],
182+
"text/plain": [
183+
" class height wkt \\\n",
184+
"0 NaN NaN POLYGON ((103.3991092 13.6076154, 103.3991148 ... \n",
185+
"1 NaN NaN POLYGON ((103.819971 13.2832912, 103.8199974 1... \n",
186+
"2 NaN NaN POLYGON ((103.8152555 13.2890315, 103.8151626 ... \n",
187+
"3 NaN NaN POLYGON ((105.5873344 12.3655821, 105.5873972 ... \n",
188+
"4 NaN NaN POLYGON ((103.0758412 13.2660819, 103.0758304 ... \n",
189+
"\n",
190+
" geometry \n",
191+
"0 POLYGON ((103.39911 13.60762, 103.39911 13.607... \n",
192+
"1 POLYGON ((103.81997 13.28329, 103.82000 13.283... \n",
193+
"2 POLYGON ((103.81526 13.28903, 103.81516 13.289... \n",
194+
"3 POLYGON ((105.58733 12.36558, 105.58740 12.365... \n",
195+
"4 POLYGON ((103.07584 13.26608, 103.07583 13.266... "
196+
]
197+
},
198+
"execution_count": 8,
199+
"metadata": {},
200+
"output_type": "execute_result"
201+
}
202+
],
203+
"source": [
204+
"# read in and process Overture buildings\n",
205+
"ob = pd.read_csv(overture_buildings)\n",
206+
"ob_geoms = ob['wkt'].apply(loads)\n",
207+
"inB = gpd.GeoDataFrame(ob, geometry=ob_geoms, crs=4326)\n",
208+
"inB.head()\n"
209+
]
210+
},
211+
{
212+
"cell_type": "code",
213+
"execution_count": 9,
214+
"metadata": {},
215+
"outputs": [],
216+
"source": [
217+
"# attempt to rasterrize the buildings as polygons\n",
218+
"if not os.path.exists(overture_raster):\n",
219+
" rasterized_buildings = rMisc.rasterizeDataFrame(inB, templateRaster=ghsl_file, mergeAlg=\"ADD\", re_proj=True, nodata=0.)\n",
220+
" with rasterio.open(overture_raster, 'w', **rasterized_buildings['meta']) as dst:\n",
221+
" dst.write_band(1, rasterized_buildings['vals'])\n",
222+
"overture_r = rasterio.open(overture_raster)"
223+
]
224+
},
225+
{
226+
"cell_type": "code",
227+
"execution_count": 10,
228+
"metadata": {},
229+
"outputs": [],
230+
"source": [
231+
"# attempt to rasterrize the buildings as points\n",
232+
"if not os.path.exists(overture_raster_points):\n",
233+
" inB_points = inB.copy()\n",
234+
" inB_points['geometry'] = inB_points['geometry'].centroid\n",
235+
" rasterized_buildings = rMisc.rasterizeDataFrame(inB_points, templateRaster=ghsl_file, mergeAlg=\"ADD\", re_proj=True, nodata=0.)\n",
236+
" with rasterio.open(overture_raster_points, 'w', **rasterized_buildings['meta']) as dst:\n",
237+
" dst.write_band(1, rasterized_buildings['vals'])\n",
238+
"overture_r_points = rasterio.open(overture_raster_points) "
239+
]
240+
},
241+
{
242+
"cell_type": "code",
243+
"execution_count": 32,
244+
"metadata": {},
245+
"outputs": [],
246+
"source": [
247+
"# Compare rasterized buildings with built area dataset\n",
248+
"### Open both datasets and threshold them to get built area data\n",
249+
"o_thresh = 1\n",
250+
"ghsl_thresh = 3000\n",
251+
"\n",
252+
"o_data = overture_r_points.read(1)\n",
253+
"o_data = (o_data > o_thresh).astype('uint8')\n",
254+
"\n",
255+
"ghsl_data = ghsl_r.read(1)\n",
256+
"ghsl_data = (ghsl_data > ghsl_thresh).astype('uint8') * 10\n",
257+
"\n",
258+
"combo_data = o_data + ghsl_data\n",
259+
"\n",
260+
"# Write out the combined data\n",
261+
"out_file = os.path.join(out_folder, \"overture_vs_ghsl.tif\")\n",
262+
"if not os.path.exists(out_file):\n",
263+
" meta = overture_r_points.meta.copy()\n",
264+
" meta.update(dtype=rasterio.uint8, nodata=0)\n",
265+
" with rasterio.open(out_file, 'w', **meta) as out_raster:\n",
266+
" out_raster.write_band(1, combo_data)\n"
267+
]
268+
},
269+
{
270+
"cell_type": "code",
271+
"execution_count": null,
272+
"metadata": {},
273+
"outputs": [],
274+
"source": []
275+
}
276+
],
277+
"metadata": {
278+
"kernelspec": {
279+
"display_name": "urban_test",
280+
"language": "python",
281+
"name": "python3"
282+
},
283+
"language_info": {
284+
"codemirror_mode": {
285+
"name": "ipython",
286+
"version": 3
287+
},
288+
"file_extension": ".py",
289+
"mimetype": "text/x-python",
290+
"name": "python",
291+
"nbconvert_exporter": "python",
292+
"pygments_lexer": "ipython3",
293+
"version": "3.12.2"
294+
}
295+
},
296+
"nbformat": 4,
297+
"nbformat_minor": 2
298+
}

0 commit comments

Comments
 (0)