Finding the Cost of Traversing Through Libya

In this notebook, we will be calculating and visulizing the cost distance of traveling from one population center to another by road while avoiding conflict zones in Libya.

Import and Setup SparkContext



In [ ]:

    
import os
import json
import requests
from functools import partial
import pyproj
import geopyspark as gps

from pyspark import SparkContext
from shapely.geometry import shape, MultiPoint, MultiLineString
from shapely.ops import transform



In [ ]:

    
conf = gps.geopyspark_conf(appName="Libya Weighted Overlay", master="local[*]")
conf.set("spark.hadoop.yarn.timeline-service.enabled", False)
pysc = SparkContext.getOrCreate(conf)

Rasterize Libya Roads to RasterLayer



In [ ]:

    
libya_roads_json = requests.get('https://s3.amazonaws.com/geopyspark-demo/libya/roads.geojson').json()
libya_roads = MultiLineString([shape(geom['geometry']) for geom in libya_roads_json['features']])



In [ ]:

    
# All execution time here is sending WKB over py4j socket
ro = gps.RasterizerOptions(includePartial=True, sampleType='PixelIsArea')

road_raster = gps.rasterize(geoms=list(libya_roads.geoms), 
                            crs="EPSG:3857",
                            zoom=8, 
                            fill_value=1,
                            cell_type=gps.CellType.FLOAT32,
                            options=ro)

road_raster.layer_metadata.bounds

Show Rasterized Roads on a Map



In [ ]:

    
# Pyramid up from base layer
road_pp = road_raster.pyramid(resample_method=gps.ResampleMethod.MAX).cache()



In [ ]:

    
# color map roads 1 to red
roads_cm = gps.ColorMap.from_colors(breaks=[1], color_list=[0xff000080])

# start JVM tile server and serve tiles to map
server = gps.TMS.build(source=road_pp, display=roads_cm)
server.bind("0.0.0.0")
server.url_pattern



In [ ]:

    
from folium import Map, TileLayer

m = Map(tiles='Stamen Toner', location=[27.7351, 17.2283], zoom_start=5)
TileLayer(tiles=server.url_pattern, attr='GeoPySpark Tiles').add_to(m)
m

Cost Distance Based on Road Network



In [ ]:

    
# road network will shape our friction layer
road_friction = road_raster.reclassify(value_map={1:1},
                                       data_type=int,
                                       replace_nodata_with=10)



In [ ]:

    
# starting points for cost distance operation

population_json = requests.get('https://s3.amazonaws.com/geopyspark-demo/libya/population.geojson').json()
population_centers = MultiPoint([shape(geom['geometry']) for geom in population_json['features']])

conflict_json = requests.get('https://s3.amazonaws.com/geopyspark-demo/libya/conflict.geojson').json()
conflict_centers = MultiPoint([shape(feature['geometry']) for feature in conflict_json['features'] if feature['geometry'] != None])

conflict_centers



In [ ]:

    
# Convert population centers data from EPSG:3857 to EPSG:4326 for display on map
project = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:3857'),
    pyproj.Proj(init='epsg:4326'))

population_4326 = transform(project, population_centers)



In [ ]:

    
# Write reprojected data to file

if 'VIRTUAL_ENV' in os.environ:
    !pip3 install geojson
else:
    !pip3 install --user geojson
    
import geojson

with open('/tmp/population-4326.geojson', 'w') as f:
    geojson.dump(geojson.Feature(geometry=population_4326, properties={}), f)
    f.flush()

Cost Distance Between Population Centers



In [ ]:

    
pop_cd = gps.cost_distance(
    friction_layer=road_friction,
    geometries=population_centers, 
    max_distance=1400000.0
)

pop_pp = pop_cd.pyramid()

Cost Distance Between Conflict Centers



In [ ]:

    
con_cd = gps.cost_distance(
    friction_layer=road_friction,
    geometries=conflict_centers, 
    max_distance=1400000.0
)

con_pp = con_cd.pyramid()

Displaying the Weighted Cost Distance Layer With Population Centers



In [ ]:

    
# prepare color map for weighted overlay based on max cost
breaks = [x for x in range(0, 1000000, 10000)]
colors = gps.get_colors_from_matplotlib(ramp_name='viridis', num_colors=len(breaks))
wo_cm = gps.ColorMap.from_colors(breaks=breaks, color_list=colors)



In [ ]:

    
# our weighted layer avoids conflict centers focusing on just population centers
weighted_overlay = (con_pp * 0.0) + (pop_pp * 1.0)

server2 = gps.TMS.build(source=weighted_overlay, display=wo_cm)
server2.bind('0.0.0.0')



In [ ]:

    
from folium import GeoJson

m2 = Map(tiles='Stamen Toner', location=[27.7351, 17.2283], zoom_start=5)
TileLayer(tiles=server2.url_pattern, attr='GeoPySpark Tiles').add_to(m2)
GeoJson("/tmp/population-4326.geojson").add_to(m2)
m2