create
or replace function st_kdensity(ids bigint [], geoms geometry [])
returns table(id bigint, geom geometry, kdensity integer) as $$ declare mc geometry;
c integer;
k numeric;
begin mc := st_centroid(st_collect(geoms));
c := array_length(ids, 1);
k := sqrt(1 / ln(2));
return query with dist as (
select
t.gid,
t.g,
st_distance(t.g, mc) as distance
from
unnest(ids, geoms) as t(gid, g)
),
md as (
select
percentile_cont(0.5) within group (
order by
distance
) as median
from
dist
),
sd as (
select
sqrt(
sum((st_x(g) - st_x(mc)) ^ 2) / c + sum((st_y(g) - st_y(mc)) ^ 2) / c
) as standard_distance
from
dist
),
sr as (
select
0.9 * least(sd.standard_distance, k * md.median) * c ^(-0.2) as search_radius
from
sd,
md
)
select
gid as id,
g as geom,
kd :: int as kdensity
from
sr,
dist as a,
lateral(
select
count(*) as kd
from
dist _b
where
st_dwithin(a.g, _b.g, sr.search_radius)
) b;
end;
$$ language plpgsql immutable parallel safe;
create table east_village_kde as WITH a AS(
SELECT
array_agg(ogc_fid) as ids,
array_agg(geom) as geoms
FROM
nyc_2015_tree_census
where
st_intersects(
geom,
(
select
geom
from
nyc_neighborhoods
where
neighborhood = 'East Village'
)
)
)
SELECT
b.*
FROM
a,
ST_KDensity(a.ids, a.geoms) b