Delete generating new lake centerlines

Delete parts for generating new lake centerlines to avoid confusion.
pull/1235/head
ache051 2021-09-18 15:29:31 +12:00 zatwierdzone przez GitHub
rodzic 7bbca699c5
commit 90495aaf28
Nie znaleziono w bazie danych klucza dla tego podpisu
ID klucza GPG: 4AEE18F83AFDEB23
8 zmienionych plików z 0 dodań i 738 usunięć

Wyświetl plik

@ -1,16 +0,0 @@
FROM python:3.8-slim-buster
RUN apt-get update && apt-get install -y --no-install-recommends \
ca-certificates \
git \
&& rm -rf /var/lib/apt/lists/*
RUN git clone https://github.com/ungarj/label_centerlines.git
COPY *.py ./label_centerlines
WORKDIR ./label_centerlines
RUN pip install -r requirements.txt
RUN python setup.py install
CMD ["./calculate_centerlines.py", "--output_driver", "GeoJSON", "/data/osm_lake_polygon.shp", "/data/osm_lake_centerline.geojson"]
#CMD ["./create_centerlines.py", "--output_driver", "GeoJSON", "/data/osm_lake_polygon.shp", "/data/osm_lake_centerline.geojson"]

Wyświetl plik

@ -1,155 +0,0 @@
#!/usr/bin/env python3
import os
import sys
import argparse
import fiona
import multiprocessing
from shapely.geometry import shape, mapping
from functools import partial
from label_centerlines import get_centerline
def main(args):
parser = argparse.ArgumentParser()
parser.add_argument(
"input_shp",
type=str,
help="input polygons"
)
parser.add_argument(
"output_geojson",
type=str,
help="output centerlines"
)
parser.add_argument(
"--segmentize_maxlen",
type=float,
help="maximum length used when segmentizing polygon borders",
default=0.5
)
parser.add_argument(
"--max_points",
type=int,
help="number of points per geometry allowed before simplifying",
default=3000
)
parser.add_argument(
"--simplification",
type=float,
help="value which increases simplification when necessary",
default=0.05
)
parser.add_argument(
"--smooth",
type=int,
help="smoothness of the output centerlines",
default=5
)
parser.add_argument(
"--output_driver",
type=str,
help="write to 'ESRI Shapefile' or 'GeoJSON' (default)",
default="GeoJSON"
)
parsed = parser.parse_args(args)
input_shp = parsed.input_shp
output_geojson = parsed.output_geojson
segmentize_maxlen = parsed.segmentize_maxlen
max_points = parsed.max_points
simplification = parsed.simplification
smooth_sigma = parsed.smooth
driver = parsed.output_driver
with fiona.open(input_shp, "r") as inp_polygons:
out_schema = inp_polygons.schema.copy()
out_schema['geometry'] = "LineString"
with fiona.open(
output_geojson,
"w",
schema=out_schema,
crs=inp_polygons.crs,
driver=driver
) as out_centerlines:
pool = multiprocessing.Pool()
func = partial(
worker,
segmentize_maxlen,
max_points,
simplification,
smooth_sigma
)
try:
feature_count = 0
for feature_name, output in pool.imap_unordered(
func,
inp_polygons
):
feature_count += 1
if output:
output["properties"]['NAME'] = output["properties"]['NAME']
out_centerlines.write(output)
print( "Written Feature %s: %s" %(
feature_count,
feature_name
))
else:
print("Invalid output for feature", feature_name)
except KeyboardInterrupt:
print("Caught KeyboardInterrupt, terminating workers")
pool.terminate()
except Exception as e:
if feature_name:
print ("%s: FAILED (%s)" %(feature_name, e))
else:
print ("feature: FAILED (%s)" %(e))
raise
finally:
pool.close()
pool.join()
def worker(
segmentize_maxlen,
max_points,
simplification,
smooth_sigma,
feature
):
geom = shape(feature['geometry'])
for name_field in ["name", "Name", "NAME"]:
if name_field in feature["properties"]:
feature_name = feature["properties"][name_field]
feature_id = feature["properties"]["OSM_ID"]
break
else:
feature_name = None
if feature_name:
print ("Processing: ", feature_name)
try:
centerlines_geom = get_centerline(
geom,
segmentize_maxlen=segmentize_maxlen,
max_points=max_points,
simplification=simplification,
smooth_sigma=smooth_sigma
)
except TypeError as e:
print (e)
except:
raise
if centerlines_geom:
return (
feature_name,
{
'properties': feature['properties'],
'geometry': mapping(centerlines_geom)
}
)
else:
return (None, None)
if __name__ == "__main__":
main(sys.argv[1:])

Wyświetl plik

@ -1,30 +0,0 @@
FROM debian:buster-slim
RUN apt-get update && apt-get install -y --no-install-recommends \
gdal-bin \
libgdal-dev \
python \
python-pip \
python-dev \
python-gdal \
git \
&& rm -rf /var/lib/apt/lists/*
RUN apt-get update && apt-get install -y --no-install-recommends \
build-essential \
&& rm -rf /var/lib/apt/lists/*
RUN apt-get update && apt-get install -y --no-install-recommends \
python-scipy \
python-setuptools \
python-shapely \
&& rm -rf /var/lib/apt/lists/*
#WORKDIR /usr/src/app
#RUN git clone https://github.com/ungarj/label_centerlines.git /usr/src/app
RUN pip install wheel
RUN pip install networkx Fiona
COPY *.py /
CMD ["./create_centerlines.py", "--output_driver", "GeoJSON", "/data/osm_lake_polygon.shp", "/data/osm_lake_centerline.geojson"]

Wyświetl plik

@ -1,186 +0,0 @@
#!/usr/bin/env python
# Author: Joachim Ungar <joachim.ungar@eox.at>
#
#-------------------------------------------------------------------------------
# Copyright (C) 2015 EOX IT Services GmbH
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies of this Software or works derived from this Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#-------------------------------------------------------------------------------
import os
import sys
import argparse
import fiona
import multiprocessing
from shapely.geometry import shape, mapping
from functools import partial
from src_create_centerlines import get_centerlines_from_geom
reload(sys)
sys.setdefaultencoding("utf-8")
print sys.version_info
def main(args):
parser = argparse.ArgumentParser()
parser.add_argument(
"input_shp",
type=str,
help="input polygons"
)
parser.add_argument(
"output_geojson",
type=str,
help="output centerlines"
)
parser.add_argument(
"--segmentize_maxlen",
type=float,
help="maximum length used when segmentizing polygon borders",
default=0.5
)
parser.add_argument(
"--max_points",
type=int,
help="number of points per geometry allowed before simplifying",
default=3000
)
parser.add_argument(
"--simplification",
type=float,
help="value which increases simplification when necessary",
default=0.05
)
parser.add_argument(
"--smooth",
type=int,
help="smoothness of the output centerlines",
default=5
)
parser.add_argument(
"--output_driver",
type=str,
help="write to 'ESRI Shapefile' or 'GeoJSON' (default)",
default="GeoJSON"
)
parsed = parser.parse_args(args)
input_shp = parsed.input_shp
output_geojson = parsed.output_geojson
segmentize_maxlen = parsed.segmentize_maxlen
max_points = parsed.max_points
simplification = parsed.simplification
smooth_sigma = parsed.smooth
driver = parsed.output_driver
with fiona.open(input_shp, "r") as inp_polygons:
out_schema = inp_polygons.schema.copy()
out_schema['geometry'] = "LineString"
with fiona.open(
output_geojson,
"w",
schema=out_schema,
crs=inp_polygons.crs,
driver=driver
) as out_centerlines:
pool = multiprocessing.Pool(1)
func = partial(
worker,
segmentize_maxlen,
max_points,
simplification,
smooth_sigma
)
try:
feature_count = 0
for feature_name, output in pool.imap_unordered(
func,
inp_polygons
):
feature_count += 1
feature_name = feature_name.encode("utf-8").strip()
if output:
output["properties"]['NAME'] = output["properties"]['NAME'].encode("utf-8").strip()
out_centerlines.write(output)
print "written feature %s: %s" %(
feature_count,
feature_name
)
else:
print "Invalid output for feature", feature_name
except KeyboardInterrupt:
print "Caught KeyboardInterrupt, terminating workers"
pool.terminate()
except Exception as e:
if feature_name:
print ("%s: FAILED (%s)" %(feature_name, e))
else:
print ("feature: FAILED (%s)" %(e))
raise
finally:
pool.close()
pool.join()
def worker(
segmentize_maxlen,
max_points,
simplification,
smooth_sigma,
feature
):
geom = shape(feature['geometry'])
for name_field in ["NAME"]:
#for name_field in ["name", "Name", "NAME"]:
if name_field in feature["properties"]:
feature_name = feature["properties"][name_field].encode('utf-8')
feature_id = feature["properties"]["OSM_ID"]
break
else:
feature_name = None
if feature_name:
print "processing", feature_name
try:
centerlines_geom = get_centerlines_from_geom(
geom,
segmentize_maxlen=segmentize_maxlen,
max_points=max_points,
simplification=simplification,
smooth_sigma=smooth_sigma
)
except TypeError as e:
print e
except:
raise
if centerlines_geom:
return (
feature_name,
{
'properties': feature['properties'],
'geometry': mapping(centerlines_geom)
}
)
else:
return (None, None)
if __name__ == "__main__":
main(sys.argv[1:])

Wyświetl plik

@ -1,304 +0,0 @@
#!/usr/bin/env python
# Author: Joachim Ungar <joachim.ungar@eox.at>
#
#-------------------------------------------------------------------------------
# Copyright (C) 2015 EOX IT Services GmbH
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies of this Software or works derived from this Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#-------------------------------------------------------------------------------
from shapely.geometry import (
shape,
LineString,
MultiLineString,
Point,
MultiPoint,
mapping
)
from shapely.wkt import loads
import ogr
from scipy.spatial import Voronoi
import networkx as nx
from itertools import combinations
import numpy as np
from scipy.ndimage import filters
from math import *
debug_output = {}
def get_centerlines_from_geom(
geometry,
segmentize_maxlen=0.5,
max_points=3000,
simplification=0.05,
smooth_sigma=5,
debug=False
):
"""
Returns centerline (for Polygon) or centerlines (for MultiPolygons) as
LineString or MultiLineString geometries.
"""
if geometry.geom_type not in ["MultiPolygon", "Polygon"]:
raise TypeError(
"Geometry type must be Polygon or MultiPolygon, not %s" %(
geometry.geom_type
)
)
if geometry.geom_type == "MultiPolygon":
out_centerlines = MultiLineString([
get_centerlines_from_geom(subgeom, segmentize_maxlen)
for subgeom in geometry
if get_centerlines_from_geom(subgeom, segmentize_maxlen) != None
])
return out_centerlines
else:
# Convert Polygon to Linestring.
if len(geometry.interiors) > 0:
boundary = geometry.exterior
else:
boundary = geometry.boundary
#print list(boundary.coords)
if debug:
debug_output['original_points'] = MultiPoint([
point
for point in list(boundary.coords)
])
# Convert to OGR object and segmentize.
ogr_boundary = ogr.CreateGeometryFromWkb(boundary.wkb)
ogr_boundary.Segmentize(segmentize_maxlen)
segmentized = loads(ogr_boundary.ExportToWkt())
# Get points.
points = segmentized.coords
# Simplify segmentized geometry if necessary. This step is required
# as huge geometries slow down the centerline extraction significantly.
tolerance = simplification
while len(points) > max_points:
# If geometry is too large, apply simplification until geometry
# is simplified enough (indicated by the "max_points" value)
tolerance += simplification
simplified = boundary.simplify(tolerance)
points = simplified.coords
if debug:
debug_output['segmentized_points'] = MultiPoint([
point
for point in points
])
# Calculate Voronoi diagram.
vor = Voronoi(points)
if debug:
debug_output['voronoi'] = multilinestring_from_voronoi(
vor,
geometry
)
# The next three steps are the most processing intensive and probably
# not the most efficient method to get the skeleton centerline. If you
# have any recommendations, I would be very happy to know.
# Convert to networkx graph.
graph = graph_from_voronoi(vor, geometry)
graph.edges(18,21)
# Get end nodes from graph.
end_nodes = get_end_nodes(graph)
if len(end_nodes) < 2:
return None
# Get longest path.
longest_paths = get_longest_paths(
end_nodes,
graph
)
# get least curved path.
best_path = get_least_curved_path(longest_paths[:5], vor.vertices)
#print (best_path == longest_paths[0])
#best_path = longest_paths[0]
centerline = LineString(vor.vertices[best_path])
if debug:
debug_output['centerline'] = centerline
# Simplify again to reduce number of points.
# simplified = centerline.simplify(tolerance)
# centerline = simplified
# Smooth out geometry.
centerline_smoothed = smooth_linestring(centerline, smooth_sigma)
out_centerline = centerline_smoothed
#print out_centerline
return out_centerline
def smooth_linestring(linestring, smooth_sigma):
"""
Uses a gauss filter to smooth out the LineString coordinates.
"""
smooth_x = np.array(filters.gaussian_filter1d(
linestring.xy[0],
smooth_sigma)
)
smooth_y = np.array(filters.gaussian_filter1d(
linestring.xy[1],
smooth_sigma)
)
smoothed_coords = np.hstack((smooth_x, smooth_y))
smoothed_coords = zip(smooth_x, smooth_y)
linestring_smoothed = LineString(smoothed_coords)
return linestring_smoothed
def get_longest_paths(nodes, graph):
"""
Returns longest path of all possible paths between a list of nodes.
"""
paths = []
distances = []
possible_paths = list(combinations(nodes, r=2))
for node1, node2 in possible_paths:
try:
path = nx.shortest_path(graph, node1, node2, "weight")
except Exception,e:
path = []
if len(path)>1:
distance = get_path_distance(path, graph)
paths.append(path)
distances.append(distance)
paths_sorted = [x for (y,x) in sorted(zip(distances, paths), reverse=True)]
# longest_path = paths_sorted[0]
# return longest_path
return paths_sorted
def get_least_curved_path(paths, vertices):
angle_sums = []
for path in paths:
path_angles = get_path_angles(path, vertices)
angle_sum = abs(sum(path_angles))
angle_sums.append(angle_sum)
paths_sorted = [x for (y,x) in sorted(zip(angle_sums, paths))]
return paths_sorted[0]
def get_path_angles(path, vertices):
angles = []
prior_line = None
next_line = None
for index, point in enumerate(path):
if index > 0 and index < len(path)-1:
prior_point = vertices[path[index-1]]
current_point = vertices[point]
next_point = vertices[path[index+1]]
angles.append(
get_angle(
(prior_point, current_point), (current_point, next_point)
)
)
return angles
def get_angle(line1, line2):
v1 = line1[0] - line1[1]
v2 = line2[0] - line2[1]
angle = np.math.atan2(np.linalg.det([v1,v2]),np.dot(v1,v2))
return np.degrees(angle)
def get_path_distance(path, graph):
"""
Returns weighted path distance.
"""
distance = 0
for i,w in enumerate(path):
j=i+1
if j<len(path):
#print graph.get_edge_data(path[i],path[j])["weight"]
#distance += round(graph.edges(path[i],path[j])["weight"], 6)
distance += round(graph.get_edge_data(path[i],path[j])["weight"], 6)
return distance
def get_end_nodes(graph):
"""
Returns list of nodes with just one neighbor node.
"""
# nodelist = []
# for i in list(graph.nodes()):
# if len(list(graph.neighbors(i)))==1:
# nodelist.append(i)
nodelist = [
i
for i in list(graph.nodes())
if len(list(graph.neighbors(i)))==1
]
return nodelist
def graph_from_voronoi(vor, geometry):
"""
Creates a networkx graph out of all Voronoi ridge vertices which are inside
the original geometry.
"""
graph = nx.Graph()
for i in vor.ridge_vertices:
if i[0]>-1 and i[1]>-1:
point1 = Point(vor.vertices[i][0])
point2 = Point(vor.vertices[i][1])
# Eliminate all points outside our geometry.
if point1.within(geometry) and point2.within(geometry):
dist = point1.distance(point2)
graph.add_nodes_from([i[0], i[1]])
graph.add_edge(i[0], i[1], weight=dist)
return graph
def multilinestring_from_voronoi(vor, geometry):
"""
Creates a MultiLineString geometry out of all Voronoi ridge vertices which
are inside the original geometry.
"""
linestrings = []
for i in vor.ridge_vertices:
if i[0]>-1 and i[1]>-1:
point1 = Point(vor.vertices[i][0])
point2 = Point(vor.vertices[i][1])
# Eliminate all points outside our geometry.
if point1.within(geometry) and point2.within(geometry):
linestring = LineString([point1, point2])
linestrings.append(linestring)
multilinestring = MultiLineString(linestrings)
return multilinestring
if __name__ == "__main__":
main(sys.argv[1:])

Wyświetl plik

@ -1,4 +0,0 @@
FROM postgis/postgis:9.6-2.5-alpine
RUN apk add libintl
COPY export-shapefile.sh /
CMD ["./export-shapefile.sh"]

Wyświetl plik

@ -1,16 +0,0 @@
#!/bin/bash
set -o errexit
set -o pipefail
set -o nounset
function export_shp() {
local lake_shapefile="/data/osm_lake_polygon.shp"
local query="SELECT osm_id, name, name_en, ST_SimplifyPreserveTopology(geometry, 100) AS geometry FROM osm_water_polygon WHERE area > 2 * 1000 * 1000 AND ST_GeometryType(geometry)='ST_Polygon' AND name <> '' ORDER BY area DESC"
pgsql2shp -f "$lake_shapefile" \
-h "$PGHOST" \
-u "$PGUSER" \
-P "$PGPASSWORD" \
"$PGDATABASE" "$query"
}
export_shp

Wyświetl plik

@ -1,27 +0,0 @@
#!/usr/bin/env bash
set -o errexit
set -o pipefail
set -o nounset
echo "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
echo " Loading OMT postgis extensions"
echo "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
for db in template_postgis "$POSTGRES_DB"; do
echo "Loading extensions into $db"
PGUSER="$POSTGRES_USER" psql --dbname="$db" <<-'EOSQL'
-- Cleanup. Ideally parent container shouldn't pre-install those.
DROP EXTENSION IF EXISTS postgis_tiger_geocoder;
DROP EXTENSION IF EXISTS postgis_topology;
-- These extensions are already loaded by the parent docker
CREATE EXTENSION IF NOT EXISTS postgis;
CREATE EXTENSION IF NOT EXISTS fuzzystrmatch;
-- Extensions needed for OpenMapTiles
CREATE EXTENSION IF NOT EXISTS hstore;
CREATE EXTENSION IF NOT EXISTS unaccent;
CREATE EXTENSION IF NOT EXISTS osml10n;
CREATE EXTENSION IF NOT EXISTS gzip;
EOSQL
done