Polygonize linestring geometries

0-polygonize.md

Ways to polygonize

Using Qgis

Manual. Uses ST_Polygonize from GDAL. I did not manage to get this to work, but it is manual, mind you.

Using Postgresql, Postgis and ST_Polygonize: elephant.sh

Fast method, but complex query. Uses ST_Polygonize from GDAL. I suggest using this method, works well with insanely large datasets. Remember to use GiST index.

Using Python3 and Shapely+Fiona

Relatively slow, but seems to work, also gives you good control over dangles and cuts.

Using node.js and Turf library: node-makepolygon.js

Works well for small files, but seems to stop on larger (8 MiB +) files. Complains if a LinearRing is less than 4 in length. Handles dangles and cuts, they say, but I have not verified. Buggy on large files.

1-py-makepolygons.py

#!/usr/bin/env python
# Not tested in a long time

import pyproj

from shapely.geometry import Point, LineString, mapping, shape, collection
from shapely.ops import cascaded_union, polygonize, polygonize_full, transform
from shapely import affinity

import pprint

import fiona
from fiona.crs import from_epsg

from argparse import ArgumentParser

yourschema =  {'geometry': 'Polygon','properties': {'test': 'int'}}

def read_shapefile(filename, expected_geometry): # Accepts "LineString" "Point" and more
  print("start read_shapefile()")
  fiona_data = fiona.open(filename, 'r')
  print(fiona_data.schema["geometry"])
  if fiona_data.schema["geometry"] != expected_geometry:
    print("Shapefile for parcel does not contain expected geometry!", expected_geometry)
    exit()
  print("end read_shapefile()")
  return fiona_data

def make_polygons(all_input_lines):
  print("start make_polygons()")
  fiona_input_parcelshapes = [shape(f['geometry']) for f in all_input_lines]
  all_input_parcellines = cascaded_union(fiona_input_parcelshapes)
  result, dangles, cuts, invalids = polygonize_full(all_input_parcellines)
  print("We made polygons, total: ", len(result))
  print("Discarded dangles, total: ", len(dangles))
  print("Discarded cuts, total: ", len(cuts))
  print("Discarded invalids, total: ", len(invalids))
  print("end make_polygons()")
  return result

def main():
  parser = ArgumentParser()
  parser.add_argument("-i", "--input-file", dest="inputfile", help="Input specified file")
  parser.add_argument("-o", "--output-file", dest="outputfile", help="Output specified file")
  parser.add_argument("-v", "--verbose", help="Enable verbose processing", action="store_true")
  args = parser.parse_args()
  inputfile = args.inputfile
  outputfile = args.outputfile
  if args.verbose:
      print("leftalign.py - Verbose processing")
  if args.verbose:
      print ('Input file is ', inputfile)
      print ('Output file is ', outputfile)
  parcels = read_shapefile(inputfile, "LineString")
  map_boundaries = parcels.bounds
  print("Map boundaries found to be (minx, miny, maxx, maxy)=", map_boundaries)

  parcel_polygons = make_polygons(parcels)
  with fiona.open(outputfile, 'w',crs=from_epsg(3857),driver='ESRI Shapefile', schema=yourschema) as output:
    for polygon in parcel_polygons:
      print("Writing polygon to file...", polygon)
      elem = {}
      elem['geometry'] = mapping(polygon) 
      elem['properties'] = {'test': 145}
      output.write(elem)

  print('finish')


if __name__== "__main__":
  main()

2-node-makepolygon.js

var turf = require('@turf/turf');
var fs = require('fs');
var path = require('path');

var inputfile = process.argv[2];
var filename = path.basename(inputfile, '.geojson');
var outputfile = filename+'.polygon';
console.log(filename);

var lines = fs.readFileSync(inputfile);
lines = JSON.parse(lines);
var polygons = turf.polygonize(lines);
fs.writeFileSync(outputfile, JSON.stringify(polygons));
console.log('Saved' + outputfile);

3-postgispolygon.sh

psql -U postgres -d mapDatabase -c "
CREATE TABLE polygons AS
SELECT (ST_Dump(foofoo.polycoll)).geom AS wkb_geometry
FROM (SELECT ST_Polygonize(wkb_geometry) AS polycoll
FROM (SELECT wkb_geometry FROM public.temp) AS foo) As foofoo;
"