Raster-based market access

import sys
import os
import rasterio

import geopandas as gpd
import skimage.graph as graph

import GOSTrocks.rasterMisc as rMisc
import GOSTrocks.mapMisc as mapMisc

sys.path.append("../../src")

import GOSTnetsraster.market_access as ma

%load_ext autoreload
%autoreload 2

iso3 = "SEN"

local_folder = f"C:/WBG/Work/Projects/HealthAccess/Tutorials/GNR/{iso3}"
results_folder = os.path.join(
    local_folder, "results"
)  # os.path.join(base_folder, iso3)
least_cost_tt_file = os.path.join(results_folder, "least_cost_travel_time.tif")

if not os.path.exists(results_folder):
    os.makedirs(results_folder)
tutorial_folder = "tutorial_data"

dests = os.path.join(tutorial_folder, f"{iso3}_cities.geojson")
friction_surface = os.path.join(tutorial_folder, f"{iso3}_friction_surface.tif")
# if you completed the friction surface generation step, you can also use:
friction_surface = os.path.join(local_folder, "friction_surface.tif")

inD = gpd.read_file(dests)
inR = rasterio.open(friction_surface)
inD = inD.to_crs(inR.crs)  # destinations and raster need to be in the same CRS

frictionD = inR.read()[0, :, :]
# convert friction surface to traversal time (lazily). Original data are
#    the original data are minutes to travel 1 m, so we will convert to
#    minutes to cross the cell
frictionD = frictionD * inR.res[0]  # cell size in meters
mcp = graph.MCP_Geometric(frictionD)

ma.calculate_travel_time?

travel_costs, traceback = ma.calculate_travel_time(inR, mcp, inD)
travel_costs = travel_costs.astype(inR.meta["dtype"])
with rasterio.open(least_cost_tt_file, "w", **inR.meta) as out_file:
    out_file.write_band(1, travel_costs)

mapMisc.static_map_raster(
    rasterio.open(least_cost_tt_file),
    thresh=[-1, 30, 60, 120, 180, 240, 10000],
    reverse_colormap=True,
)

ma.generate_feature_vectors?

drive_time_thresholds = [60, 120, 180, 240]  # minutes
drive_vectors = ma.generate_feature_vectors(inR, mcp, inD, drive_time_thresholds)
drive_vectors.to_file(os.path.join(results_folder, "drive_vectors.shp"))

drive_vectors.head()

mapMisc.static_map_vector(
    drive_vectors.loc[drive_vectors["threshold"] == 60], map_column="IDX"
)

ma.calculate_gravity?

# calculate gravity to nearest cities, using the population column GC_POP_TOT_2025
gravity_col = "GC_POP_TOT_2025"
gravity_file = os.path.join(results_folder, f"gravity_cities_{gravity_col}.tif")
gravity = ma.calculate_gravity(inR, mcp, inD, gravity_col, outfile=gravity_file)

with rMisc.create_rasterio_inmemory(
    inR.meta, gravity["gravity"][:, :, 0]
) as gravity_raster:
    mapMisc.static_map_raster(
        gravity_raster,
        thresh=[0, 1000, 5000, 10000, 50000, 100000, 500000],
        reverse_colormap=True,
    )