# -*- coding: utf-8 -*-
"""
The code in this module is partly based on third party code which has been
licensed under GNU-GPL3. The following functions are copied and adapted from:
https://github.com/FRESNA/vresutils,
Copyright 2015-2017 Frankfurt Institute for Advanced Studies
* _shape2poly()
* simplify_poly()
* nuts()
"""
import os
from collections import OrderedDict
from functools import partial
from itertools import product
from itertools import takewhile
from operator import attrgetter
from operator import itemgetter
import pandas as pd
try:
from shapely.geometry import LinearRing
from shapely.geometry import MultiPolygon
from shapely.geometry import Polygon
from shapely.geometry import shape
from shapely.ops import transform
from shapely.prepared import prep
except ImportError:
raise ImportError("Need to install shapely to use geometry module!")
try:
import pyproj
except ImportError:
raise ImportError("Need to install pyproj to use geometry module!")
try:
from geojson import Feature
from geojson import FeatureCollection
from geojson import dump
from geojson import load
except ImportError:
raise ImportError("Need to install geojson to use geometry module!")
try:
import scipy.sparse as sparse
except ImportError:
raise ImportError("Need to install scipy to use geometry module!")
try:
import shapefile
except ImportError:
raise ImportError("Need to install pyshp to use geometry module!")
import numpy as np
[docs]
def read_geometries(filename, directory="data/geometries"):
"""
Reads geometry resources from the datapackage. Data may either be stored
in geojson format or as WKT representation in CSV-files.
Parameters
----------
filename: string
Name of the elements to be read, for example `buses.geojson`
directory: string
Directory where the file is located. Default: `data/geometries`
Returns
-------
pd.Series
"""
path = os.path.join(directory, filename)
if os.path.splitext(filename)[1] == ".geojson":
if os.path.exists(path):
with open(path, "r") as infile:
features = load(infile)["features"]
names = [f["properties"]["name"] for f in features]
geometries = [shape(f["geometry"]) for f in features]
geometries = pd.Series(dict(zip(names, geometries)))
if os.path.splitext(filename)[1] == ".csv":
if os.path.exists(path):
geometries = pd.read_csv(path, sep=";", index_col=["name"])
else:
geometries = pd.Series(name="geometry")
geometries.index.name = "name"
return geometries
[docs]
def write_geometries(filename, geometries, directory="data/geometries"):
"""Writes geometries to filesystem.
Parameters
----------
filename: string
Name of the geometries stored, for example `buses.geojson`
geometries: pd.Series
Index entries become name fields in GeoJSON properties.
directory: string
Directory where the file is stored. Default: `data/geometries`
Returns
-------
path: string
Returns the path where the file has been stored.
"""
path = os.path.join(directory, filename)
if os.path.splitext(filename)[1] == ".geojson":
features = FeatureCollection(
[
Feature(geometry=v, properties={"name": k})
for k, v in geometries.iteritems()
]
)
if os.path.exists(path):
with open(path) as infile:
existing_features = load(infile)["features"]
names = [f["properties"]["name"] for f in existing_features]
assert all(i not in names for i in geometries.index), (
"Cannot " "create duplicate entries in %s." % filename
)
features["features"] += existing_features
with open(path, "w") as outfile:
dump(features, outfile)
if os.path.splitext(filename)[1] == ".csv":
if os.path.exists(path):
existing_geometries = read_geometries(filename, directory)
geometries = pd.concat(
[existing_geometries, geometries], verify_integrity=True
)
geometries.index.name = "name"
geometries.to_csv(path, sep=";", header=True)
return path
def _shape2poly(sh, tolerance=0.03, minarea=0.03, projection=None):
"""
Notes
-----
Copied from: https://github.com/FRESNA/vresutils,
Copyright 2015-2017 Frankfurt Institute for Advanced Studies
"""
if len(sh.points) == 0:
return None
if projection is None:
pts = sh.points
elif projection == "invwgs":
pts = np.asarray(
_shape2poly.wgs(*np.asarray(sh.points).T, inverse=True)
).T
else:
raise TypeError("Unknown projection {}".format(projection))
minlength = 2 * np.pi * np.sqrt(minarea / np.pi)
def parts2polys(parts):
rings = list(map(LinearRing, parts))
while rings:
exterior = rings.pop(0)
interiors = list(takewhile(attrgetter("is_ccw"), rings))
rings = rings[len(interiors) :]
yield Polygon(
exterior, [x for x in interiors if x.length > minlength]
)
polys = sorted(
parts2polys(np.split(pts, sh.parts[1:])),
key=attrgetter("area"),
reverse=True,
)
mainpoly = polys[0]
mainlength = np.sqrt(mainpoly.area / (2.0 * np.pi))
if mainpoly.area > minarea:
mpoly = MultiPolygon(
[
p
for p in takewhile(lambda p: p.area > minarea, polys)
if mainpoly.distance(p) < mainlength
]
)
else:
mpoly = mainpoly
return simplify_poly(mpoly, tolerance)
_shape2poly.wgs = pyproj.Proj(
"+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84" " +units=m +no_defs"
)
[docs]
def simplify_poly(poly, tolerance):
"""
Notes
-----
Copied from: https://github.com/FRESNA/vresutils,
Copyright 2015-2017 Frankfurt Institute for Advanced Studies
"""
if tolerance is None:
return poly
else:
return poly.simplify(tolerance, preserve_topology=True)
[docs]
def nuts(filepath=None, nuts=0, subset=None, tolerance=0.03, minarea=1.0):
"""
Reads shapefile with nuts regions and converts to polygons
Returns
-------
OrderedDict
Country keys as keys of dict and shapely polygons as corresponding
values
Notes
-----
Copied from: https://github.com/FRESNA/vresutils,
Copyright 2015-2017 Frankfurt Institute for Advanced Studies
"""
sf = shapefile.Reader(filepath)
nuts = OrderedDict(
sorted(
[
(rec[0], _shape2poly(sh, tolerance, minarea))
for rec, sh in zip(sf.iterRecords(), sf.iterShapes())
if rec[1] == nuts
],
key=itemgetter(0),
)
)
return nuts
[docs]
def reproject(
geom, fr=pyproj.Proj(init="EPSG:4326"), to=pyproj.Proj(init="EPSG:25832")
):
"""
Notes
-----
Copied and adapted from: https://github.com/FRESNA/vresutils,
Copyright 2015-2017 Frankfurt Institute for Advanced Studies
"""
reproject_pts = partial(pyproj.transform, fr, to)
return transform(reproject_pts, geom)
[docs]
def Shapes2Shapes(
orig, dest, normed=True, equalarea=False, prep_first=True, **kwargs
):
"""
Notes
-----
Copied from: https://github.com/FRESNA/vresutils,
Copyright 2015-2017 Frankfurt Institute for Advanced Studies
"""
if equalarea:
dest = list(map(reproject, dest))
orig = list(map(reproject, orig))
if prep_first:
orig_prepped = list(map(prep, orig))
else:
orig_prepped = orig
transfer = sparse.lil_matrix((len(dest), len(orig)), dtype=np.float)
for i, j in product(range(len(dest)), range(len(orig))):
if orig_prepped[j].intersects(dest[i]):
area = orig[j].intersection(dest[i]).area
transfer[i, j] = area / dest[i].area
# sum of input vectors must be preserved
if normed:
ssum = np.squeeze(np.asarray(transfer.sum(axis=0)), axis=0)
for i, j in zip(*transfer.nonzero()):
transfer[i, j] /= ssum[j]
return transfer
[docs]
def intersects(geom, labels, geometries):
""" """
for label, geom_to_check in zip(labels, geometries):
if geom.intersects(geom_to_check):
return label
return float("NaN")