<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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
Step1:  Rasterize Libya Roads to RasterLayer
Step2:  Show Rasterized Roads on a Map
Step3:  Cost Distance Based on Road Network
Step4:  Cost Distance Between Population Centers
Step5:  Cost Distance Between Conflict Centers
Step6:  Displaying the Weighted Cost Distance Layer With Population Centers
<ASSISTANT_TASK:>
Python Code:
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
conf = gps.geopyspark_conf(appName="Libya Weighted Overlay", master="local[*]")
conf.set("spark.hadoop.yarn.timeline-service.enabled", False)
pysc = SparkContext.getOrCreate(conf)
Explanation: 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
End of explanation
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']])
# 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
Explanation: Rasterize Libya Roads to RasterLayer
End of explanation
# Pyramid up from base layer
road_pp = road_raster.pyramid(resample_method=gps.ResampleMethod.MAX).cache()
# 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
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
Explanation: Show Rasterized Roads on a Map
End of explanation
# road network will shape our friction layer
road_friction = road_raster.reclassify(value_map={1:1},
                                       data_type=int,
                                       replace_nodata_with=10)
# 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
# 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)
# 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()
Explanation: Cost Distance Based on Road Network
End of explanation
pop_cd = gps.cost_distance(
    friction_layer=road_friction,
    geometries=population_centers, 
    max_distance=1400000.0
)
pop_pp = pop_cd.pyramid()
Explanation: Cost Distance Between Population Centers
End of explanation
con_cd = gps.cost_distance(
    friction_layer=road_friction,
    geometries=conflict_centers, 
    max_distance=1400000.0
)
con_pp = con_cd.pyramid()
Explanation: Cost Distance Between Conflict Centers
End of explanation
# 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)
# 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')
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
Explanation: Displaying the Weighted Cost Distance Layer With Population Centers
End of explanation
<END_TASK>