ogr2ogr \
-f PostgreSQL PG:"host=localhost user=docker password=docker dbname=routing port=25432" \
FDNY_Firehouse_Listing.csv \
-dialect sqlite -sql \
"select *, makepoint(cast(Longitude as float), cast(Latitude as float), 4326) as geom from FDNY_Firehouse_Listing" \
-nln nyc_fire_stations -lco GEOMETRY_NAME=geom
ogr2ogr \
-f PostgreSQL PG:"host=localhost user=docker password=docker \
dbname=routing port=25432" nyc_mappluto_22v3_shp/MapPLUTO.shp \
-nln building_footprints -lco GEOMETRY_NAME=geom \
-nlt MULTIPOLYGON -mapFieldType Real=String
ogr2ogr \
-f PostgreSQL PG:"host=localhost user=docker password=docker \
dbname=routing port=25432" \
-dialect SQLite -sql "SELECT neighborhood, boroughcode, borough, \
ST_Area(st_collect(st_transform(geometry, 3857))) AS area, \
st_collect(geometry) as geom FROM nyc_hoods \
group by neighborhood, boroughcode, borough" \
nyc_hoods.geojson \
-nln nyc_neighborhoods -lco GEOMETRY_NAME=geom -nlt PROMOTE_TO_MULTI
create table new_stations (name text, geom geometry)
insert into
new_stations (name, geom)
values
(
'Withers and Woodpoint',
st_setsrid(st_makepoint(-73.941455, 40.717628), 4326)
),
(
'Grand and Graham',
st_setsrid(st_makepoint(-73.943791, 40.711584), 4326)
),
(
'Berry and N 11th',
st_setsrid(st_makepoint(-73.9566498, 40.7207395), 4326)
)
create table all_stations as
select
facilityname as name,
geom
from
nyc_fire_stations
where
-- Finding all the stations in our three target neighborhoods
st_intersects(
geom,
(
select
st_union(st_transform(geom, 4326))
from
nyc_neighborhoods
where
neighborhood in (
'Williamsburg',
'East Williamsburg',
'Greenpoint'
)
)
)
union
select
*
from
new_stations
create table bk_blgds as
select
*
from
building_footprints
where
st_intersects(
st_transform(geom, 4326),
-- Find all the buildings that are in the three neighborhoods
-- and within 100 meters
(
select
st_buffer(st_union(geom) :: geography, 100) :: geometry
from
nyc_neighborhoods
where
neighborhood in (
'Williamsburg',
'East Williamsburg',
'Greenpoint'
)
)
)
create table fire_odm as
with starts as (
-- Aggregate the source way ID into an array for
-- the closest way ID to each station
select
array_agg(z.source)
from
all_stations
cross join lateral (
select
source
from
ways
join car_config c using (tag_id)
where
-- Here we are picking all the way tags that
-- are not in this list
c.tag_value not in (
'track',
'bridleway',
'bus_guideway',
'byway',
'cycleway',
'path',
'track',
'grade1',
'grade2',
'grade3',
'grade4',
'grade5',
'unclassified',
'footway',
'pedestrian',
'steps'
)
order by
the_geom <-> all_stations.geom
limit
1
) z
),
-- Here we do the same with the destinations for the buildings
destinations as (
select
array_agg(z.source)
from
bk_blgds
cross join lateral (
select
source
from
ways
join car_config c using (tag_id)
where
c.tag_value not in (
'track',
'bridleway',
'bus_guideway',
'byway',
'cycleway',
'path',
'track',
'grade1',
'grade2',
'grade3',
'grade4',
'grade5',
'unclassified',
'footway',
'pedestrian',
'steps'
)
order by
ways.the_geom <-> st_transform(st_centroid(bk_blgds.geom), 4326)
limit
1
) z
-- We only select the buildings in the neighborhood boundaries
where
st_intersects(
st_transform(geom, 4326),
(
select
st_buffer(st_union(geom) :: geography, 100) :: geometry
from
nyc_neighborhoods
where
neighborhood in (
'williamsburg',
'east williamsburg',
'greenpoint'
)
)
)
)
select
*
from
-- Here we pass these arguments to the pgr_bdDijkstraCost function
pgr_bddijkstracost(
$$
select
id,
source,
target,
cost_s as cost,
reverse_cost_s as reverse_cost
from
ways
join car_config using (tag_id)
where
st_intersects(
st_transform(the_geom, 4326),
(
select
st_buffer(st_union(geom) :: geography, 100) :: geometry
from
nyc_neighborhoods
where
neighborhood in (
'williamsburg',
'east williamsburg',
'greenpoint'
)
)
) $$,
-- We pass the arrays from above as the two arguments
(
select
*
from
starts
),
(
select
*
from
destinations
),
true
);
create
or replace function spopt_pmedian(
odm_table TEXT,
optimal_facilities INT,
clients_table TEXT,
facilities_table TEXT
)
returns table (
facility TEXT,
end_vid INT,
ogc_fid INT,
geom GEOMETRY
) AS $$
import spopt
from spopt.locate import PMedian
import numpy
import pulp
import pandas as pd
odm = plpy.execute(f''' select * from {odm_table} ''')
odm_new = []
for i in odm:
odm_new.append({'start_vid': i['start_vid'], 'end_vid': i['end_vid'], 'agg_cost': i['agg_cost']})
solver = pulp.PULP_CBC_CMD(keepFiles=True)
df = pd.DataFrame.from_dict(odm_new)
data = df.pivot_table(index='end_vid', columns='start_vid', values='agg_cost', fill_value=1000000).values
vals = df.pivot_table(index='end_vid', columns='start_vid', values='agg_cost', fill_value=1000000)
clients = data.shape[0]
weights = [1] * clients
pmedian_from_cm = PMedian.from_cost_matrix(
numpy.array(data),
numpy.array(weights),
p_facilities=int(optimal_facilities),
name="p-median-network-distance" )
pmedian_from_cm = pmedian_from_cm.solve(solver)
station_ids = plpy.execute(f'''
select
z.source as end_vid,
{facilities_table}.name
from
{facilities_table}
cross join lateral (
SELECT
source
FROM
ways
join configuration c using (tag_id)
where
c.tag_value not in (
'track',
'bridleway',
'bus_guideway',
'byway',
'cycleway',
'path',
'track',
'grade1',
'grade2',
'grade3',
'grade4',
'grade5',
'unclassified',
'footway',
'pedestrian',
'steps'
)
ORDER BY
ways.the_geom <-> st_transform(st_centroid({facilities_table}.geom), 4326)
limit
1
) z
''')
stations = []
for i in station_ids:
stations.append(i)
stations_df = pd.DataFrame.from_dict(stations)
cleaned_points = []
for i in pmedian_from_cm.fac2cli:
if len(i) > 0:
z = stations_df[stations_df['end_vid'] == vals.columns[pmedian_from_cm.fac2cli.index(i)]]['name'].values[0]
for j in i:
struct = {'facility': z, 'end_vid': list(vals.index)[j]}
cleaned_points.append(struct)
df_startids = pd.DataFrame.from_dict(cleaned_points)
orig_data = plpy.execute(f'''
select
z.source as end_vid,
{clients_table}.ogc_fid,
st_transform(st_centroid({clients_table}.geom), 4326) as geom
from
{clients_table}
cross join lateral (
SELECT
source
FROM
ways
join car_config c using (tag_id)
where
c.tag_value not in (
'track',
'bridleway',
'bus_guideway',
'byway',
'cycleway',
'path',
'track',
'grade1',
'grade2',
'grade3',
'grade4',
'grade5',
'unclassified',
'footway',
'pedestrian',
'steps'
)
ORDER BY
ways.the_geom <-> st_transform(st_centroid({clients_table}.geom), 4326)
limit
1
) z
''')
orig_formatted = []
for i in orig_data:
orig_formatted.append(i)
orig_df = pd.DataFrame.from_dict(orig_formatted)
final_df = orig_df.merge(df_startids, how='left', on='end_vid')
final_df = final_df.replace(numpy.nan, None)
return final_df.to_dict(orient='records')
$$
LANGUAGE 'plpython3u';
create table bk_final as
select
*
from
spopt_pmedian('fire_odm', 5, 'bk_blgds', 'all_stations');
create table final_stations as
select
*
from
all_stations
where
name in (
select
distinct facility
from
bk_final
)