create table pharmacy_stats as
select
p.id,
p.amenity,
p.brand,
p.name,
p.geom,
-- Add a buffer of 800 meters for later use
st_transform(st_buffer(st_transform(p.geom, 3857), 800), 4326) as buffer,
-- Get the total proportional population for each buffer
sum(
bgs.population * (
st_area(
st_intersection(
st_transform(st_buffer(st_transform(p.geom, 3857), 800), 4326),
bgs.geom
)
) / st_area(bgs.geom)
)
) as pop
from
nyc_pharmacies p
join nyc_2021_census_block_groups bgs on st_intersects(p.geom, st_transform(bgs.geom, 4326))
-- Get just the stores for CVS and Duane Reade
where
p.name ilike '%duane%'
or p.name ilike '%cvs%'
group by
p.id,
p.amenity,
p.brand,
p.name,
p.geom
select
*
from
nyc_pharmacies
where
name ilike '%duane%'
or name ilike '%cvs%'
-- CTE with all Duane Reade locations
with dr as (
select
id,
amenity,
brand,
name,
geom,
st_transform(geom, 3857) as geom_3857
from
nyc_pharmacies
where
name ilike '%duane%'
),
-- CTE with all CVS locations
cvs as (
select
id,
amenity,
brand,
name,
geom,
st_transform(geom, 3857) as geom_3857
from
nyc_pharmacies
where
name ilike '%cvs%'
),
-- Find all Duane Reade locations within 200 meters of a CVS
remove_nearest as (
select
dr.*
from
dr,
cvs
where
st_dwithin(dr.geom_3857, cvs.geom_3857, 200)
)
select
count(*)
from
remove_nearest
with dr as (
select
*
from
pharmacy_stats
where
name ilike '%duane%'
)
select
dr.*,
-- Find the area of overlap between the two buffer groups
st_area(
st_intersection(pharmacy_stats.buffer, dr.buffer)
) / st_area(dr.buffer)
from
dr
join pharmacy_stats on st_intersects(dr.buffer, pharmacy_stats.buffer)
where
-- Find the number of pharmacies that have an overlap greater than 75%
st_area(
st_intersection(pharmacy_stats.buffer, dr.buffer)
) / st_area(dr.buffer) >.75
and pharmacy_stats.name ilike '%cvs%'
with a as (
select
-- Union all buffers to find the "before" scenario total area
st_union(
st_transform(st_buffer(st_transform(p.geom, 3857), 800), 4326)
) as buffer,
st_area(
st_union(st_buffer(st_transform(p.geom, 3857), 800))
) as area
from
nyc_pharmacies p
join nyc_2021_census_block_groups bgs on st_intersects(
st_transform(st_buffer(st_transform(p.geom, 3857), 800), 4326),
st_transform(bgs.geom, 4326)
)
where
p.name ilike '%duane%'
or p.name ilike '%cvs%'
)
select
a.*,
-- Find the total populatino of all combined buffers
sum(
bgs.population * (
st_area(st_intersection(a.buffer, bgs.geom)) / st_area(bgs.geom)
)
) as pop
from
a
join nyc_2021_census_block_groups bgs on st_intersects(a.buffer, st_transform(bgs.geom, 4326))
group by
1,
2
create table duane_reade_ids as with overlap as (
with dr as (
select
*
from
pharmacy_stats
where
name ilike '%duane%'
)
-- Find all the Duane Reade locations that have an overlap of over 75%
select
dr.id
from
dr
join pharmacy_stats on st_intersects(dr.buffer, pharmacy_stats.buffer)
where
st_area(
st_intersection(pharmacy_stats.buffer, dr.buffer)
) / st_area(dr.buffer) >.75
and pharmacy_stats.name ilike '%cvs%'
),
-- Select all Duane Reade stores
dr as (
select
id,
amenity,
brand,
name,
geom,
st_transform(geom, 3857) as geom_3857
from
nyc_pharmacies
where
name ilike '%duane%'
),
-- Select all CVS stores
cvs as (
select
id,
amenity,
brand,
name,
geom,
st_transform(geom, 3857) as geom_3857
from
nyc_pharmacies
where
name ilike '%cvs%'
),
-- Remove all the Duane Reade locations within 200 meters of a CVS
remove_nearest as (
select
dr.id
from
dr,
cvs
where
st_dwithin(dr.geom_3857, cvs.geom_3857, 200)
)
select
id
from
remove_nearest
union
select
id
from
overlap
-- We run the same query but we remove IDs not in the table we just created
with p as (
select
*
from
nyc_pharmacies
where
id not in (
select
id
from
duane_reade_ids
)
),
a as (
select
st_union(
st_transform(st_buffer(st_transform(p.geom, 3857), 800), 4326)
) as buffer,
st_area(
st_union(st_buffer(st_transform(p.geom, 3857), 800))
) as area
from
nyc_pharmacies p
join nyc_2021_census_block_groups bgs on st_intersects(
st_transform(st_buffer(st_transform(p.geom, 3857), 800), 4326),
st_transform(bgs.geom, 4326)
)
where
p.name ilike '%duane%'
or p.name ilike '%cvs%'
)
select
a.*,
sum(
bgs.population * (
st_area(st_intersection(a.buffer, bgs.geom)) / st_area(bgs.geom)
)
) as pop
from
a
join nyc_2021_census_block_groups bgs on st_intersects(a.buffer, st_transform(bgs.geom, 4326))
group by
1,
2