diff --git a/data/chessboard/building_to_edge.py b/data/chessboard/building_to_edge.py
index 191bc7b..d4857d3 100644
--- a/data/chessboard/building_to_edge.py
+++ b/data/chessboard/building_to_edge.py
@@ -1,7 +1,5 @@
import geopandas
-import pandas as pd
import pandashp as pdshp
-import shutil
import os
diff --git a/rivus/utils/pandashp.py b/rivus/utils/pandashp.py
index 3294c57..34e0ecb 100644
--- a/rivus/utils/pandashp.py
+++ b/rivus/utils/pandashp.py
@@ -1,143 +1,12 @@
-""" pandashp: read/write shapefiles to/from special DataFrames
-
-Offers two functions read_shp and write_shp that convert ESRI shapefiles to
-pandas DataFrames that can be manipulated at will and then written back to
-shapefiles. Opens up data manipulation capabilities beyond a simple GIS field
-calculator.
-
-Usage:
- import pandashp as pdshp
- # calculate population density from shapefile of cities (stupid, I know)
- cities = pdshp.read_shp('cities_germany_projected')
- cities['popdens'] = cities['population'] / cities['area']
- pdshp.write_shp(cities, 'cities_germany_projected_popdens')
-
-"""
-
-try:
- from itertools import izip as zip
-except ImportError: # zip is a builtin in Python 3.x
- pass
-import numpy as np
import pandas as pd
-import shapefile
from . import shapelytools
import warnings
-from shapely.geometry import LineString, Point, Polygon
-
-
-def read_shp(filename):
- """Read shapefile to dataframe w/ geometry.
-
- Args:
- filename: ESRI shapefile name to be read (without .shp extension)
-
- Returns:
- pandas DataFrame with column geometry, containing individual shapely
- Geometry objects (i.e. Point, LineString, Polygon) depending on
- the shapefiles original shape type
-
- """
- sr = shapefile.Reader(filename)
-
- cols = sr.fields[:] # [:] = duplicate field list
- if cols[0][0] == 'DeletionFlag':
- cols.pop(0)
- cols = [col[0] for col in cols] # extract field name only
- cols.append('geometry')
-
- records = [row for row in sr.iterRecords()]
-
- if sr.shapeType == shapefile.POLYGON:
- geometries = [Polygon(shape.points)
- if len(shape.points) > 2 else np.NaN # invalid geometry
- for shape in sr.iterShapes()]
- elif sr.shapeType == shapefile.POLYLINE:
- geometries = [LineString(shape.points) for shape in sr.iterShapes()]
- elif sr.shapeType == shapefile.POINT:
- geometries = [Point(*shape.points[0]) for shape in sr.iterShapes()]
- else:
- raise NotImplementedError
-
- data = [r + [g] for r, g in zip(records, geometries)]
-
- df = pd.DataFrame(data, columns=cols)
- df = df.convert_objects(convert_numeric=True)
-
- if np.NaN in geometries:
- # drop invalid geometries
- df = df.dropna(subset=['geometry'])
- num_skipped = len(geometries) - len(df)
- warnings.warn('Skipped {} invalid geometrie(s).'.format(num_skipped))
- return df
-
-
-def write_shp(filename, dataframe, write_index=True):
- """Write dataframe w/ geometry to shapefile.
-
- Args:
- filename: ESRI shapefile name to be written (without .shp extension)
- dataframe: a pandas DataFrame with column geometry and homogenous
- shape types (Point, LineString, or Polygon)
- write_index: add index as column to attribute tabel (default: true)
-
- Returns:
- Nothing.
-
- """
-
- df = dataframe.copy()
- if write_index:
- df.reset_index(inplace=True)
-
- # split geometry column from dataframe
- geometry = df.pop('geometry')
-
- # write geometries to shp/shx, according to geometry type
- if isinstance(geometry.iloc[0], Point):
- sw = shapefile.Writer(shapefile.POINT)
- for point in geometry:
- sw.point(point.x, point.y)
-
- elif isinstance(geometry.iloc[0], LineString):
- sw = shapefile.Writer(shapefile.POLYLINE)
- for line in geometry:
- sw.line([list(line.coords)])
-
- elif isinstance(geometry.iloc[0], Polygon):
- sw = shapefile.Writer(shapefile.POLYGON)
- for polygon in geometry:
- sw.poly([list(polygon.exterior.coords)])
- else:
- raise NotImplementedError
-
- # add fields for dbf
- for k, column in enumerate(df.columns):
- # unicode strings freak out pyshp, so remove u'..'
- column = str(column)
-
- if np.issubdtype(df.dtypes[k], np.number):
- # detect and convert integer-only columns
- if (df[column] % 1 == 0).all():
- df[column] = df[column].astype(np.integer)
-
- # now create the appropriate fieldtype
- if np.issubdtype(df.dtypes[k], np.floating):
- sw.field(column, 'N', decimal=5)
- else:
- sw.field(column, 'I', decimal=0)
- else:
- sw.field(column)
-
- # add records to dbf
- for record in df.itertuples():
- sw.record(*record[1:]) # drop first tuple element (=index)
-
- sw.save(filename)
+from shapely.geometry import LineString
+from geopandas import GeoDataFrame
def match_vertices_and_edges(vertices, edges, vertex_cols=('Vertex1', 'Vertex2')):
- """Adds unique IDs to vertices and corresponding edges.
+ """Add columns to edges DataFrame, with the corresponding vertex IDS per edge.
Identifies, which nodes coincide with the endpoints of edges and creates
matching IDs for matching points, thus creating a node-edge graph whose
@@ -146,25 +15,30 @@ def match_vertices_and_edges(vertices, edges, vertex_cols=('Vertex1', 'Vertex2')
is 'Vertex1' and 'Vertex2'.
Args:
- vertices: pandas DataFrame with geometry column of type Point
+ vertices (DataFrame): pandas DataFrame with geometry column of type Point
edges: pandas DataFrame with geometry column of type LineString
- vertex_cols: tuple of 2 strings for the IDs numbers
+ vertex_cols (indexable): tuple of 2 strings for the IDs numbers
Returns:
- Nothing, the mathing IDs are added to the columns vertex_cols in
+ Nothing, the matching IDs are added to the columns vertex_cols in
argument edges
+
+ Note:
+ Although different column names than ``Vertex1`` and ``Vertex2`` are possible,
+ if you want to work with ``rivus.main`` then you should use these as these
+ are pretty much hart coded into the other functions...
"""
vertex_indices = []
for e, line in enumerate(edges.geometry):
edge_endpoints = []
for k, vertex in enumerate(vertices.geometry):
- # change start -------
- # if line.touches(vertex) or line.intersects(vertex):
- # edge_endpoints.append(vertices.index[k])
- if vertex.buffer(0.001).intersects(line) or line.intersects(vertex):
+ if line.touches(vertex) or line.intersects(vertex):
edge_endpoints.append(vertices.index[k])
- # change stop -------
+ # For imperfectly drawn vector layers: TODO CRS-relevant radius
+ # if (vertex.buffer(0.001).intersects(line) or
+ # line.intersects(vertex)):
+ # edge_endpoints.append(vertices.index[k])
if len(edge_endpoints) == 0:
warnings.warn("edge " + str(e) +
@@ -185,15 +59,15 @@ def find_closest_edge(polygons, edges, to_attr='index', column='nearest'):
"""Find closest edge for centroid of polygons.
Args:
- polygons: a pandas DataFrame with geometry column of Polygons
- edges: a pandas DataFrame with geometry column of LineStrings
- to_attr: a column name in DataFrame edges (default: index)
- column: a column name to be added/overwrite in DataFrame polygons with
- the value of column to_attr from the nearest edge in edges
+ polygons (GeoDataFrame): GeoDataFrame with Polygons
+ edges (GeoDataFrame): GeoDataFrame with LineStrings
+ to_attr (str): a column name in GeoDataFrame edges (default: index)
+ column (str): a column name to be added/overwrite in ``polygons`` with
+ the value of column ``to_attr`` from the nearest edge in edges
Returns:
a list of LineStrings connecting polygons' centroids with the nearest
- point in in edges. Side effect: polygons recieves new column with the
+ point in in edges. Side effect: polygons receives new column with the
attribute value of nearest edge. Warning: if column exists, it is
overwritten.
"""
@@ -208,27 +82,30 @@ def find_closest_edge(polygons, edges, to_attr='index', column='nearest'):
nearest_point = shapelytools.project_point_to_object(
centroid, nearest_edge)
- connecting_lines.append(LineString(tuple(centroid.coords) +
- tuple(nearest_point.coords)))
+ connecting_lines.append(
+ LineString(tuple(centroid.coords) + tuple(nearest_point.coords)))
nearest_indices.append(edges[to_attr][nearest_index])
polygons[column] = pd.Series(nearest_indices, index=polygons.index)
- return pd.DataFrame({'geometry': connecting_lines})
+ return GeoDataFrame(geometry=connecting_lines, crs=polygons.crs)
-def bounds(df):
- """Return a DataFrame of minx, miny, maxx, maxy of each geometry."""
- bounds = np.array([geom.bounds for geom in df.geometry])
- return pd.DataFrame(bounds,
- columns=['minx', 'miny', 'maxx', 'maxy'],
- index=df.index)
+def total_bounds(gdf):
+ """Return bounding box (minx, miny, maxx, maxy) of all geometries.
+ Parameters
+ ----------
+ gdf : GeoDataFrame
+ Any GeoDataFrame, in rivus usually vertex or edge.
-def total_bounds(df):
- """Return bounding box (minx, miny, maxx, maxy) of all geometries. """
- b = bounds(df)
+ Returns
+ -------
+ Tuple
+ global (minx, miny, maxx, maxy)
+ """
+ b = gdf.bounds
return (b['minx'].min(),
b['miny'].min(),
b['maxx'].max(),
diff --git a/rivus/utils/pandaspyomo.py b/rivus/utils/pandaspyomo.py
deleted file mode 100644
index 7cbe129..0000000
--- a/rivus/utils/pandaspyomo.py
+++ /dev/null
@@ -1,208 +0,0 @@
-""" pandaspyomo: read data from coopr.pyomo models to pandas DataFrames
-
-Pyomo is a GAMS-like model description language for mathematical
-optimization problems. This module provides functions to read data from
-Pyomo model instances and result objects. Use list_entities to get a list
-of all entities (sets, params, variables, objectives or constraints) inside a
-pyomo instance, before get its contents by get_entity (or get_entities).
-
-Usage:
- import pandaspyomo as pdpo
- pdpo.list_entities(instance, 'var')
- [('EprOut', ['time', 'process', 'commodity', 'commodity']), ...
- ('EprIn', ['time', 'process', 'commodity', 'commodity'])]
- epr = pdpo.get_entities(instance, ['EprOut', 'EprInt'])
- ...
-
-"""
-
-import coopr.pyomo as pyomo
-import pandas as pd
-
-def get_entity(instance, name):
- """ Return a DataFrame for an entity in model instance.
-
- Args:
- instance: a Pyomo ConcreteModel instance
- name: name of a Set, Param, Var, Constraint or Objective
-
- Returns:
- a single-columned Pandas DataFrame with domain as index
- """
-
- # retrieve entity, its type and its onset names
- entity = instance.__getattribute__(name)
- labels = _get_onset_names(entity)
-
- # extract values
- if isinstance(entity, pyomo.Set):
- # Pyomo sets don't have values, only elements
- results = pd.DataFrame([(v, 1) for v in entity.value])
-
- # for unconstrained sets, the column label is identical to their index
- # hence, make index equal to entity name and append underscore to name
- # (=the later column title) to preserve identical index names for both
- # unconstrained supersets
- if not labels:
- labels = [name]
- name = name+'_'
-
- elif isinstance(entity, pyomo.Param):
- if entity.dim() > 1:
- results = pd.DataFrame([v[0]+(v[1],) for v in entity.iteritems()])
- else:
- results = pd.DataFrame(entity.iteritems())
- else:
- # create DataFrame
- if entity.dim() > 1:
- # concatenate index tuples with value if entity has
- # multidimensional indices v[0]
- results = pd.DataFrame(
- [v[0]+(v[1].value,) for v in entity.iteritems()])
- else:
- # otherwise, create tuple from scalar index v[0]
- results = pd.DataFrame(
- [(v[0], v[1].value) for v in entity.iteritems()])
-
- # check for duplicate onset names and append one to several "_" to make
- # them unique, e.g. ['sit', 'sit', 'com'] becomes ['sit', 'sit_', 'com']
- for k, label in enumerate(labels):
- if label in labels[:k]:
- labels[k] = labels[k] + "_"
-
- if not results.empty:
- # name columns according to labels + entity name
- results.columns = labels + [name]
- results.set_index(labels, inplace=True)
-
- return results
-
-
-def get_entities(instance, names):
- """ Return one DataFrame with entities in columns and a common index.
-
- Works only on entities that share a common domain (set or set_tuple), which
- is used as index of the returned DataFrame.
-
- Args:
- instance: a Pyomo ConcreteModel instance
- names: list of entity names (as returned by list_entities)
-
- Returns:
- a Pandas DataFrame with entities as columns and domains as index
- """
-
- df = pd.DataFrame()
- for name in names:
- other = get_entity(instance, name)
-
- if df.empty:
- df = other
- else:
- index_names_before = df.index.names
-
- df = df.join(other, how='outer')
-
- if index_names_before != df.index.names:
- df.index.names = index_names_before
-
- return df
-
-
-def list_entities(instance, entity_type):
- """ Return list of sets, params, variables, constraints or objectives
-
- Args:
- instance: a Pyomo ConcreteModel object
- entity_type: "set", "par", "var", "con" or "obj"
-
- Returns:
- DataFrame of entities
-
- Example:
- >>> data = read_excel('mimo-example.xlsx')
- >>> model = create_model(data, range(1,25))
- >>> list_entities(model, 'obj') #doctest: +NORMALIZE_WHITESPACE
- Description Domain
- Name
- obj minimize(cost = sum of all cost types) []
-
- """
-
- # helper function to discern entities by type
- def filter_by_type(entity, entity_type):
- if entity_type == 'set':
- return isinstance(entity, pyomo.Set) and not entity.virtual
- elif entity_type == 'par':
- return isinstance(entity, pyomo.Param)
- elif entity_type == 'var':
- return isinstance(entity, pyomo.Var)
- elif entity_type == 'con':
- return isinstance(entity, pyomo.Constraint)
- elif entity_type == 'obj':
- return isinstance(entity, pyomo.Objective)
- else:
- raise ValueError("Unknown entity_type '{}'".format(entity_type))
-
- # iterate through all model components and keep only
- iter_entities = instance.__dict__.iteritems()
- entities = sorted(
- (name, entity.doc, _get_onset_names(entity))
- for (name, entity) in iter_entities
- if filter_by_type(entity, entity_type))
-
- # if something was found, wrap tuples in DataFrame, otherwise return empty
- if entities:
- entities = pd.DataFrame(entities,
- columns=['Name', 'Description', 'Domain'])
- entities.set_index('Name', inplace=True)
- else:
- entities = pd.DataFrame()
- return entities
-
-
-def _get_onset_names(entity):
- """
- Example:
- >>> data = read_excel('mimo-example.xlsx')
- >>> model = create_model(data, range(1,25))
- >>> _get_onset_names(model.e_co_stock)
- ['t', 'sit', 'com', 'com_type']
- """
- # get column titles for entities from domain set names
- labels = []
-
- if isinstance(entity, pyomo.Set):
- if entity.dimen > 1:
- # N-dimensional set tuples, possibly with nested set tuples within
- if entity.domain:
- domains = entity.domain.set_tuple
- else:
- domains = entity.set_tuple
-
- for domain_set in domains:
- labels.extend(_get_onset_names(domain_set))
-
- elif entity.dimen == 1:
- if entity.domain:
- # 1D subset; add domain name
- labels.append(entity.domain.name)
- else:
- # unrestricted set; add entity name
- labels.append(entity.name)
- else:
- # no domain, so no labels needed
- pass
-
- elif isinstance(entity, (pyomo.Param, pyomo.Var, pyomo.Constraint,
- pyomo.Objective)):
- if entity.dim() > 0 and entity._index:
- labels = _get_onset_names(entity._index)
- else:
- # zero dimensions, so no onset labels
- pass
-
- else:
- raise ValueError("Unknown entity type!")
-
- return labels
\ No newline at end of file
diff --git a/rivus/utils/pyomotools.py b/rivus/utils/pyomotools.py
deleted file mode 100644
index e28de43..0000000
--- a/rivus/utils/pyomotools.py
+++ /dev/null
@@ -1,77 +0,0 @@
-""" pyomotools: common helper functions for pyomo model creation """
-from datetime import datetime
-import pandas as pd
-import xlrd
-
-__all__ = ["now", "read_xls"]
-
-def now(mydateformat='%Y%m%dT%H%M%S'):
- """ Return current datetime as string.
-
- Just a shorthand to abbreviate the common task to obtain the current
- datetime as a string, e.g. for result versioning.
-
- Args:
- mydateformat: optional format string (default: '%Y%m%dT%H%M%S')
-
- Returns:
- datetime.now(), formated to string with argument mydateformat, e.g.
- YYYYMMDDThhmmss ==> 20131007H123456
- """
- return datetime.now().strftime(mydateformat)
-
-
-def read_xls(filename, sheets=[]):
- """ Convert Excel file to dict of pandas DataFrames.
-
- Parses all spreadsheets within an Excel file using pandas.ExcelFile.parse,
- if its top left cell is not empty. The first row is expected to contain
- column titles. Titles starting with uppercase lettres are used as index
- columns in the resulting DataFrame. Here is a short example summarizing
- these specifications:
-
- Process CoIn CoOut | cap eff ... avail
- -------------------------------------------
- PP Coal Elec | 100 0.90 ... 24
- WT Wind Elec | 300 0.95 ... 10
- PV Solar Elec | 200 0.92 ... 8
-
- A spreadsheet is skipped if a) it is completely empty or b) has an empty
- first row.
-
- Args:
- filename: an Excel spreadsheet filename
-
- Returns:
- dict of pandas DataFrames with sheet names as keys
- """
-
- dfs = {}
- xls = pd.ExcelFile(filename)
- for sheet in xls.book.sheets():
- # skip sheet if list of sheets was specified
- if sheets and sheet.name not in sheets:
- continue
-
- # extract the sheet's first row to check for emptiness
- first_row = sheet.row_slice(0)
-
- # skip a spreadsheet if completely empty or its first cell is blank
- if not first_row \
- or first_row[0].ctype in (xlrd.XL_CELL_BLANK, xlrd.XL_CELL_EMPTY):
- continue
-
- # otherwise determine column numbers of titles starting with an
- # uppercase lettre while skipping empty columns
- uppercase_columns = [k for k, column_title in enumerate(first_row)
- if column_title.value
- and column_title.value[0].isupper()]
-
- # parse those columns to a pandas DataFrame
- df = xls.parse(sheet.name, index_col=uppercase_columns)
-
- # and prune any columns with only NaN values
- # these are mainly empty columns
- dfs[sheet.name] = df.dropna(axis=1, how='all')
-
- return dfs
\ No newline at end of file
diff --git a/rivus/utils/shapelytools.py b/rivus/utils/shapelytools.py
index 460f0bc..805757c 100644
--- a/rivus/utils/shapelytools.py
+++ b/rivus/utils/shapelytools.py
@@ -1,161 +1,226 @@
-from shapely.geometry import (box, LineString, MultiLineString, MultiPoint,
- Point, Polygon)
+from shapely.geometry import (box, LineString, MultiLineString, MultiPoint,
+ Point, Polygon)
import shapely.ops
+
def endpoints_from_lines(lines):
- """Return list of terminal points from list of LineStrings."""
-
+ """Return list of terminal points from list of LineStrings.
+ Where in case of common end points, no duplicate is generated.
+
+ Parameters
+ ----------
+ lines : list(LineStrings)
+
+ Returns
+ -------
+ list of Shapely.Point
+ """
+
all_points = []
for line in lines:
- for i in [0, -1]: # start and end point
+ for i in [0, -1]: # start and end point
all_points.append(line.coords[i])
-
+
unique_points = set(all_points)
-
+
return [Point(p) for p in unique_points]
-
+
+
def vertices_from_lines(lines):
- """Return list of unique vertices from list of LineStrings."""
-
+ """Return list of unique vertices from list of LineStrings.
+
+ Parameters
+ ----------
+ lines : list(LineStrings)
+
+ Returns
+ -------
+ list of Shapely.Point
+ """
+
vertices = []
for line in lines:
vertices.extend(list(line.coords))
return [Point(p) for p in set(vertices)]
-#
-def snapping_vertexis_from_lines(lines, closeness_limit):
+def snapping_vertices_from_lines(lines, closeness_limit):
"""" Return a list of edge endpoints, where close, but not identical
- points are filtered
+ points are unified
+
+ Note
+ -----
+ When working with projected coordinates, the unit of closeness_limit is meters.
+ But do not forget, that the majority of functions in rivus returns
+ data in WGS84 as CRS. (LatLon, not projected coordinates.)
+ In this case, the unit of closeness_limit is a good question. But the meaning
+ of meters would be well beyond the decimal point.
- Args:
- lines: list of LineStrings
- closeness_limit: critical distance for selecting uniqe verteces
+ Suggestion: use projected crs for this calculation.
+
+ Parameters
+ ----------
+ lines
+ list of LineStrings
+ closeness_limit
+ critical distance for selecting unique vertices.
+ From which vertices will be united with their neighbours in their proximity.
+
+ Returns
+ -------
+ TYPE
+ Description
"""
uni_verts = []
for line in lines:
for end in line.coords:
if len(uni_verts):
- inter = Point(end).buffer(closeness_limit).intersection(MultiPoint(uni_verts))
+ inter = (Point(end).buffer(closeness_limit)
+ .intersection(MultiPoint(uni_verts)))
if inter.type != 'Point':
uni_verts.append(end)
else:
uni_verts.append(end)
return [Point(p) for p in uni_verts]
-#
def prune_short_lines(lines, min_length):
"""Remove lines from a LineString DataFrame shorter than min_length.
-
+
Deletes all lines from a list of LineStrings or a MultiLineString
- that have a total length of less than min_length. Vertices of touching
+ that have a total length of less than min_length. Vertices of touching
lines are contracted towards the centroid of the removed line.
-
- Args:
- lines: list of LineStrings or a MultiLineString
- min_length: minimum length of a single LineString to be preserved
-
- Returns:
- the pruned pandas DataFrame
- """
- pruned_lines = [line for line in lines] # converts MultiLineString to list
+
+ Parameters
+ ----------
+ lines
+ list of LineStrings or a MultiLineString
+ min_length
+ minimum length of a single LineString to be preserved
+
+ Returns
+ -------
+ list
+ including not pruned lines
+ """
+ pruned_lines = [line for line in lines] # converts MultiLineString to list
to_prune = []
-
+
for i, line in enumerate(pruned_lines):
if line.length < min_length:
to_prune.append(i)
for n in neighbors(pruned_lines, line):
contact_point = line.intersection(pruned_lines[n])
- pruned_lines[n] = bend_towards(pruned_lines[n],
+ pruned_lines[n] = bend_towards(pruned_lines[n],
where=contact_point,
to=line.centroid)
-
- return [line for i, line in enumerate(pruned_lines) if i not in to_prune]
+
+ return [line for i, line in enumerate(pruned_lines) if i not in to_prune]
def neighbors(lines, of):
"""Find the indices in a list of LineStrings that touch a given LineString.
-
- Args:
- lines: list of LineStrings in which to search for neighbors
- of: the LineString which must be touched
-
- Returns:
- list of indices, so that all lines[indices] touch the LineString of
+
+ Parameters
+ ----------
+ lines
+ list of LineStrings in which to search for neighbours
+ of
+ the LineString which must be touched
+
+ Returns
+ -------
+ list of indices, so that all lines[indices] touch the LineString of
"""
return [k for k, line in enumerate(lines) if line.touches(of)]
-
+
def bend_towards(line, where, to):
- """Move the point where along a line to the point at location to.
-
- Args:
- line: a LineString
- where: a point ON the line (not necessarily a vertex)
- to: a point NOT on the line where the nearest vertex will be moved to
-
- Returns:
- the modified (bent) line
+ """Move the point ``where`` on the ``line`` to the point ``to``.
+
+ Parameters
+ ----------
+ line : Shapely.LineString
+ a LineString
+ where : Shapely.Point
+ a point ON the line (not necessarily a vertex)
+ to : Shapely.Point
+ a point NOT on the line where the nearest vertex will be moved to
+
+ Returns
+ -------
+ LineString
+ the modified (bent) line.
+
+ Raises
+ ------
+ ValueError
+ ``line`` does not contain the point ``where``.
"""
-
+
if not line.contains(where) and not line.touches(where):
raise ValueError('line does not contain the point where.')
-
+
coords = line.coords[:]
- # easy case: where is (within numeric precision) a vertex of line
+ # easy case: ``where`` is (within numeric precision) a vertex of line
for k, vertex in enumerate(coords):
if where.almost_equals(Point(vertex)):
# move coordinates of the vertex to destination
coords[k] = to.coords[0]
return LineString(coords)
-
- # hard case: where lies between vertices of line, so
+
+ # hard case: ``where`` lies between vertices of line, so
# find nearest vertex and move that one to point to
- _, min_k = min((where.distance(Point(vertex)), k)
- for k, vertex in enumerate(coords))
+ _, min_k = min((where.distance(Point(vertex)), k)
+ for k, vertex in enumerate(coords))
coords[min_k] = to.coords[0]
return LineString(coords)
def snappy_endings(lines, max_distance):
"""Snap endpoints of lines together if they are at most max_length apart.
-
- Args:
- lines: a list of LineStrings or a MultiLineString
- max_distance: maximum distance two endpoints may be joined together
+
+ Parameters
+ ----------
+ lines
+ a list of LineStrings or a MultiLineString
+ max_distance
+ maximum distance two endpoints may be joined together
+
+ Returns
+ -------
+ list(LineString)
+ snapped lines
"""
-
+
# initialize snapped lines with list of original lines
# snapping points is a MultiPoint object of all vertices
snapped_lines = [line for line in lines]
- #
+ # If the ends of lines truly meet at an intersection
snapping_points = vertices_from_lines(snapped_lines)
- #
- #
- # snapping_points = snapping_vertexis_from_lines(snapped_lines, max_distance)
- #
-
+ # If the ends of (self drawn) lines may not snap together perfectly
+ # snapping_points = snapping_vertices_from_lines(snapped_lines, max_distance)
+
# isolated endpoints are going to snap to the closest vertex
- isolated_endpoints = find_isolated_endpoints(snapped_lines)
-
+ isolated_endpoints = find_isolated_endpoints(snapped_lines)
+
# only move isolated endpoints, one by one
for endpoint in isolated_endpoints:
# find all vertices within a radius of max_distance as possible
- target = nearest_neighbor_within(snapping_points, endpoint,
+ target = nearest_neighbor_within(snapping_points, endpoint,
max_distance)
-
+
# do nothing if no target point to snap to is found
if not target:
- continue
-
- # find the LineString to modify within snapped_lines and update it
+ continue
+
+ # find the LineString to modify within snapped_lines and update it
for i, snapped_line in enumerate(snapped_lines):
if endpoint.touches(snapped_line):
- snapped_lines[i] = bend_towards(snapped_line, where=endpoint,
+ snapped_lines[i] = bend_towards(snapped_line, where=endpoint,
to=target)
break
-
+
# also update the corresponding snapping_points
for i, snapping_point in enumerate(snapping_points):
if endpoint.equals(snapping_point):
@@ -166,94 +231,116 @@ def snappy_endings(lines, max_distance):
snapped_lines = [s for s in snapped_lines if s.length > 0]
return snapped_lines
-
-
+
+
def nearest_neighbor_within(others, point, max_distance):
"""Find nearest point among others up to a maximum distance.
-
- Args:
- others: a list of Points or a MultiPoint
- point: a Point
- max_distance: maximum distance to search for the nearest neighbor
-
- Returns:
+
+ Parameters
+ ----------
+ others
+ a list of Points or a MultiPoint
+ point
+ a Point
+ max_distance
+ maximum distance to search for the nearest neighbour
+
+ Returns
+ -------
+ Shapely.Point
A shapely Point if one is within max_distance, None otherwise
"""
search_region = point.buffer(max_distance)
interesting_points = search_region.intersection(MultiPoint(others))
-
+
if not interesting_points:
closest_point = None
elif isinstance(interesting_points, Point):
closest_point = interesting_points
- else:
+ else:
distances = [point.distance(ip) for ip in interesting_points
if point.distance(ip) > 0]
closest_point = interesting_points[distances.index(min(distances))]
-
+
return closest_point
def find_isolated_endpoints(lines):
"""Find endpoints of lines that don't touch another line.
-
- Args:
- lines: a list of LineStrings or a MultiLineString
-
- Returns:
- A list of line end Points that don't touch any other line of lines
+
+ Parameters
+ ----------
+ lines
+ a list of LineStrings or a MultiLineString
+
+ Returns
+ -------
+ A list of line end Points that don't touch any other line of lines
"""
-
+
isolated_endpoints = []
for i, line in enumerate(lines):
other_lines = lines[:i] + lines[i+1:]
- for q in [0,-1]:
+ for q in [0, -1]:
endpoint = Point(line.coords[q])
- if any(endpoint.touches(another_line)
+ if any(endpoint.touches(another_line)
for another_line in other_lines):
continue
else:
isolated_endpoints.append(endpoint)
return isolated_endpoints
-
+
+
def closest_object(geometries, point):
"""Find the nearest geometry among a list, measured from fixed point.
-
- Args:
- geometries: a list of shapely geometry objects
- point: a shapely Point
-
- Returns:
- Tuple (geom, min_dist, min_index) of the geometry with minimum distance
- to point, its distance min_dist and the list index of geom, so that
- geom = geometries[min_index].
- """
- min_dist, min_index = min((point.distance(geom), k)
+
+ Parameters
+ ----------
+ geometries
+ a list of shapely geometry objects
+ point
+ a shapely Point
+
+ Returns
+ -------
+ Tuple (geom, min_dist, min_index) of the geometry with minimum distance
+ to point, its distance min_dist and the list index of geom, so that
+ geom = geometries[min_index].
+ """
+ min_dist, min_index = min((point.distance(geom), k)
for (k, geom) in enumerate(geometries))
-
+
return geometries[min_index], min_dist, min_index
-
-
+
+
def project_point_to_line(point, line_start, line_end):
"""Find nearest point on a straight line, measured from given point.
-
- Args:
- point: a shapely Point object
- line_start: the line starting point as a shapely Point
- line_end: the line end point as a shapely Point
-
- Returns:
- a shapely Point that lies on the straight line closest to point
-
- Source: http://gis.stackexchange.com/a/438/19627
+
+ Parameters
+ ----------
+ point
+ a shapely Point object
+ line_start
+ the line starting point as a shapely Point
+ line_end
+ the line end point as a shapely Point
+
+ Returns
+ -------
+ Shapely.Point
+ The projected point
+
+ References
+ -----------
+ `Source `_
"""
line_magnitude = line_start.distance(line_end)
-
+
u = ((point.x - line_start.x) * (line_end.x - line_start.x) +
(point.y - line_start.y) * (line_end.y - line_start.y)) \
- / (line_magnitude ** 2)
+ / (line_magnitude ** 2)
- # closest point does not fall within the line segment,
+ # closest point does not fall within the line segment,
# take the shorter distance to an endpoint
if u < 0.00001 or u > 1:
ix = point.distance(line_start)
@@ -266,25 +353,31 @@ def project_point_to_line(point, line_start, line_end):
ix = line_start.x + u * (line_end.x - line_start.x)
iy = line_start.y + u * (line_end.y - line_start.y)
return Point([ix, iy])
-
-def pairs(lst):
+
+
+def pairs(a_list):
"""Iterate over a list in overlapping pairs.
-
- Args:
- lst: an iterable/list
-
- Returns:
- Yields a pair of consecutive elements (lst[k], lst[k+1]) of lst. Last
- call yields (lst[-2], lst[-1]).
-
- Example:
- lst = [4, 7, 11, 2]
- pairs(lst) yields (4, 7), (7, 11), (11, 2)
-
- Source:
- http://stackoverflow.com/questions/1257413/1257446#1257446
+
+ Parameters
+ ----------
+ a_list
+ an iterable/list
+
+ Example
+ -------
+ a_list = [4, 7, 11, 2]
+ pairs(a_list) yields (4, 7), (7, 11), (11, 2)
+
+ References
+ -----------
+ `Source `_
+
+ Yields
+ ------
+ Yields a pair of consecutive elements (a_list[k], a_list[k+1]) of a_list.
+ Last call yields (a_list[-2], a_list[-1]).
"""
- i = iter(lst)
+ i = iter(a_list)
prev = next(i)
for item in i:
yield prev, item
@@ -293,58 +386,72 @@ def pairs(lst):
def project_point_to_object(point, geometry):
"""Find nearest point in geometry, measured from given point.
-
- Args:
- point: a shapely Point
- geometry: a shapely geometry object (LineString, Polygon)
-
- Returns:
- a shapely Point that lies on geometry closest to point
+
+ Parameters
+ ----------
+ point
+ a shapely Point
+ geometry
+ a shapely geometry object (LineString, Polygon)
+
+ Returns
+ -------
+ Shapely.Point
+ that lies on geometry closest to point
+
+ Raises
+ ------
+ NotImplementedError
+ project_point_to_object not implemented for geometry type
"""
nearest_point = None
min_dist = float("inf")
-
+
if isinstance(geometry, Polygon):
for seg_start, seg_end in pairs(list(geometry.exterior.coords)):
line_start = Point(seg_start)
line_end = Point(seg_end)
-
+
intersection_point = project_point_to_line(point, line_start, line_end)
- cur_dist = point.distance(intersection_point)
-
+ cur_dist = point.distance(intersection_point)
+
if cur_dist < min_dist:
min_dist = cur_dist
nearest_point = intersection_point
-
+
elif isinstance(geometry, LineString):
for seg_start, seg_end in pairs(list(geometry.coords)):
line_start = Point(seg_start)
line_end = Point(seg_end)
-
+
intersection_point = project_point_to_line(point, line_start, line_end)
- cur_dist = point.distance(intersection_point)
-
+ cur_dist = point.distance(intersection_point)
+
if cur_dist < min_dist:
min_dist = cur_dist
nearest_point = intersection_point
else:
- raise NotImplementedError("project_point_to_object not implemented for"+
+ raise NotImplementedError("project_point_to_object not implemented for" +
" geometry type '" + geometry.type + "'.")
return nearest_point
-
+
def one_linestring_per_intersection(lines):
- """ Move line endpoints to intersections of line segments.
-
+ """Move line endpoints to intersections of line segments.
+
Given a list of touching or possibly intersecting LineStrings, return a
- list LineStrings that have their endpoints at all crossings and
+ list of LineStrings that have their endpoints at all crossings and
intersecting points and ONLY there.
-
- Args:
+
+ Parameters
+ ----------
+ lines : list(LineString)
a list of LineStrings or a MultiLineString
-
- Returns:
- a list of LineStrings
+ Possible complex collection of line strings.
+
+ Returns
+ -------
+ list of LineStrings
"""
lines_merged = shapely.ops.linemerge(lines)
@@ -360,16 +467,19 @@ def one_linestring_per_intersection(lines):
def linemerge(linestrings_or_multilinestrings):
- """ Merge list of LineStrings and/or MultiLineStrings.
-
+ """Merge list of LineStrings and/or MultiLineStrings.
+
Given a list of LineStrings and possibly MultiLineStrings, merge all of
them to a single MultiLineString.
-
- Args:
+
+ Parameters
+ ----------
+ linestrings_or_multilinestrings
list of LineStrings and/or MultiLineStrings
-
- Returns:
- a merged LineString or MultiLineString
+
+ Returns
+ -------
+ a merged LineString or MultiLineString
"""
lines = []
for line in linestrings_or_multilinestrings:
@@ -379,6 +489,5 @@ def linemerge(linestrings_or_multilinestrings):
else:
# line is a line, so simply append it
lines.append(line)
-
+
return shapely.ops.linemerge(lines)
-
diff --git a/rivus/utils/shptools.py b/rivus/utils/shptools.py
deleted file mode 100644
index 1b8b34c..0000000
--- a/rivus/utils/shptools.py
+++ /dev/null
@@ -1,218 +0,0 @@
-import itertools
-import shapefile
-from shapely.geometry import Polygon, MultiLineString, LineString, Point
-import pdb
-
-def read_shp(filename):
- """Read contents of a shapefile to a shapely geometry object.
-
- Usage:
- geometries = read_shp(filename)
- (geometries, records, fields) = read_shp(filename)
-
- Arguments:
- filename shapefile name
-
- Returns:
- geometries list of shapely geometries (Polygon, LineString, ...)
- records list of records (list of values), one per geometry
- fields list of fieldnames of a record
- """
- sr = shapefile.Reader(filename)
-
- if sr.shapeType == shapefile.POLYGON:
- shapes = sr.shapes()
- geometries = [Polygon(shape.points) for shape in shapes]
-
- fields = sr.fields[:]
- if fields[0][0] == 'DeletionFlag':
- fields.pop(0)
- fields = [field[0] for field in fields] # extract field name only
-
- records = []
- for record in sr.records():
- for i, value in enumerate(record):
- try:
- record[i] = float(value) # convert record values to numeric...
- except ValueError:
- pass # ... if possible
-
- records.append(record)
-
- return (geometries, records, fields)
-
- elif sr.shapeType == shapefile.POLYLINE:
- shapes = sr.shapes()
- geometries = [LineString(shape.points) for shape in shapes]
-
- fields = sr.fields[:] # [:] = duplicate field list
- if fields[0][0] == 'DeletionFlag':
- fields.pop(0)
- fields = [field[0] for field in fields] # extract field name only
-
- records = []
- for record in sr.records():
- for i, value in enumerate(record):
- try:
- record[i] = float(value) # convert record values to numeric...
- except ValueError:
- pass # ... if possible
-
- records.append(record)
-
- return (geometries, records, fields)
-
-
- elif sr.shapeType == shapefile.MULTIPOINT:
- raise NotImplementedError
-
- else:
- raise NotImplementedError
-
-
-
-
-
-def write_shp(filename, geometry, records=[], fields=[]):
- """Write a single shapely MultiLineString or Polygon to a shapefile.
-
- Argument geometry may also be a list of LineString objects. In that case,
- each entry may have associated data in the optional records list. If
- records is given, fields is the list of fieldnames for the value list in
- each record.
-
- Usage:
- write_shp(filename, geometry, records=[], fields=[])
-
- Arguments:
- filename filename of shapefile
- geometry a shapely geometry (MultiLineString, Polygon, ...)
- records optional list of list of values, one per geometry
- fields optional (implied by records) list of fieldnames
- """
-
- # SINGLE MULTILINESTRING
- if isinstance(geometry, MultiLineString):
- sw = shapefile.Writer(shapefile.POLYLINE)
-
- # fields
- sw.field("length")
- sw.field("start-x")
- sw.field("start-y")
- sw.field("end-x")
- sw.field("end-y")
- sw.field("npoints")
-
- # geometry and record
- for line in geometry:
- sw.line([list(line.coords)])
- sw.record(line.length,
- line.coords[0][0],
- line.coords[0][1],
- line.coords[-1][0],
- line.coords[-1][1],
- len(line.coords))
- sw.save(filename)
-
- # SINGLE POLYGON
- elif isinstance(geometry, Polygon):
- # data
- parts = [list(geometry.exterior.coords)]
- parts.extend(list(interior.coords) for interior in geometry.interiors)
-
- sw = shapefile.Writer(shapefile.POLYGON)
- sw.field("area")
- sw.poly(parts)
- sw.record(geometry.area)
- sw.save(filename)
-
- # LISTS
- elif isinstance(geometry, list):
-
- # check for fields and records
- if (fields and not records) or (not fields and records):
- raise ValueError('Arguments records and fields must both be provided.')
-
- if records and (len(fields) != len(records[0])):
- raise ValueError('Length of fields and records do not match.')
-
- # derive field types based on record values
- field_type = {}
- precision = {}
- for i, field in enumerate(fields):
- if all(type(record[i]) in [int, long] for record in records):
- field_type[field] = 'N' # integers
- precision[field] = 0
- elif all(type(record[i]) in [int, long, float] for record in records):
- field_type[field] = 'N' # numeric
- precision[field] = 5
- else:
- field_type[field] = 'C' # string (characters)
- precision[field] = 0
-
-
- # LIST OF LINESTRINGS
- if isinstance(geometry[0], LineString):
- sw = shapefile.Writer(shapefile.POLYLINE)
-
- # fields
- for field in fields:
- sw.field(field, field_type[field], decimal=precision[field])
-
- if not fields:
- # add dummy length field and prepare records for it
- records = [[line.length] for line in geometry]
- sw.field('length', 'N', decimal=5)
-
- # geometry and records
- for line, record in itertools.izip(geometry, records):
- sw.line([list(line.coords)])
- sw.record(*record)
-
- sw.save(filename)
-
- # LIST OF POLYGONS
- elif isinstance(geometry[0], Polygon):
- sw = shapefile.Writer(shapefile.POLYGON)
-
- # fields
- for field in fields:
- sw.field(field, field_type[field], decimal=precision[field])
-
- # geometry and records
- for polygon, record in itertools.izip(geometry, records):
- sw.poly([list(polygon.exterior.coords)])
- sw.record(*record)
-
- sw.save(filename)
-
- # LIST OF POINTS
- elif isinstance(geometry[0], Point):
- sw = shapefile.Writer(shapefile.POINT)
-
- # fields
- for field in fields:
- sw.field(field, field_type[field], decimal=precision[field])
-
- if not fields:
- # add dummy length field and prepare records for it
- records = [[point.x, point.y] for point in geometry]
- sw.field('x', 'N', decimal=5)
- sw.field('y', 'N', decimal=5)
-
- # shapes and records
- for point, record in itertools.izip(geometry, records):
- sw.point(point.x, point.y)
- sw.record(*record)
-
- # save
- sw.save(filename)
-
-
- else:
- raise NotImplementedError
-
- else:
- raise NotImplementedError
-
-
diff --git a/rivus/utils/skeletrontools.py b/rivus/utils/skeletrontools.py
index 639f5d6..0870af0 100644
--- a/rivus/utils/skeletrontools.py
+++ b/rivus/utils/skeletrontools.py
@@ -1,7 +1,5 @@
-from shapely.geometry import Polygon, LineString, Point, box
-
+from shapely.geometry import Polygon, box
import Skeletron
-import pandashp
import shapely.ops
@@ -25,20 +23,28 @@ def extract_lines_from_graph(graphs):
return lines
-def skeletonize(roads, buffer_length=60,
- dissolve_length=30,
- simplify_length=30,
- buffer_resolution=2,
- psg_length=150):
+def skeletonize(roads, buffer_length=60, dissolve_length=30,
+ simplify_length=30, buffer_resolution=2, psg_length=150):
"""Uses qhull to find simplified road network for given DataFrame of roads.
- Args:
- roads pandashp DataFrame of shapely LINESTRINGs (projected)
- buffer_length optional roughly equivalent to amount of generalization
- dissolve_length optional (def: 30) considerably smaller than buffer_length
- simplify_length optional
- buffer_resolution optional
- psg_length optional (default: 150) Skeletron algorithm length
+ Parameters
+ ----------
+ roads : GeoDataFrame
+ pandashp DataFrame of shapely LINESTRINGs (projected)
+ buffer_length : int, optional
+ optional roughly equivalent to amount of generalization
+ dissolve_length : int, optional
+ optional (def: 30) considerably smaller than buffer_length
+ simplify_length : int, optional
+ meters?
+ buffer_resolution : int, optional
+ TODO
+ psg_length : int, optional
+ optional (default: 150) Skeletron algorithm length
+
+ Returns
+ -------
+ LineString or MultiLineString when lines are not contiguous
"""
# buffer and merge streets
diff --git a/runhaag.py b/runhaag.py
index 3970376..5777563 100644
--- a/runhaag.py
+++ b/runhaag.py
@@ -1,7 +1,7 @@
import matplotlib.pyplot as plt
import os
import pandas as pd
-from rivus.utils import pandashp as pdshp
+import geopandas
from rivus.main import rivus
try:
import pyomo.environ
@@ -66,8 +66,8 @@ def prepare_edge(edge_shapefile, building_shapefile):
# 2. group DataFrame by columns 'nearest' (ID of nearest edge) and 'type'
# (residential, commercial, industrial, other)
# 3. sum by group and unstack, i.e. convert secondary index 'type' to columns
- buildings = pdshp.read_shp(building_shapefile)
- building_type_mapping = {
+ buildings = geopandas.read_file(building_shapefile + '.shp')
+ building_type_mapping = {
'church': 'other',
'farm': 'other',
'hospital': 'residential',
@@ -85,7 +85,7 @@ def prepare_edge(edge_shapefile, building_shapefile):
# 1. read shapefile to DataFrame (with geometry column)
# 2. join DataFrame total_area on index (=ID)
# 3. fill missing values with 0
- edge = pdshp.read_shp(edge_shapefile)
+ edge = geopandas.read_file(edge_shapefile + '.shp')
edge = edge.set_index('Edge')
edge = edge.join(total_area)
edge = edge.fillna(0)
@@ -100,7 +100,7 @@ def run_scenario(scenario):
# prepare input data
data = rivus.read_excel(data_spreadsheet)
- vertex = pdshp.read_shp(vertex_shapefile)
+ vertex = geopandas.read_file(vertex_shapefile + '.shp')
edge = prepare_edge(edge_shapefile, building_shapefile)
# apply scenario function to input data
diff --git a/runhg15.py b/runhg15.py
index 7c5c5c4..e341e77 100644
--- a/runhg15.py
+++ b/runhg15.py
@@ -8,7 +8,6 @@
PYOMO3 = True
import geopandas
import os
-from rivus.utils import pandashp as pdshp
from rivus.main import rivus
from datetime import datetime
@@ -197,7 +196,7 @@ def prepare_edge(edge_shapefile, building_shapefile):
# 1. read shapefile to DataFrame (with geometry column)
# 2. join DataFrame total_area on index (=ID)
# 3. fill missing values with 0
- edge = pdshp.read_shp(edge_shapefile)
+ edge = geopandas.read_file(edge_shapefile + '.shp')
edge = edge.set_index('Edge')
edge = edge.join(total_area)
edge = edge.fillna(0)
@@ -211,7 +210,7 @@ def run_scenario(scenario, result_dir):
# prepare input data
data = rivus.read_excel(data_spreadsheet)
- vertex = pdshp.read_shp(vertex_shapefile)
+ vertex = geopandas.read_file(vertex_shapefile + '.shp')
edge = prepare_edge(edge_shapefile, building_shapefile)
# apply scenario function to input data
diff --git a/runmin.py b/runmin.py
index 4e2dabb..8cbed8a 100644
--- a/runmin.py
+++ b/runmin.py
@@ -8,8 +8,7 @@
PYOMO3 = True
import matplotlib.pyplot as plt
import os
-import pandas as pd
-from rivus.utils import pandashp as pdshp
+import geopandas
from rivus.main import rivus
base_directory = os.path.join('data', 'mnl')
@@ -24,21 +23,21 @@
# 2. group DataFrame by columns 'nearest' (ID of nearest edge) and 'type'
# (residential, commercial, industrial, other)
# 3. sum by group and unstack, i.e. convert secondary index 'type' to columns
-buildings = pdshp.read_shp(building_shapefile)
+buildings = geopandas.read_file(building_shapefile + '.shp')
buildings_grouped = buildings.groupby(['nearest', 'type'])
total_area = buildings_grouped.sum()['total_area'].unstack()
-# load edges (streets) and join with summed areas
+# load edges (streets) and join with summed areas
# 1. read shapefile to DataFrame (with geometry column)
# 2. join DataFrame total_area on index (=ID)
# 3. fill missing values with 0
-edge = pdshp.read_shp(edge_shapefile)
+edge = geopandas.read_file(edge_shapefile + '.shp')
edge = edge.set_index('Edge')
edge = edge.join(total_area)
edge = edge.fillna(0)
# load nodes
-vertex = pdshp.read_shp(vertex_shapefile)
+vertex = geopandas.read_file(vertex_shapefile + '.shp')
# load spreadsheet data
data = rivus.read_excel(data_spreadsheet)
@@ -46,11 +45,11 @@
# create and solve model
prob = rivus.create_model(data, vertex, edge)
if PYOMO3:
- prob = prob.create() # no longer needed in Pyomo 4<
+ prob = prob.create() # no longer needed in Pyomo 4<
solver = SolverFactory('glpk')
result = solver.solve(prob, tee=True)
if PYOMO3:
- prob.load(result) #no longer needed in Pyomo 4<
+ prob.load(result) # no longer needed in Pyomo 4<
# load results
costs, Pmax, Kappa_hub, Kappa_process = rivus.get_constants(prob)
@@ -62,23 +61,23 @@
# create result directory if not existing already
if not os.path.exists(result_dir):
os.makedirs(result_dir)
-
+
rivus.save(prob, os.path.join(result_dir, 'prob.pgz'))
rivus.report(prob, os.path.join(result_dir, 'report.xlsx'))
# plot all caps (and demands if existing)
for com, plot_type in [('Elec', 'caps'), ('Heat', 'caps'), ('Gas', 'caps'),
('Elec', 'peak'), ('Heat', 'peak')]:
-
+
# create plot
- fig = rivus.plot(prob, com, mapscale=False, tick_labels=False,
- plot_demand=(plot_type == 'peak'))
+ fig = rivus.plot(prob, com, mapscale=False, tick_labels=False,
+ plot_demand=(plot_type == 'peak'))
plt.title('')
# save to file
for ext in ['png', 'pdf']:
-
+
# determine figure filename from plot type, commodity and extension
fig_filename = os.path.join(
result_dir, '{}-{}.{}').format(plot_type, com, ext)
- fig.savefig(fig_filename, dpi=300, bbox_inches='tight',
+ fig.savefig(fig_filename, dpi=300, bbox_inches='tight',
transparent=(ext=='pdf'))
diff --git a/runmoosh.py b/runmoosh.py
index 1b97bb7..554613d 100644
--- a/runmoosh.py
+++ b/runmoosh.py
@@ -9,8 +9,8 @@
import matplotlib.pyplot as plt
import os
import pandas as pd
+import geopandas
from operator import itemgetter
-from rivus.utils import pandashp as pdshp
from rivus.main import rivus
base_directory = os.path.join('data', 'moosh')
@@ -44,7 +44,7 @@ def setup_solver(optim):
# 2. group DataFrame by columns 'nearest' (ID of nearest edge) and 'type'
# (residential, commercial, industrial, other)
# 3. sum by group and unstack, i.e. convert secondary index 'type' to columns
-buildings = pdshp.read_shp(building_shapefile)
+buildings = geopandas.read_file(building_shapefile + '.shp')
building_type_mapping = {
'church': 'other',
'farm': 'other',
@@ -64,13 +64,13 @@ def setup_solver(optim):
# 1. read shapefile to DataFrame (with geometry column)
# 2. join DataFrame total_area on index (=ID)
# 3. fill missing values with 0
-edge = pdshp.read_shp(edge_shapefile)
+edge = geopandas.read_file(edge_shapefile + '.shp')
edge = edge.set_index('Edge')
edge = edge.join(total_area)
edge = edge.fillna(0)
# load nodes
-vertex = pdshp.read_shp(vertex_shapefile)
+vertex = geopandas.read_file(vertex_shapefile + '.shp')
# load spreadsheet data
data = rivus.read_excel(data_spreadsheet)