Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merge polygons with no features with neighboring polygons #53

Merged
merged 4 commits into from
Sep 21, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 72 additions & 5 deletions fmtm_splitter/fmtm_algorithm.sql
Original file line number Diff line number Diff line change
Expand Up @@ -229,10 +229,10 @@ CREATE TABLE clusteredbuildings AS (
),

-- Cluster the buildings within each splitpolygon. The second term in the
-- call to the ST_ClusterKMeans function is the number of clusters to create,
-- so we're dividing the number of features by a constant (10 in this case)
-- to get the number of clusters required to get close to the right number
-- of features per cluster.
-- call to the ST_ClusterKMeans function is the number of clusters to
-- create, so we're dividing the number of features by a constant
-- (10 in this case) to get the number of clusters required to get close
-- to the right number of features per cluster.
-- TODO: This should certainly not be a hardcoded, the number of features
-- per cluster should come from a project configuration table
buildingstocluster AS (
Expand Down Expand Up @@ -380,13 +380,80 @@ CREATE TABLE taskpolygons AS (
);

-- ALTER TABLE taskpolygons ADD PRIMARY KEY(taskid);

-- Step 1: Identify polygons without any buildings
DROP TABLE IF EXISTS no_building_polygons;
CREATE TABLE no_building_polygons AS (
SELECT
tp.*,
tp.geom AS no_building_geom
FROM taskpolygons AS tp
LEFT JOIN buildings AS b ON ST_INTERSECTS(tp.geom, b.geom)
WHERE b.geom IS NULL
);

-- Step 2: Identify neighboring polygons
DROP TABLE IF EXISTS neighboring_polygons;
CREATE TABLE neighboring_polygons AS (
SELECT
nb.*,
nb.geom AS neighbor_geom,
nb_building_count,
nbp.no_building_geom
FROM taskpolygons AS nb
INNER JOIN no_building_polygons AS nbp
-- Finds polygons that touch each other
ON ST_TOUCHES(nbp.no_building_geom, nb.geom)
LEFT JOIN (
-- Step 3: Count buildings in the neighboring polygons
SELECT
nb.geom,
COUNT(b.geom) AS nb_building_count
FROM taskpolygons AS nb
LEFT JOIN buildings AS b ON ST_INTERSECTS(nb.geom, b.geom)
GROUP BY nb.geom
) AS building_counts
ON nb.geom = building_counts.geom
);

-- Step 4: Find the optimal neighboring polygon to avoid,
-- same polygons with the smallest number of buildings merging into
-- multiple neighboring polygons
DROP TABLE IF EXISTS optimal_neighbors;
CREATE TABLE optimal_neighbors AS (
SELECT
nbp.no_building_geom,
nbp.neighbor_geom
FROM neighboring_polygons AS nbp
WHERE nbp.nb_building_count = (
SELECT MIN(nb_building_count)
FROM neighboring_polygons AS np
WHERE np.no_building_geom = nbp.no_building_geom
)
);

-- Step 5: Merge the small polygons with their optimal neighboring polygons
UPDATE taskpolygons tp
SET geom = ST_UNION(tp.geom, nbp.no_building_geom)
FROM optimal_neighbors AS nbp
WHERE tp.geom = nbp.neighbor_geom;
DELETE FROM taskpolygons
WHERE geom IN (
SELECT no_building_geom
FROM no_building_polygons
);


DROP TABLE IF EXISTS no_building_polygons;
DROP TABLE IF EXISTS neighboring_polygons;

SELECT POPULATE_GEOMETRY_COLUMNS('public.taskpolygons'::regclass);
CREATE INDEX taskpolygons_idx
ON taskpolygons
USING gist (geom);
-- VACUUM ANALYZE taskpolygons;


-- Generate GeoJSON output
SELECT
JSONB_BUILD_OBJECT(
'type', 'FeatureCollection',
Expand Down
2 changes: 1 addition & 1 deletion fmtm_splitter/splitter.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
import geojson
import numpy as np
from geojson import Feature, FeatureCollection, GeoJSON
from osm_rawdata.postgres import PostgresClient
from psycopg2.extensions import connection
from shapely.geometry import Polygon, shape
from shapely.geometry.geo import mapping
Expand All @@ -42,7 +43,6 @@
drop_tables,
insert_geom,
)
from osm_rawdata.postgres import PostgresClient

# Instantiate logger
log = logging.getLogger(__name__)
Expand Down
Loading