"""
Container objects for working with geometric information
"""
import numpy as np
from ..utils import import_optional_dependency
[docs]class Shape:
"""
Parent class for handling geo interfacing, do not instantiate directly
Parameters
----------
type : str
shapetype string
coordinates : list or tuple
list of tuple of point or linestring coordinates
exterior : list or tuple
2d list of polygon coordinates
interiors : list or tuple
2d or 3d list of polygon interiors
"""
def __init__(
self,
shapetype,
coordinates=None,
exterior=None,
interiors=None,
):
self.__type = shapetype
if shapetype == "Polygon":
self.exterior = tuple(map(tuple, exterior))
self.interiors = (
() if interiors is None else (tuple(map(tuple, i)) for i in interiors)
)
self.interiors = tuple(self.interiors)
elif shapetype == "LineString":
self.coords = list(map(tuple, coordinates))
if len(self.coords[0]) == 3:
self.has_z = True
elif shapetype == "Point":
while len(coordinates) == 1:
coordinates = coordinates[0]
self.coords = coordinates
if len(coordinates) == 3:
self.has_z = True
else:
err = (
"Supported shape types are Polygon, LineString, "
f"and Point: Supplied shape type {shapetype}"
)
raise TypeError(err)
@property
def __geo_interface__(self):
"""
Creates the geojson standard representation of a shape
Returns
-------
dict
"""
geo_interface = {}
if self.__type == "Polygon":
geo_interface = {
"coordinates": tuple([self.exterior] + list(self.interiors)),
"type": self.__type,
}
elif self.__type == "LineString":
geo_interface = {"coordinates": tuple(self.coords), "type": self.__type}
elif self.__type == "Point":
geo_interface = {"coordinates": tuple(self.coords), "type": self.__type}
return geo_interface
@property
def geojson(self):
return self.__geo_interface__
[docs] @staticmethod
def from_geojson(geo_interface):
"""
Method to load from geojson
Parameters
----------
geo_interface : geojson, dict
geojson compliant representation of a linestring
Returns
-------
Polygon, LineString, or Point
"""
if geo_interface["type"] in {"Polygon", "MultiPolygon"}:
coord_list = geo_interface["coordinates"]
if geo_interface["type"] == "Polygon":
coord_list = [coord_list]
geoms = []
for coords in coord_list:
exteriors = coords[0]
interiors = None
if len(coords) > 1:
interiors = coords[1:]
geoms.append(Polygon(exteriors, interiors))
if len(geoms) == 1:
shape = geoms[0]
else:
shape = MultiPolygon(geoms)
elif geo_interface["type"] == "LineString":
shape = LineString(geo_interface["coordinates"])
elif geo_interface["type"] == "MultiLineString":
geoms = [LineString(coords) for coords in geo_interface["coordinates"]]
shape = MultiLineString(geoms)
elif geo_interface["type"] == "Point":
shape = Point(geo_interface["coordinates"])
elif geo_interface["type"] == "MultiPoint":
geoms = [Point(coords) for coords in geo_interface["coordinates"]]
shape = MultiPoint(geoms)
else:
err = (
"Supported shape types are Polygon, LineString, and "
"Point: Supplied shape type {}".format(geo_interface["type"])
)
raise TypeError(err)
return shape
[docs]class Collection(list):
"""
The collection object is container for a group of flopy geometries
This class acts as a base class for MultiPoint, MultiLineString, and
MultiPolygon classes. This class can also accept a mix of geometries
and act as a stand alone container.
Parameters
----------
geometries : list
list of flopy.util.geometry objects
"""
def __init__(self, geometries=()):
super().__init__(geometries)
def __repr__(self):
return f"Shapes: {list(self)}"
@property
def __geo_interface__(self):
return {
"type": "GeometryCollection",
"geometries": [g.__geo_interface__ for g in self],
}
@property
def bounds(self):
"""
Method to calculate the bounding box of the collection
Returns
-------
tuple (xmin, ymin, xmax, ymax)
"""
bbox = [geom.bounds for geom in self]
xmin, ymin = np.min(bbox, axis=0)[0:2]
xmax, ymax = np.max(bbox, axis=0)[2:]
return xmin, ymin, xmax, ymax
def __reversed__(self):
for shp in self:
reversed(shp)
return self
[docs] def plot(self, ax=None, **kwargs):
"""
Plotting method for collection
Parameters
----------
ax : matplotlib.axes object
kwargs : keyword arguments
matplotlib keyword arguments
Returns
-------
matplotlib.axes object
"""
for g in self:
ax = g.plot(ax=ax, **kwargs)
xmin, ymin, xmax, ymax = self.bounds
ax.set_ylim([ymin - 0.005, ymax + 0.005])
ax.set_xlim([xmin - 0.005, xmax + 0.005])
return ax
[docs]class MultiPolygon(Collection):
"""
Container for housing and describing multipolygon geometries (e.g. to be
read or written to shapefiles or other geographic data formats)
Parameters
----------
polygons : list, tuple, default ()
list of flopy.utils.geometry.Polygon objects
"""
def __init__(self, polygons=()):
for p in polygons:
if not isinstance(p, Polygon):
raise TypeError("Only Polygon instances are supported")
super().__init__(polygons)
def __repr__(self):
return f"MultiPolygon: {list(self)}"
@property
def __geo_interface__(self):
return {
"type": "MultiPolygon",
"coordinates": [g.__geo_interface__["coordinates"] for g in self],
}
[docs]class MultiLineString(Collection):
"""
Container for housing and describing multilinestring geometries (e.g. to be
read or written to shapefiles or other geographic data formats)
Parameters
----------
linestrings : list, tuple, default ()
list of flopy.utils.geometry.LineString objects
"""
def __init__(self, linestrings=()):
for l in linestrings:
if not isinstance(l, LineString):
raise TypeError("Only LineString instances are supported")
super().__init__(linestrings)
def __repr__(self):
return f"LineString: {list(self)}"
@property
def __geo_interface__(self):
return {
"type": "MultiLineString",
"coordinates": [g.__geo_interface__["coordinates"] for g in self],
}
[docs]class MultiPoint(Collection):
"""
Container for housing and describing multipoint geometries (e.g. to be
read or written to shapefiles or other geographic data formats)
Parameters
----------
points : list, tuple, default ()
list of flopy.utils.geometry.Point objects
"""
def __init__(self, points=()):
for p in points:
if not isinstance(p, Point):
raise TypeError("Only Point instances are supported")
super().__init__(points)
def __repr__(self):
return f"MultiPoint: {list(self)}"
@property
def __geo_interface__(self):
return {
"type": "MultiPoint",
"coordinates": [g.__geo_interface__["coordinates"] for g in self],
}
[docs]class Polygon(Shape):
type = "Polygon"
shapeType = 5 # pyshp
def __init__(self, exterior, interiors=None):
"""
Container for housing and describing polygon geometries (e.g. to be
read or written to shapefiles or other geographic data formats)
Parameters
----------
exterior : sequence
Sequence of coordinates describing the outer ring of the polygon.
interiors : sequence of sequences
Describes one or more holes within the polygon
Attributes
----------
exterior : (x, y, z) coordinates of exterior
interiors : tuple of (x, y, z) coordinates of each interior polygon
patch : descartes.PolygonPatch representation
bounds : (xmin, ymin, xmax, ymax)
Tuple describing bounding box for polygon
geojson : dict
Returns a geojson representation of the feature
pyshp_parts : list of lists
Returns a list of all parts (each an individual polygon).
Can be used as input for the shapefile.Writer.poly method
(pyshp package)
Methods
-------
get_patch
Returns a descartes PolygonPatch object representation of the
polygon. Accepts keyword arguments to descartes.PolygonPatch.
Requires the descartes package (pip install descartes).
plot
Plots the feature using descartes (via get_patch) and
matplotlib.pyplot. Accepts keyword arguments to
descartes.PolygonPatch. Requires the descartes package
(pip install descartes).
Notes
-----
Multi-polygons not yet supported.
z information is only stored if it was entered.
"""
super().__init__(
self.type, coordinates=None, exterior=exterior, interiors=interiors
)
def __eq__(self, other):
if not isinstance(other, Polygon):
return False
if other.exterior != self.exterior:
return False
if other.interiors != self.interiors:
return False
return True
@property
def _exterior_x(self):
return [x for x, y in self.exterior]
@property
def _exterior_y(self):
return [y for x, y in self.exterior]
@property
def bounds(self):
ymin = np.min(self._exterior_y)
ymax = np.max(self._exterior_y)
xmin = np.min(self._exterior_x)
xmax = np.max(self._exterior_x)
return xmin, ymin, xmax, ymax
@property
def pyshp_parts(self):
# exterior ring must be clockwise (negative area)
# interiors rings must be counter-clockwise (positive area)
shapefile = import_optional_dependency("shapefile")
exterior = list(self.exterior)
if shapefile.signed_area(exterior) > 0:
exterior.reverse()
interiors = []
for i in self.interiors:
il = list(i)
if shapefile.signed_area(il) < 0:
il.reverse()
interiors.append(il)
result = [exterior]
for i in interiors:
result.append(i)
return result
@property
def patch(self):
return self.get_patch()
def __reversed__(self):
# method to reverse the sorting on polygon points, patch for
# differences in pyshp 2.2.0 vs pyshp < 2.2.0
if self.exterior:
self.exterior = self.exterior[::-1]
if self.interiors:
interiors = []
for i in self.interiors:
interiors.append(i[::-1])
self.interiors = interiors
return self
[docs] def get_patch(self, **kwargs):
descartes = import_optional_dependency("descartes")
from descartes import PolygonPatch
return PolygonPatch(self.geojson, **kwargs)
[docs] def plot(self, ax=None, **kwargs):
"""
Plot the feature.
Parameters
----------
ax : matplotlib.pyplot axes instance
Accepts keyword arguments to descartes.PolygonPatch. Requires the
descartes package (pip install descartes).
"""
import matplotlib.pyplot as plt
if ax is None:
ax = plt.gca()
try:
ax.add_patch(self.get_patch(**kwargs))
xmin, ymin, xmax, ymax = self.bounds
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
except:
print("could not plot polygon feature")
return ax
[docs]class LineString(Shape):
type = "LineString"
shapeType = 3
has_z = False
def __init__(self, coordinates):
"""
Container for housing and describing linestring geometries (e.g. to be
read or written to shapefiles or other geographic data formats)
Parameters
----------
coordinates : sequence
Sequence of coordinates describing a line
Attributes
----------
coords : list of (x, y, z) coordinates
x : list of x coordinates
y : list of y coordinates
z : list of z coordinates
bounds : (xmin, ymin, xmax, ymax)
Tuple describing bounding box for linestring
geojson : dict
Returns a geojson representation of the feature
pyshp_parts : list of lists
Returns a list of all parts (each an individual linestring).
Can be used as input for the shapefile.Writer.line method (pyshp package)
Methods
-------
plot
Plots the feature using matplotlib.pyplot.
Accepts keyword arguments to pyplot.plot.
Notes
-----
Multi-linestrings not yet supported.
z information is only stored if it was entered.
"""
super().__init__(self.type, coordinates)
def __eq__(self, other):
if not isinstance(other, LineString):
return False
if other.x != self.x:
return False
if other.y != self.y:
return False
if other.z != self.z:
return False
return True
@property
def x(self):
return [c[0] for c in self.coords]
@property
def y(self):
return [c[1] for c in self.coords]
@property
def z(self):
return 0 if not self.has_z else [c[2] for c in self.coords]
@property
def bounds(self):
ymin = np.min(self.y)
ymax = np.max(self.y)
xmin = np.min(self.x)
xmax = np.max(self.x)
return xmin, ymin, xmax, ymax
@property
def pyshp_parts(self):
return [self.coords]
def __reversed__(self):
self.coords = self.coords[::-1]
return self
[docs] def plot(self, ax=None, **kwargs):
import matplotlib.pyplot as plt
if ax is None:
ax = plt.gca()
ax.plot(self.x, self.y, **kwargs)
xmin, ymin, xmax, ymax = self.bounds
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
return ax
[docs]class Point(Shape):
type = "Point"
shapeType = 1
has_z = False
def __init__(self, *coordinates):
"""
Container for housing and describing point geometries (e.g. to be read
or written to shapefiles or other geographic data formats)
Parameters
----------
coordinates : tuple
x, y or x, y, z
Attributes
----------
coords : x, y, z coordinates
x : x coordinate
y : y coordinate
z : z coordinate
bounds : (xmin, ymin, xmax, ymax)
Tuple describing bounding box
geojson : dict
Returns a geojson representation of the feature
pyshp_parts : list of tuples
Can be used as input for the shapefile.Writer.line method
(pyshp package)
Methods
-------
plot
Plots the feature using matplotlib.pyplot.
Accepts keyword arguments to pyplot.scatter.
Notes
-----
z information is only stored if it was entered.
"""
super().__init__(self.type, coordinates)
def __eq__(self, other):
if not isinstance(other, Point):
return False
if other.x != self.x:
return False
if other.y != self.y:
return False
if other.z != self.z:
return False
return True
@property
def x(self):
return self.coords[0]
@property
def y(self):
return self.coords[1]
@property
def z(self):
return 0 if not self.has_z else self.coords[2]
@property
def bounds(self):
ymin = np.min(self.y)
ymax = np.max(self.y)
xmin = np.min(self.x)
xmax = np.max(self.x)
return xmin, ymin, xmax, ymax
@property
def pyshp_parts(self):
return self.coords
def __reversed__(self):
return self
[docs] def plot(self, ax=None, **kwargs):
import matplotlib.pyplot as plt
if ax is None:
ax = plt.gca()
ax.scatter(self.x, self.y, **kwargs)
xmin, ymin, xmax, ymax = self.bounds
ax.set_xlim(xmin - 1, xmax + 1) # singular bounds otherwise
ax.set_ylim(ymin - 1, ymax + 1)
return ax
[docs]def rotate(x, y, xoff, yoff, angrot_radians):
"""
Given x and y array-like values calculate the rotation about an
arbitrary origin and then return the rotated coordinates.
"""
if isinstance(x, list):
x = np.array(x)
if isinstance(y, list):
y = np.array(y)
xrot = (
xoff + np.cos(angrot_radians) * (x - xoff) - np.sin(angrot_radians) * (y - yoff)
)
yrot = (
yoff + np.sin(angrot_radians) * (x - xoff) + np.cos(angrot_radians) * (y - yoff)
)
return xrot, yrot
[docs]def get_polygon_area(geom):
"""
Calculate the area of a closed polygon
Parameters
----------
geom : geospatial representation of polygon
accepted types:
vertices np.array([(x, y),....])
geojson.Polygon
shapely.Polygon
shapefile.Shape
Returns
-------
area : float
area of polygon centroid
"""
from .geospatial_utils import GeoSpatialUtil
if isinstance(geom, (list, tuple, np.ndarray)):
geom = [geom]
geom = GeoSpatialUtil(geom, shapetype="Polygon")
verts = np.array(geom.points[0])
nverts = verts.shape[0]
a = 0.0
for iv in range(nverts - 1):
x = verts[iv, 0]
y = verts[iv, 1]
xp1 = verts[iv + 1, 0]
yp1 = verts[iv + 1, 1]
a += x * yp1 - xp1 * y
a = abs(a * 0.5)
return a
[docs]def get_polygon_centroid(geom):
"""
Calculate the centroid of a closed polygon
Parameters
----------
geom : geospatial representation of polygon
accepted types:
vertices np.array([(x, y),....])
geojson.Polygon
shapely.Polygon
shapefile.Shape
Returns
-------
centroid : tuple
(x, y) of polygon centroid
"""
from .geospatial_utils import GeoSpatialUtil
if isinstance(geom, (list, tuple, np.ndarray)):
geom = [geom]
geom = GeoSpatialUtil(geom, shapetype="Polygon")
verts = np.array(geom.points[0])
nverts = verts.shape[0]
cx = 0.0
cy = 0.0
for i in range(nverts - 1):
x = verts[i, 0]
y = verts[i, 1]
xp1 = verts[i + 1, 0]
yp1 = verts[i + 1, 1]
cx += (x + xp1) * (x * yp1 - xp1 * y)
cy += (y + yp1) * (x * yp1 - xp1 * y)
a = get_polygon_area(verts)
cx = cx * 1.0 / 6.0 / a
cy = cy * 1.0 / 6.0 / a
return cx, cy
[docs]def is_clockwise(*geom):
"""
Determine if a ring is defined clockwise
Parameters
----------
*geom : geospatial representation of polygon
accepted types:
vertices [(x, y),....]
geojson.Polygon
shapely.Polygon
shapefile.Shape
x and y vertices: [x1, x2, x3], [y1, y2, y3]
Returns
-------
clockwise : bool
True when the ring is defined clockwise, False otherwise
"""
from .geospatial_utils import GeoSpatialUtil
if len(geom) == 2:
x, y = geom
else:
geom = GeoSpatialUtil(geom, shapetype="Polygon")
x, y = np.array(geom.points[0]).T
if not ((x[0] == x[-1]) and (y[0] == y[-1])):
# close the ring if needed
x = np.append(x, x[0])
y = np.append(y, y[0])
return np.sum((np.diff(x)) * (y[1:] + y[:-1])) > 0
[docs]def point_in_polygon(xc, yc, polygon):
"""
Use the ray casting algorithm to determine if a point
is within a polygon. Enables very fast
intersection calculations!
Parameters
----------
xc : np.ndarray
2d array of xpoints
yc : np.ndarray
2d array of ypoints
polygon : iterable (list)
polygon vertices [(x0, y0),....(xn, yn)]
note: polygon can be open or closed
Returns
-------
mask: np.array
True value means point is in polygon!
"""
x0, y0 = polygon[0]
xt, yt = polygon[-1]
# close polygon if it isn't already
if (x0, y0) != (xt, yt):
polygon.append((x0, y0))
ray_count = np.zeros(xc.shape, dtype=int)
num = len(polygon)
j = num - 1
for i in range(num):
tmp = polygon[i][0] + (polygon[j][0] - polygon[i][0]) * (yc - polygon[i][1]) / (
polygon[j][1] - polygon[i][1]
)
comp = np.asarray(
((polygon[i][1] > yc) ^ (polygon[j][1] > yc)) & (xc < tmp)
).nonzero()
j = i
if len(comp[0]) > 0:
ray_count[comp[0], comp[1]] += 1
mask = np.ones(xc.shape, dtype=bool)
mask[ray_count % 2 == 0] = False
return mask
[docs]def project_point_onto_xc_line(line, pts, d0=0, calc_dist=False):
"""
Method to project points onto a cross sectional line
that is defined by distance. Used for plotting MODPATH results
on to a cross section!
line : list or np.ndarray
numpy array of [(x0, y0), (x1, y1)] that defines the line
to project on to
pts : list or np.ndarray
numpy array of [(x, y),] points to be projected
d0 : distance offset along line of min(xl)
calc_dist : bool
boolean flag to indicate that the return type is distance offset
by d0 (calculated from line[0]). If false this will return the
"projected" xy location along the cross-sectional line
Returns:
tuple of (x , y) or distance
"""
if isinstance(line, list):
line = np.array(line)
if isinstance(pts, list):
pts = np.array(pts)
if pts.ndim == 1:
pts = np.expand_dims(pts, axis=0)
x0, x1 = line.T[0, :]
y0, y1 = line.T[1, :]
dx = x1 - x0
dy = y1 - y0
if dx == 0:
dx = 1e-10 # trap the vertical line condition to avoid inf slope
m = dy / dx
b = y0 - (m * x0)
x = pts.T[0]
y = pts.T[1]
bx = ((m * y) + (x - m * b)) / (1 + m**2)
by = ((m**2 * y) + (m * x) + b) / (1 + m**2)
if calc_dist:
# get distance between bxy, xy0
dist = distance(bx, by, x0, y0)
# get distance between xy0, xy1
dist0 = distance(x0, y0, x1, y1)
# get distance between bxy, xy1
dist1 = distance(bx, by, x1, y1)
adj = np.full((dist.size,), 1)
for ix, d in enumerate(dist):
if d <= dist0:
if dist1[ix] > dist0:
adj[ix] = -1
else:
if dist1[ix] > d:
adj[ix] = -1
dist *= adj
dist += d0
if len(dist) == 1:
dist = dist[0]
return dist
else:
return bx, by
[docs]def distance(x0, y0, x1, y1):
"""
General distance equation
Parameters
----------
x0 : float
y0 : float
x1 : np.array or float
y1 : np.array or float
Returns
-------
distance
"""
asq = (x0 - x1) ** 2
bsq = (y0 - y1) ** 2
dist = np.sqrt(asq + bsq)
return dist