Components assessment

[iferencik]

Built env. component currently consists of:

• buildings
  • source: vida google+ms+osm
  • dissaggregated by country
  • cloud source
  • variables: list of var, to do
  • asssess
    1. download (async)
    2. compute (reproject)
    3. no of buildings
       4. evaluate (converting building centroids into raster indices and count buildings)
• roads
• el. grid

[Thuhaa]

The Built Environment

Related to issue #18

Comprised of:

• Roads
• Electricity Grid
• Buildings

1. Roads

Download

The Roads dataset can be downloaded from Overture Maps
Implemetation of download need to consider the following:

1. Country Identifier - ISO3 code
2. Bounding Box - Download for a specific bounding box

Processing

Variables:

• Length of Roads Affected (Issue roads_length_affected #50 )
• Percentage of roads affected by specific crisis (length) dissagregated by divisions within a country (admin1, admin2... etc) (Issue percentage_roads_length_affected #51)

Steps

a). Identify the divisions within the country that have been affectred by a crisis.
b). Clip the roads to divisions within the country that have been affected by that crisis.
c). Create a mask of the affected "area" such that

True	False	False	True	True
True	False	True	True	False
True	True	True	True	True
True	True	True	False	False
False	True	True	True	True

Where the True would be the affected are

a.
To get the affected roads:
affected_roads = roads * affected_mask

This process require a rasterization of the roads to produce a raster of the roads footprint with burn value 1. And the resolution to be similar to the mask, and have the exact same extent so as to produce proper results when running the zonal stats.

The rasterization can be done by gdal.RasterizeLayer
You can get the total length of roads by running zonal stats with sum then multiplying the result by the given resolution in gdal.RasterizeLayer

💡 Note: This is pseudocode. Not actual implementation of api

length = ZonalStats(operation="sum", input_dataset="rasterized_roads.tif", polygon_dataset="kenya_admin1.fgb") * resolution

• Where the roads intersect the affected area to have the affected_roads
• Find out where the roads intersect the non_affected area to have non_affected_roads

affected_roads_percentage = affected_roads_length/total_length_of_roads

affected_roads_percentage would be a decimal value between 0 and 1 for each admin division

2. Electricity Grid

Download

The current known source of electricity grid data (February 14th 2025) is GridFinder
Similat to roads, the download would need to consider:

1. Country Identifier - ISO3 code
2. Bounding Box - Download for a specified bounding box

Processing

Variables

1. Percentage of the electricity grid that have been affected by the specific crisis. (Issue percentage_elgrid_length_affected #54 )
2. Length of the electricity grid affected (Issue electricity_grid_length #52 )
3. Total length of grid (Issue electricity_grid_length #52 )

Steps

a). Identify divisions within the country/bounding box that have been affected by the crisis.
b). Clip the electricity grid to those specific divisions.
c). Create a mask of the affected area as in the roads such that:

True	False	False	True	True
True	False	True	True	False
True	True	True	True	True
True	True	True	False	False
False	True	True	True	True

Where True would be the affected parts

To get the affected electricity grid:
affected_grid = grid * affected_mask

This process require a rasterization of the grid to produce a raster of the grid footprint with burn value 1. And the resolution to be similar to the mask, and have the exact same extent so as to produce proper results when running the zonal stats.

The rasterization can be done by gdal.RasterizeLayer
You can get the length of grid by running zonal stats with sum then multiplying the result by the given resolution in gdal.RasterizeLayer

💡 Note: This is pseudocode. Not actual implementation of api

length = ZonalStats(operation="sum", input_dataset="rasterized_grid.tif", polygon_dataset="kenya_admin1.fgb") * resolution

• where the electricity grid intersect the affected area - you have the affected grid
• Where you have the electricity grid intersect the non affected area, you have the non affected grid

affected_grid_percentage = affected_grid_length / total_grid_length

affected_grid_percentage would be a decimal value between 0 and 1.

3. Buildings

Variables:

1. Percentage of affected buildings (by area_in_square_kilometers of affected buildings) (Issue percentage_buildings_surface_area_affected #48 )
2. Total building area (Issue buildings_surface_area #46 )
3. Affected building area (Issue buildings_surface_area_affected #47 )
4. Total number of buildings (Issue number_of_buildings #43 )
5. Percentage of buildings Affected (by count of buildings) (Issue percentage_buildings_affected #45 )
6. Affected building count (Issue number_of_buildings_affected #44 )

Download buildings from vida which has buildings from google, buildings, and open street map

1. Provide bounding box or admin units for area of interest (AOI) to download buildings within that area.
   rapida tool currently has implementation of this in https://github.com/UNDP-Data/geo-cb-surge/blob/main/cbsurge/components/builtenv/buildings/fgb.py
2. The building dataset is massive - sometimes even for reasonably small sized administrative units/bounding boxes. This makes computation resource intensive. It is therefore imperative to make it as computationally easy as possible. For this reason. Get the centroids within each buiding polygon: Centroids are simple points that are located at each of the polygons center and would have the exact same attributes as the buildings, only that the geometries are simple enough to make it less computationally intensive.

💡 Note: This is an example of python function to generate the centroid dataset, and maintain all the attributes of the original dataset

from osgeo import ogr, osr
import os

def compute_centroids(input_file: str, output_file: str, driver_name: str = "FlatGeobuf"):
    """
    Computes centroids of polygons/multipolygons from an input vector dataset
    and saves them in a new dataset while preserving attributes.

    :param input_file: Path to the input vector file (.fgb, .gpkg, .shp, etc.).
    :param output_file: Path to the output file where centroids will be saved.
    :param driver_name: OGR driver name (default: "FlatGeobuf"). Can be "GPKG" or "ESRI Shapefile".
    """
    dataset = ogr.Open(input_file, 0)  # Read-only mode
    if not dataset:
        raise RuntimeError(f"Could not open {input_file}")
    input_layer = dataset.GetLayer()
    spatial_ref = input_layer.GetSpatialRef()  # Keep original spatial reference
    driver = ogr.GetDriverByName(driver_name)
    if not driver:
        raise RuntimeError(f"OGR driver '{driver_name}' not found")
    if os.path.exists(output_file):
        driver.DeleteDataSource(output_file)
    output_dataset = driver.CreateDataSource(output_file)
    output_layer = output_dataset.CreateLayer("centroids", srs=spatial_ref, geom_type=ogr.wkbPoint)
    input_layer_def = input_layer.GetLayerDefn()
    for i in range(input_layer_def.GetFieldCount()):
        output_layer.CreateField(input_layer_def.GetFieldDefn(i))
    for feature in input_layer:
        geom = feature.GetGeometryRef()
        if geom and geom.GetGeometryName() in ["POLYGON", "MULTIPOLYGON"]:
            centroid = geom.Centroid()
            new_feature = ogr.Feature(output_layer.GetLayerDefn())
            new_feature.SetGeometry(centroid)
            for i in range(feature.GetFieldCount()):
                new_feature.SetField(i, feature.GetField(i))
            output_layer.CreateFeature(new_feature)
            new_feature = None  # Cleanup

    # Close datasets
    dataset = None
    output_dataset = None

if __name__ == "__main__":
    compute_centroids("input_building_dataset.fgb", "buildings_centroids.fgb", driver_name="FlatGeobuf")

The following ogr2ogr has been previously used:

#!/bin/bash
ogr2ogr -sql "SELECT ST_Centroid(geometry) as geometry, fid, area_in_meters, confidence FROM bangladesh_buildings" -dialect sqlite bangladesh_buildings_centroids.fgb bangladesh_buildings.fgb -nln bangladesh_buildings_centroids -overwrite -progress 

8. The next step would be to polygonize the mask of the area affected

from osgeo import gdal, ogr, osr

def polygonize_raster(input_raster, output_shapefile):
    """
    Converts a raster file into a vector polygon layer.
    
    Parameters:
    - input_raster: str, path to the input raster file
    - output_shapefile: str, path to the output shapefile
    """
    dataset = gdal.Open(input_raster)
    if dataset is None:
        raise ValueError(f"Could not open raster file: {input_raster}")
    
    band = dataset.GetRasterBand(1)
    driver = ogr.GetDriverByName("ESRI Shapefile")
    out_ds = driver.CreateDataSource(output_shapefile)
    out_layer = out_ds.CreateLayer("polygonized", geom_type=ogr.wkbPolygon)
    field = ogr.FieldDefn("DN", ogr.OFTInteger)
    out_layer.CreateField(field)

    gdal.Polygonize(band, None, out_layer, 0, [], callback=None)

    # Cleanup
    dataset = None
    out_ds = None

polygonize_raster("affected_areas.tif", "affected_areas.gpkg")

9. compute intersection of the centroids and the affected area polygon. And with this produce counts of affected buildings.

ogr2ogr -f GPKG intersection_output.gpkg \
  -dialect SQLite \
  -sql "SELECT ST_Intersection(a.geometry, b.geometry) AS geometry, a.*, b.* FROM building_centroids a JOIN affected_areas b ON ST_Intersects(a.geometry, b.geometry)" \
  building_centroids.gpkg

10. Run the following computations.

percentage_of_affected_buildings_by_count = affected_buildings_count/total_buildings_in_admin

percentage_of_area_of_affected_buildings = sum(area_of_addected_buildings)/sum(total_area_covered_by_buildings)