Using ArcGIS
Import the ArcGIS API for Python and other packages
[1]:
import warnings
warnings.filterwarnings('ignore')
from arcgis.features import GeoAccessor
from arcgis.gis import GIS
import pulp
from allagash import Coverage, Problem
Load the data into spatially enabled dataframes
[2]:
d = GeoAccessor.from_featureclass("sample_data/demand_point.shp")
s = GeoAccessor.from_featureclass("sample_data/facility_service_areas.shp")
s2 = GeoAccessor.from_featureclass("sample_data/facility2_service_areas.shp")
d.spatial.sr = s.spatial.sr = s2.spatial.sr = "32612"
Plot the data
[3]:
# NBVAL_IGNORE_OUTPUT
gis = GIS()
m = gis.map("Millcreek, UT", zoomlevel=10)
m.add_layer(d.spatial.to_featureset())
m.add_layer(s.spatial.to_featureset())
m.add_layer(s2.spatial.to_featureset())
m
Create the coverage
[4]:
coverage1 = Coverage.from_spatially_enabled_dataframes(d, s, "GEOID10", "ORIG_ID")
coverage2 = Coverage.from_spatially_enabled_dataframes(d, s2, "GEOID10", "ORIG_ID")
Create the problem
[5]:
problem = Problem.lscp([coverage1, coverage2])
Solve the problem
[6]:
problem.solve(pulp.GLPK(msg=False));
Plot the selected locations
[7]:
# NBVAL_IGNORE_OUTPUT
selected_locations = s.query(f"ORIG_ID in ({','.join(problem.selected_supply(coverage1))})")
selected_locations2 = s2.query(f"ORIG_ID in ({','.join(problem.selected_supply(coverage2))})")
m = gis.map("Millcreek, UT", zoomlevel=10)
m.add_layer(d.spatial.to_featureset())
m.add_layer(selected_locations.spatial.to_featureset())
m.add_layer(selected_locations2.spatial.to_featureset())
m
Find the total number of locations required
[8]:
print(problem.pulp_problem.objective.value())
24