3.2 Extract blocks from road network - orbisgis/h2gis GitHub Wiki

In this use case, we present a robust and efficient method to extract blocks from a road network.

Here, we are using the ROUTE (ROADS in french) layer from BD TOPO, made by IGN.

The study area is the department of Loire-Atlantique, in France.

Input/output tables in numbers:

The input table ROUTE is composed of 235 392 LINESTRINGs.
The resulting table BLOCKS is composed of 54 029 POLYGONs.

SQL instructions


-- Create the grid (cells of 2km). Will be used to clip the roads
DROP TABLE IF EXISTS grid;
CREATE TABLE grid AS SELECT * FROM st_makegrid('ROUTE', 2000, 2000);
CREATE SPATIAL INDEX ON grid(the_geom);

-- Create polygon in cells from cutting roads
DROP TABLE IF EXISTS polygons_cell;
CREATE TABLE polygons_cell(the_geom geometry, index_i int, index_j int, id int) 
  AS (SELECT ST_POLYGONIZE(ST_UNION(ST_COLLECTIONEXTRACT(ST_INTERSECTION(ST_ACCUM(a.the_geom), b.the_geom),2), ST_EXTERIORRING(ST_GEOMETRYN(b.the_geom,1)))) as the_geom, b.ID_COL, b.ID_ROW, b.ID 
FROM ROUTE a, grid b 
  WHERE a.the_geom && b.the_geom GROUP BY b.the_geom);
CREATE INDEX ON polygons_cell(id);
CREATE INDEX ON grid(id);
DROP TABLE IF EXISTS explode_polygons_cell;

-- Explode multi-polygons cell into simple polygon
CREATE TABLE explode_polygons_cell(gid serial, THE_GEOM POLYGON, INDEX_I int,INDEX_J int) 
  AS (SELECT null, the_geom, index_i, index_j FROM ST_EXPLODE('polygons_cell')) 
UNION ALL (SELECT null, a.the_geom,a.ID_COL, a.ID_ROW FROM grid a LEFT JOIN polygons_cell b ON a.ID = b.ID GROUP BY A.ID, b.ID, a.the_geom HAVING b.ID IS NULL);
CREATE INDEX ON explode_polygons_cell(index_i,index_j);
CREATE SPATIAL INDEX ON explode_polygons_cell(the_geom);

-- Find polygon with common segment between adjacent cells
DROP TABLE IF EXISTS polygon_edges;
CREATE TABLE polygon_edges (EDGE_ID SERIAL, START_NODE INT, END_NODE INT) 
  AS SELECT null, A.GID as  START_NODE, B.GID as END_NODE 
FROM explode_polygons_cell A, explode_polygons_cell B 
  WHERE (A.GID<B.GID AND A.THE_GEOM && B.THE_GEOM AND (ABS(A.index_i-B.index_i) = 1 OR ABS(A.index_j-B.index_j) = 1) AND ST_DIMENSION(ST_INTERSECTION(A.THE_GEOM, B.THE_GEOM)) = 1);

DROP TABLE IF EXISTS polygon_edges_EDGE_CC, polygon_edges_NODE_CC;

CALL ST_CONNECTEDCOMPONENTS('polygon_edges', 'undirected');

-- Unify polygons that share a boundary
DROP TABLE IF EXISTS polygons_cell_union;
CREATE TABLE polygons_cell_union 
  AS SELECT ST_UNION(ST_ACCUM(A.THE_GEOM)) AS THE_GEOM 
FROM explode_polygons_cell A, polygon_edges_NODE_CC B 
  WHERE A.GID=B.NODE_ID GROUP BY B.CONNECTED_COMPONENT;

-- Create the blocks
DROP TABLE IF EXISTS blocks;
CREATE TABLE blocks (GID SERIAL, THE_GEOM GEOMETRY, SURFACE DOUBLE) 
  AS (SELECT null, THE_GEOM, ST_AREA(THE_GEOM) FROM polygons_cell_union) 
UNION ALL (SELECT null, a.the_geom, ST_AREA(THE_GEOM) FROM explode_polygons_cell a LEFT JOIN polygon_edges_NODE_CC b ON a.GID = b.NODE_ID GROUP BY A.GID, b.NODE_ID, a.the_geom HAVING NODE_ID IS NULL);

3.2 Extract blocks from road network - orbisgis/h2gis GitHub Wiki

SQL instructions

Illustrations

Input data

Blocks

Zoom on Blocks