<jupyter_start><jupyter_text>Hack the Continent Open Buildings Challenge A study on the use of google open building data to extrapolate Gauteng's quality of life data to the rest of South Africa.This notebook is a study of how to extrapolate quality of living survey (QOLS) data using Google open buildings (GOB) data and the relative wealth index (RWI) produced from data made available by Meta. Google used image recognition on satellite images to produce the open building data set.I investigate how the average building size, as well as Google's trust in the data, correlated to different fields in the QOLS data set.In the learning and forecasting model, I also included a relative wealth index that was calculated previously for the entire country.  The quality of life survey only covers Gauteng, the model is thus trained using Gauteng data and then used to predict values for the rest of South Africa.The RWI data is calculated for a grid size of 2.4km. The GOB data has a latitude and longitude for each building detected.As the QOLS data had an aim of gathering a minimum amount of data points per ward, it was deemed that using wards as the granularity of calculation was the most practical approach to split the area. Descriptive statistics were calculated for each of the data set, per ward, and then the data sets were merged to generate a learning data set with ward level granularity.      The data sets used are:  QOLS - Quality of Life Survey (https://www.gcro.ac.za/research/project/detail/quality-of-life-survey-v-201718/)  GOB - Google Open Buildings](https://sites.research.google/open-buildings/)  RWI - Relative Wealth Index - Microestimates of wealth for all low- and middle-income countries Guanghua Chi, <NAME>, Sourav Chatterjee, <NAME>stock Proceedings of the National Academy of Sciences Jan 2022, 119 (3) e2113658119; DOI: 10.1073/pnas.2113658119 (https://data.humdata.org/dataset/76f2a2ea-ba50-40f5-b79c-db95d668b843/resource/66567fcf-98a0-4246-a799-6d4616633228/download/zaf_relative_wealth_index.csv)   Municipal ward boundaries: https://dataportal-mdb-sa.opendata.arcgis.com/datasets/e0223a825ea2481fa72220ad3204276b/about  This whole project will be made available on github.<jupyter_code># Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License<jupyter_output><empty_output><jupyter_text>I have downloaded the data sets to my personal Google drive for quicker access. The below simply connects to my google drive, list the files, and creates a variable to represent the path to input and output data.<jupyter_code>from google.colab import drive
drive.mount('/content/drive')
#gauteng = gpd.read_file()
!ls /content/drive
import os
rawdatadir = "/content/drive/MyDrive/Colab Notebooks/ZindiHackAfrica/01_raw_data"
preprocesseddir = "/content/drive/MyDrive/Colab Notebooks/ZindiHackAfrica/03_preprocessed_data/"
os.listdir(rawdatadir)<jupyter_output><empty_output><jupyter_text>Install and import the required libraries.<jupyter_code>!pip install s2geometry pygeos geopandas

import functools
import glob
import gzip
import multiprocessing
import os
import shutil
import tempfile
from typing import List, Optional, Tuple


import gdal
import geopandas as gpd
from google.colab import files
from IPython import display
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import s2geometry as s2
import shapely
import tensorflow as tf
import tqdm.notebook

import xgboost as xgb
from sklearn.metrics import mean_squared_error
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
### I want to be able to scroll across all the field, particularly in the qols data set so I remove the limit on the number of columns displayed.
pd.options.display.max_columns = None<jupyter_output><empty_output><jupyter_text>Preprocessing opening the data sets for pre-processing, The raw GOB data is not stored locally, so not in the list below but we'll get to it further along.<jupyter_code>#rsa = gpd.read_file(rawdatadir + '/stanford-nt680nw4634-geojson.json')  # geometries for south African municipal areas... seems kinda old, eg reference to "northern province"
rwi = gpd.read_file(rawdatadir + '/zaf_relative_wealth_index.csv')  # relative wealth index
#allpoi = gpd.read_file(f'zip://{rawdatadir}/hotosm_zaf_points_of_interest_points_shp.zip')  # points of interest schools etc
#financialpoi = gpd.read_file(f'zip://{rawdatadir}/hotosm_zaf_financial_services_points_shp.zip') # points of interest financial 
qols = gpd.read_file(f'{rawdatadir}/qols-v-2017-2018-v1.1.csv') # quality of life standards from 2017 - 2018
wards = gpd.read_file(f'zip://{rawdatadir}/SA_Wards2020.zip') # South African Ward boudaries<jupyter_output><empty_output><jupyter_text>Have a look at the data gathered in the QOLS survey, its pretty awesome, there's a lot here that could be invetigated.<jupyter_code>qols.head(3)<jupyter_output><empty_output><jupyter_text>Just a sanity check that the QOLS data has the same number of wards for Gauteng as the wards data set.  Interesting to note that there are 529 wards in our training set that more than I expected, and it give a better training set size are more granularity than I thought. Nice!<jupyter_code>qolsnumwards = len(qols['ward'].unique())
numwards = len(wards[wards['Province'] == 'Gauteng']['WardID'].unique())
print(f'number of wards in the quality of life survey data set: {qolsnumwards}')
print(f'number of wards according to arcgis 2020: {numwards}')<jupyter_output><empty_output><jupyter_text>I chose two fields to focus on for this work. One that I expect should correlate well to the GOB data and one that I expect would not correlate well. They are the number of rooms in the household and the length of stay in the current household.  below I check the unique values in each of the two files.<jupyter_code>print('Values in the rooms field: ', (qols['Q1_02_rooms'].unique()))
print('Values in the length of stay field: ', qols['Q4_04_lengthof_stay'].unique())<jupyter_output><empty_output><jupyter_text>The number of rooms is numeric except the 10+, it is thus replaced by '10' and then cast to an int.  The length of stay cannot be cast to a numeric value easily and thus is mapped with a dictionary lookup.<jupyter_code>maplengthof_stay = {'More than 10 years' : 10, '3-4 years': 4, "I've always lived here": 25,
       'Less than 1 year':1, '1-2 years':2, '5-10 years':8}

qolstoy = qols[['ward', 'Q4_04_lengthof_stay', 'Q1_02_rooms']].copy()
qolstoy.replace(maplengthof_stay, inplace = True)
qolstoy.replace({'10+':'10'}, inplace = True)
qolstoy['Q1_02_rooms']=qolstoy['Q1_02_rooms'].astype(int)

qolstoy = qolstoy.groupby('ward').agg(['median','min', 'max', 'std'])
wards.set_geometry('geometry', drop=True, inplace=True, crs='EPSG:4326')
wards<jupyter_output><empty_output><jupyter_text>The RWI data needs to be reformatted and aggregated per ward.<jupyter_code>rwi['geometry'] =gpd.points_from_xy(rwi['longitude'], rwi['latitude'])
rwi = gpd.GeoDataFrame(rwi, crs="EPSG:4326")

points = wards.sjoin(rwi)

rwiward = points[['WardID','rwi','error']].groupby('WardID').agg(['median','min', 'max', 'std'])
rwiward<jupyter_output><empty_output><jupyter_text>the RWI data can then be joined onto the QOLS data set per ward.<jupyter_code>rwiqolward = rwiward.merge(qolstoy, left_index = True, right_index = True, how = 'left')<jupyter_output><empty_output><jupyter_text>Plotting the boundaries of the wards data. just to get a feel for what they look like. We can see that wards seem to be smaller where there are more dense populations, also there are a lot of wards in the South Africa.<jupyter_code>wards.boundary.plot()<jupyter_output><empty_output><jupyter_text>Below is the code that is used to get the data from Googles Open building data, I used it as supplied in the example notebook with some minor changes.  Download buildings data for a region in Africa [takes up to 15 minutes for large countries]<jupyter_code>#@markdown Select a region from either the [Natural Earth low res](https://www.naturalearthdata.com/downloads/110m-cultural-vectors/110m-admin-0-countries/) (fastest), [Natural Earth high res](https://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-0-countries/) or [World Bank high res](https://datacatalog.worldbank.org/dataset/world-bank-official-boundaries) shapefiles:
region_border_source = 'Natural Earth (Low Res 110m)'  #@param ["Natural Earth (Low Res 110m)", "Natural Earth (High Res 10m)", "World Bank (High Res 10m)"]
region = 'ZAF (South Africa)'  #@param ["", "AGO (Angola)", "BDI (Burundi)", "BEN (Benin)", "BFA (Burkina Faso)", "BWA (Botswana)", "CAF (Central African Republic)", "CIV (Côte d'Ivoire)", "COD (Democratic Republic of the Congo)", "COG (Republic of the Congo)", "DJI (Djibouti)", "DZA (Algeria)", "EGY (Egypt)", "ERI (Eritrea)", "ETH (Ethiopia)", "GAB (Gabon)", "GHA (Ghana)", "GIN (Guinea)", "GMB (The Gambia)", "GNB (Guinea-Bissau)", "GNQ (Equatorial Guinea)", "KEN (Kenya)", "LBR (Liberia)", "LSO (Lesotho)", "MDG (Madagascar)", "MOZ (Mozambique)", "MRT (Mauritania)", "MWI (Malawi)", "NAM (Namibia)", "NER (Niger)", "NGA (Nigeria)", "RWA (Rwanda)", "SDN (Sudan)", "SEN (Senegal)", "SLE (Sierra Leone)", "SOM (Somalia)", "SWZ (eSwatini)", "TGO (Togo)", "TUN (Tunisia)", "TZA (Tanzania)", "UGA (Uganda)", "ZAF (South Africa)", "ZMB (Zambia)", "ZWE (Zimbabwe)"]
province = "" #@param ["NORTHERN CAPE", "NORTHERN PROVINCE", "MPUMALANGA", "NORTH WEST", "GAUTENG", "FREE STATE", "KWAZULU-NATAL", "EASTERN CAPE", "WESTERN CAPE", ""]

# @markdown Alternatively, specify an area of interest in [WKT format](https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry) (assumes crs='EPSG:4326'); this [tool](https://arthur-e.github.io/Wicket/sandbox-gmaps3.html) might be useful.
your_own_wkt_polygon = ''  #@param {type:"string"}


BUILDING_DOWNLOAD_PATH = ('gs://open-buildings-data/v1/'
                          'polygons_s2_level_6_gzip_no_header')

def get_filename_and_region_dataframe(
    region_border_source: str, region: str,
    your_own_wkt_polygon: str) -> Tuple[str, gpd.geodataframe.GeoDataFrame]:
  """Returns output filename and a geopandas dataframe with one region row."""

  if your_own_wkt_polygon:
    filename = 'open_buildings_v1_polygons_your_own_wkt_polygon.csv.gz'
    region_df = gpd.GeoDataFrame(
        geometry=gpd.GeoSeries.from_wkt([your_own_wkt_polygon]),
        crs='EPSG:4326')
    if not isinstance(region_df.iloc[0].geometry,
                      shapely.geometry.polygon.Polygon) and not isinstance(
                          region_df.iloc[0].geometry,
                          shapely.geometry.multipolygon.MultiPolygon):
      raise ValueError("`your_own_wkt_polygon` must be a POLYGON or "
                      "MULTIPOLYGON.")
  elif province:
    region_df = gpd.GeoDataFrame(geometry=[rsa[(rsa['nam'] == province)].geometry.unary_union], crs=rsa.crs)
    filename = f'open_buildings_v1_polygons_{province}.csv.gz'
    jhbpoly.plot()

    print(f'Preparing your_own_wkt_polygon.')
    return filename, region_df

  if not region:
    raise ValueError('Please select a region or set your_own_wkt_polygon.')

  if region_border_source == 'Natural Earth (Low Res 110m)':
    url = ('https://www.naturalearthdata.com/http//www.naturalearthdata.com/'
           'download/110m/cultural/ne_110m_admin_0_countries.zip')
    !wget -N {url}
    display.clear_output()
    region_shapefile_path = os.path.basename(url)
    source_name = 'ne_110m'
  elif region_border_source == 'Natural Earth (High Res 10m)':
    url = ('https://www.naturalearthdata.com/http//www.naturalearthdata.com/'
           'download/10m/cultural/ne_10m_admin_0_countries.zip')
    !wget -N {url}
    display.clear_output()
    region_shapefile_path = os.path.basename(url)
    source_name = 'ne_10m'
  elif region_border_source == 'World Bank (High Res 10m)':
    url = ('https://development-data-hub-s3-public.s3.amazonaws.com/ddhfiles/'
           '779551/wb_countries_admin0_10m.zip')
    !wget -N {url}
    !unzip -o {os.path.basename(url)}
    display.clear_output()
    region_shapefile_path = 'WB_countries_Admin0_10m'
    source_name = 'wb_10m'

  region_iso_a3 = region.split(' ')[0]
  filename = f'open_buildings_v1_polygons_{source_name}_{region_iso_a3}.csv.gz'
  region_df = gpd.read_file(region_shapefile_path).query(
      f'ISO_A3 == "{region_iso_a3}"').dissolve(by='ISO_A3')[['geometry']]
  print(f'Preparing {region} from {region_border_source}.')
  return filename, region_df


def get_bounding_box_s2_covering_tokens(
    region_geometry: shapely.geometry.base.BaseGeometry) -> List[str]:
  region_bounds = region_geometry.bounds
  s2_lat_lng_rect = s2.S2LatLngRect_FromPointPair(
      s2.S2LatLng_FromDegrees(region_bounds[1], region_bounds[0]),
      s2.S2LatLng_FromDegrees(region_bounds[3], region_bounds[2]))
  coverer = s2.S2RegionCoverer()
  # NOTE: Should be kept in-sync with s2 level in BUILDING_DOWNLOAD_PATH.
  coverer.set_fixed_level(6)
  coverer.set_max_cells(1000000)
  return [cell.ToToken() for cell in coverer.GetCovering(s2_lat_lng_rect)]


def s2_token_to_shapely_polygon(
    s2_token: str) -> shapely.geometry.polygon.Polygon:
  s2_cell = s2.S2Cell(s2.S2CellId_FromToken(s2_token, len(s2_token)))
  coords = []
  for i in range(4):
    s2_lat_lng = s2.S2LatLng(s2_cell.GetVertex(i))
    coords.append((s2_lat_lng.lng().degrees(), s2_lat_lng.lat().degrees()))
  return shapely.geometry.Polygon(coords)


def download_s2_token(
    s2_token: str, region_df: gpd.geodataframe.GeoDataFrame) -> Optional[str]:
  """Downloads the matching CSV file with polygons for the `s2_token`.

  NOTE: Only polygons inside the region are kept.
  NOTE: Passing output via a temporary file to reduce memory usage.

  Args:
    s2_token: S2 token for which to download the CSV file with building
      polygons. The S2 token should be at the same level as the files in
      BUILDING_DOWNLOAD_PATH.
    region_df: A geopandas dataframe with only one row that contains the region
      for which to keep polygons.

  Returns:
    Either filepath which contains a gzipped CSV without header for the
    `s2_token` subfiltered to only contain building polygons inside the region
    or None which means that there were no polygons inside the region for this
    `s2_token`.
  """
  s2_cell_geometry = s2_token_to_shapely_polygon(s2_token)
  region_geometry = region_df.iloc[0].geometry
  prepared_region_geometry = shapely.prepared.prep(region_geometry)
  # If the s2 cell doesn't intersect the country geometry at all then we can
  # know that all rows would be dropped so instead we can just return early.
  if not prepared_region_geometry.intersects(s2_cell_geometry):
    return None
  try:
    # Using tf.io.gfile.GFile gives better performance than passing the GCS path
    # directly to pd.read_csv.
    with tf.io.gfile.GFile(
        os.path.join(BUILDING_DOWNLOAD_PATH, f'{s2_token}_buildings.csv.gz'),
        'rb') as gf:
      # If the s2 cell is fully covered by country geometry then can skip
      # filtering as we need all rows.
      if prepared_region_geometry.covers(s2_cell_geometry):
        with tempfile.NamedTemporaryFile(mode='w+b', delete=False) as tmp_f:
          shutil.copyfileobj(gf, tmp_f)
          return tmp_f.name
      # Else take the slow path.
      # NOTE: We read in chunks to save memory.
      csv_chunks = pd.read_csv(
          gf, chunksize=2000000, dtype=object, compression='gzip', header=None)
      tmp_f = tempfile.NamedTemporaryFile(mode='w+b', delete=False)
      tmp_f.close()
      for csv_chunk in csv_chunks:
        points = gpd.GeoDataFrame(
            geometry=gpd.points_from_xy(csv_chunk[1], csv_chunk[0]),
            crs='EPSG:4326')
        # sjoin 'within' was faster than using shapely's 'within' directly.
        points = gpd.sjoin(points, region_df, predicate='within')
        csv_chunk = csv_chunk.iloc[points.index]
        csv_chunk.to_csv(
            tmp_f.name,
            mode='ab',
            index=False,
            header=False,
            compression={
                'method': 'gzip',
                'compresslevel': 1
            })
      return tmp_f.name
  except tf.errors.NotFoundError:
    return None

# Clear output after pip install.
display.clear_output()
filename, region_df = get_filename_and_region_dataframe(region_border_source,
                                                        region,
                                                        your_own_wkt_polygon)
# Remove any old outputs to not run out of disk.
for f in glob.glob('/tmp/open_buildings_*'):
  os.remove(f)
# Write header to the compressed CSV file.
with gzip.open(f'/tmp/{filename}', 'wt') as merged:
  merged.write(','.join([
      'latitude', 'longitude', 'area_in_meters', 'confidence', 'geometry',
      'full_plus_code'
  ]) + '\n')
download_s2_token_fn = functools.partial(download_s2_token, region_df=region_df)
s2_tokens = get_bounding_box_s2_covering_tokens(region_df.iloc[0].geometry)
# Downloads CSV files for relevant S2 tokens and after filtering appends them
# to the compressed output CSV file. Relies on the fact that concatenating
# gzipped files produces a valid gzip file.
# NOTE: Uses a pool to speed up output preparation.
with open(f'/tmp/{filename}', 'ab') as merged:
  with multiprocessing.Pool(4) as e:
    for fname in tqdm.notebook.tqdm(
        e.imap_unordered(download_s2_token_fn, s2_tokens),
        total=len(s2_tokens)):
      if fname:
        with open(fname, 'rb') as tmp_f:
          shutil.copyfileobj(tmp_f, merged)
        os.unlink(fname)<jupyter_output><empty_output><jupyter_text>Visualise the dataFirst we convert the CSV file into a GeoDataFrame. The CSV files can be quite large because they include the polygon outline of every building. For this example we only need longitude and latitude, so we only process those columns to save memory.<jupyter_code>buildings = pd.read_csv(
    f"/tmp/{filename}", engine="c",
    usecols=['latitude', 'longitude', 'area_in_meters', 'confidence'])

print(f"Read {len(buildings):,} records.")<jupyter_output><empty_output><jupyter_text>Here I write the GOB data and the rwi+qols data to file so they can be opened more quickly in future. They are commented out for now as I have the data on file and I do not want to accidentally overwrite it.<jupyter_code>#rwiqolward.to_csv(preprocesseddir + 'rwiqolward.csv')  # geometries for south African municipal areas... seems kinda old, eg reference to "northern province"
#buildings.to_csv(preprocesseddir + 'buildings.csv')  # relative wealth index<jupyter_output><empty_output><jupyter_text>Open preprocessed data and do learning<jupyter_code>rwiqolward = pd.read_csv(preprocesseddir + 'rwiqolward.csv', header=[0,1,2])<jupyter_output><empty_output><jupyter_text>The Buildings data set is really big, and trying to aggregate it per ward means my free colab machine runs out of ram and falls over.  Here I find the extremes of the latitude and longitude for each ward and filter the building data on that before filtering it on the geometry of the ward.I then aggregate the building data for the ward.  This is pretty time consuming and I do not want to do it every time so I write the output to file.<jupyter_code># get the descriptive statistics for building area per ward.
buildings = pd.read_csv(preprocesseddir + 'buildings.csv')
outlist = []
for wardno in wards['WardID'].unique():
  eachward = wards[wards['WardID'] == wardno]
  minx,miny,maxx,maxy = eachward.bounds.iloc[0]
  buildingswards = buildings[(buildings['longitude'] < maxx) &(buildings['longitude'] > minx) &(buildings['latitude'] > miny) &(buildings['latitude'] < maxy)]
  buildingswards = gpd.GeoDataFrame(buildingswards,
        geometry=gpd.points_from_xy(buildingswards.longitude,
                                    buildingswards.latitude),crs = "EPSG:4326")
  buildingswards = gpd.sjoin(eachward, buildingswards, how="left")
  bwaggs = buildingswards[['area_in_meters',	'confidence']].agg(['median','min','max','std','count'])
  dfout = pd.DataFrame([{'wardID': wardno}])
  dfout[['Amedian','Amin','Amax','Astd','Acount','Cmedian','Cmin','Cmax','Cstd','Ccount']] = bwaggs['area_in_meters'].append(bwaggs['confidence'])#.append(wardno)
  outlist.append(dfout)
dfdone = pd.concat(outlist)
dfdone.to_csv(f'{preprocesseddir}dfbuildingstats.csv')

dfdone= pd.read_csv(f'{preprocesseddir}dfbuildingstats.csv')
dfdone['wardID'] = dfdone['wardID'].astype(int)
print(rwiqolward.shape)
rwiqolward.head()<jupyter_output><empty_output><jupyter_text>From the above we can see that the headings from my rwiqolsward data set are a mess so I cleaned them up a little.<jupyter_code>cols = ['wardID'] + [' '.join(col).strip().split(' Unnamed:')[0] for col in rwiqolward.columns.values][1:]
rwiqolward.columns = cols<jupyter_output><empty_output><jupyter_text>Merging the rwiqolsward data with the building data(dfdone) gets me the cleaned data set that I want to use for learning and prediction. Model1: Learn the number of rooms in a house hold.<jupyter_code>dfrooms = dfdone.merge(rwiqolward[['wardID','rwi median','rwi min','rwi max','rwi std','error median','error min','error max','error std','Q1_02_rooms median']], on = 'wardID', how = 'left').copy()
dfrooms.head()<jupyter_output><empty_output><jupyter_text>Since this data set is for all of South Africa, but I only have qols data for Gauteng, I split this data in two: the learning data for Gauteng, and the predicting data for the rest of the country.<jupyter_code>dfqolspred = dfrooms[dfrooms['Q1_02_rooms median'].isna()].copy()
dfqolslearn = dfrooms[~dfrooms['Q1_02_rooms median'].isna()]<jupyter_output><empty_output><jupyter_text>Set up the data for a little machine learning.<jupyter_code>data = dfqolslearn#.values
# split data into input and output columns
X, y = data.iloc[:, :-1], data.iloc[:, -1]
data_dmatrix = xgb.DMatrix(data=X,label=y)

#from sklearn.model_selection import train_test_split
#
#X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=123)<jupyter_output><empty_output><jupyter_text>Do a quick cross validation XGboost training exercise to see how well the data can be learned, we can see how many rounds are needed to have good learning. And give us a good idea if the model is fitting well or if its overfitting.This graph shows that the model learns well on this data the mean root mean square error is less that one, and the test and train plots track well together, we can see some over fitting as the rounds go on, buy at around 40 rounds the fit looks good and the over fitting looks minimal.<jupyter_code>params = {"objective":"reg:squarederror",'colsample_bytree': 0.1,'learning_rate': 0.08,
                'max_depth': 8, 'alpha': 8}
cv_results = xgb.cv(dtrain=data_dmatrix, params=params, nfold=3,
                    num_boost_round=100,early_stopping_rounds=10,metrics="rmse", as_pandas=True, seed=123)
cv_results[['train-rmse-mean','test-rmse-mean']].plot()<jupyter_output><empty_output><jupyter_text>Cross validation does not produce a model though so we need to train the model on our training data, we can train it on all the training data, as we are happy with the parameters from the cross validation output.<jupyter_code>xg_reg = xgb.train(params=params, dtrain=data_dmatrix, num_boost_round=40)<jupyter_output><empty_output><jupyter_text>XGboots has a cool functions that shows us the importance of the most important features in the data set. Here we see that the standard deviation on the area of the building is the most important feature to the model, this is unexpected as I would have expected the median area to be the most useful. I am glad I checked, because this could mean that the learning was not a good as I thought when I looked at the test/train plot above.<jupyter_code>xgb.plot_importance(xg_reg)
plt.rcParams['figure.figsize'] = [5, 5]
plt.show()<jupyter_output><empty_output><jupyter_text>Just for completeness I then use this model to predict the median number of rooms for every ward across the country.<jupyter_code>pred_dmatrix = xgb.DMatrix(data=dfqolspred.iloc[:, :-1])
preds = xg_reg.predict(pred_dmatrix)
dfqolspred.loc[:,'Q1_02_rooms median'] = preds
dfqolspred.head()<jupyter_output><empty_output><jupyter_text>Model2: Learn the length of stay.<jupyter_code>rwiqolward.head()
dflos = dfdone.merge(rwiqolward[['wardID','rwi median','rwi min','rwi max','rwi std','error median','error min','error max','error std','Q4_04_lengthof_stay median']], on = 'wardID', how = 'left')
dflos.drop('wardID',axis = 1, inplace = True)<jupyter_output><empty_output><jupyter_text>Since this data set is for all of South Africa, but I only have qols data for Gauteng, I split this data in two: the learning data for Gauteng, and the predicting data for the rest of the country.<jupyter_code>dfqolspred = dflos[dflos['Q4_04_lengthof_stay median'].isna()]
dfqolslearn = dflos[~dflos['Q4_04_lengthof_stay median'].isna()]<jupyter_output><empty_output><jupyter_text>Format the data for consumption by xgboost library<jupyter_code>data = dfqolslearn#.values
# split data into input and output columns
X, y = data.iloc[:, :-1], data.iloc[:, -1]
data_dmatrix = xgb.DMatrix(data=X,label=y)<jupyter_output><empty_output><jupyter_text>The train and test plot below shows that after about 10 rounds the model starts to over fit very badly, at the point of divergence of the plots the Root Mean Square Error is almost 5. Which is pretty bad considering our input range is 2-25.<jupyter_code>params = {"objective":"reg:squarederror",'colsample_bytree': 0.1,'learning_rate': 0.08,
                'max_depth': 8, 'alpha': 8}
cv_results = xgb.cv(dtrain=data_dmatrix, params=params, nfold=3,
                    num_boost_round=100,early_stopping_rounds=10,metrics="rmse", as_pandas=True, seed=123)
cv_results[['train-rmse-mean','test-rmse-mean']].plot()<jupyter_output><empty_output><jupyter_text>I keep the same parameters used in the cross validation but limit the number of rounds to 12 for the training, to avoid over fitting.<jupyter_code>xg_reg = xgb.train(params=params, dtrain=data_dmatrix, num_boost_round=12)<jupyter_output><empty_output><jupyter_text>Having a look at the feature importance we see that the RWI max is the most important for this model, I can understand that there is a correlation between the wealth of an area and how long people have lived there. There is some learning happening here but I am not confident that the predictions created by this model would actually be good enough to be valuable, so I am not going to extrapolate this for the rest of South Africa.<jupyter_code>xgb.plot_importance(xg_reg)
plt.rcParams['figure.figsize'] = [5, 5]
plt.show()<jupyter_output><empty_output><jupyter_text>Conclusion The QOLS Data is awesome, its a good sample size ~23000 and covers a lot of detail. I have just scratched the surface on what is in it and there is a lot of scope to improve on what I have done here. One could look at other fields on the QOLS data, bring in other data sets to enrich the learning data set, and try different models.  In the set up for this project I have found many other data sets that may be useful to expand this project, but in the spirit of the Zindi Competition I wanted to use the GOB data and one other as a counter example.my list of other resource will be uploaded to my github in due course https://github.com/g-webberThanks for reading.<jupyter_code><jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10195924485207523, 'fraction_numerical': 0.01705863459548196, 'mean_word_length': 4.869710934227063, 'pattern_counts': {'":': 3, '<': 129, '<?xml version=': 0, '>': 133, 'https://': 14, 'lorem ipsum': 0, 'www.': 8, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '31', 'dir': 'jupyter-structured-clean-dedup', 'id': '221676', 'license': 'MIT', 'n_tokens_mistral': 9475, 'n_tokens_neox': 8985, 'n_words': 2862, 'path': 'Hack the Continent Open Buildings Challenge/Solution 3/Open_buildings_spatial_analysis_examples.ipynb', 'repo_name': 'ZindiAfrica/Data-collection-and-visualisation', 'size': '28033'}