Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Quick fixes for gridintersect using linestrings and structured grid #996

Merged
merged 6 commits into from
Sep 25, 2020
Merged
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
52 changes: 30 additions & 22 deletions flopy/utils/gridintersect.py
Original file line number Diff line number Diff line change
Expand Up @@ -710,7 +710,7 @@ def _intersect_linestring_structured(self, shp, keepzerolengths=False):
nodelist, lengths, vertices = [], [], []
ixshapes = []
for ls in lineclip:
n, l, v, ix = self._get_nodes_intersecting_linestring(ls)
n, l, v, ixs = self._get_nodes_intersecting_linestring(ls)
nodelist += n
lengths += l
# if necessary, transform coordinates back to real
Expand All @@ -722,26 +722,31 @@ def _intersect_linestring_structured(self, shp, keepzerolengths=False):
):
v_realworld = []
for pt in v:
pt = np.array(pt)
rx, ry = transform(
[pt[0]],
[pt[1]],
pt[:, 0],
pt[:, 1],
self.mfgrid.xoffset,
self.mfgrid.yoffset,
self.mfgrid.angrot_radians,
inverse=False,
)
v_realworld.append([rx, ry])
ix_realworld = rotate(
ix, self.mfgrid.angrot, origin=(0.0, 0.0)
)
ix_realworld = translate(
ix_realworld, self.mfgrid.xoffset, self.mfgrid.yoffset
)
ixs_realworld = []
for ix in ixs:
ix_realworld = rotate(
ix, self.mfgrid.angrot, origin=(0.0, 0.0)
)
ix_realworld = translate(
ix_realworld,
self.mfgrid.xoffset, self.mfgrid.yoffset
)
ixs_realworld.append(ix_realworld)
else:
v_realworld = v
ix_realworld = ix
ixs_realworld = ixs
vertices += v_realworld
ixshapes += ix_realworld
ixshapes += ixs_realworld
else: # linestring is fully within grid
(
nodelist,
Expand All @@ -758,9 +763,10 @@ def _intersect_linestring_structured(self, shp, keepzerolengths=False):
):
v_realworld = []
for pt in vertices:
pt = np.array(pt)
rx, ry = transform(
[pt[0]],
[pt[1]],
pt[:, 0],
pt[:, 1],
self.mfgrid.xoffset,
self.mfgrid.yoffset,
self.mfgrid.angrot_radians,
Expand Down Expand Up @@ -887,7 +893,7 @@ def _get_nodes_intersecting_linestring(self, linestring):
ixshapes.append(intersect)
length = intersect.length
lengths.append(length)
if intersect.geom_type == "MultiLineString":
if hasattr(intersect, "geoms"):
x, y = [], []
for igeom in intersect.geoms:
x.append(igeom.xy[0])
Expand Down Expand Up @@ -972,7 +978,7 @@ def _check_adjacent_cells_intersecting_line(
intersect = linestring.intersection(pl)
ixshape.append(intersect)
length.append(intersect.length)
if intersect.geom_type == "MultiLineString":
if hasattr(intersect, "geoms"):
x, y = [], []
for igeom in intersect.geoms:
x.append(igeom.xy[0])
Expand All @@ -999,7 +1005,7 @@ def _check_adjacent_cells_intersecting_line(
intersect = linestring.intersection(pl)
ixshape.append(intersect)
length.append(intersect.length)
if intersect.geom_type == "MultiLineString":
if hasattr(intersect, "geoms"):
x, y = [], []
for igeom in intersect.geoms:
x.append(igeom.xy[0])
Expand All @@ -1026,7 +1032,7 @@ def _check_adjacent_cells_intersecting_line(
intersect = linestring.intersection(pl)
ixshape.append(intersect)
length.append(intersect.length)
if intersect.geom_type == "MultiLineString":
if hasattr(intersect, "geoms"):
x, y = [], []
for igeom in intersect.geoms:
x.append(igeom.xy[0])
Expand All @@ -1053,7 +1059,7 @@ def _check_adjacent_cells_intersecting_line(
intersect = linestring.intersection(pl)
ixshape.append(intersect)
length.append(intersect.length)
if intersect.geom_type == "MultiLineString":
if hasattr(intersect, "geoms"):
x, y = [], []
for igeom in intersect.geoms:
x.append(igeom.xy[0])
Expand Down Expand Up @@ -1302,9 +1308,10 @@ def _transform_geo_interface_polygon(self, polygon):
# transform shell coordinates
shell_pts = []
for pt in shell:
pt = np.array(pt)
rx, ry = transform(
[pt[0]],
[pt[1]],
pt[:, 0],
pt[:, 1],
self.mfgrid.xoffset,
self.mfgrid.yoffset,
self.mfgrid.angrot_radians,
Expand All @@ -1316,9 +1323,10 @@ def _transform_geo_interface_polygon(self, polygon):
if holes:
holes_pts = []
for pt in holes:
pt = np.array(pt)
rx, ry = transform(
[pt[0]],
[pt[1]],
pt[:, 0],
pt[:, 1],
self.mfgrid.xoffset,
self.mfgrid.yoffset,
self.mfgrid.angrot_radians,
Expand Down