Gael-de-Sailly
/
flopy
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

feat(geospatial_util): geospatial consolidation utilities (#1002) 

* *Added GeoSpatialUtil for geospatial conversions
*Added GeoSpatialCollection for batch conversion of "Collection" objects

*Updated geometry.py and added additional geometry types to geometry.py

*updates to gridgen, shapefile_utils, triangle, rasters, and gridintersect modules that allow the user to supply any supported shape object to these geospatial methods. Supported shape objects are, geojson, shapely, shapefile, flopy.utils.geometry, and lists of vertices.

*Added t074_test_geospatial_util.py

issue #772

* remove intersect_linestring, intersect_point, and intersect_polygon from GridIntersect.

*Refactored t065_test_gridintersect.py for intersect() method
*Updated flopy3_grid_intersection_demo.ipynb for intersect() method
develop
Joshua Larsen 2020-10-19 10:40:49 -07:00 committed by GitHub
parent 698b4ecfc2
commit 9ad6732c4b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
10 changed files with 1677 additions and 320 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

136
autotest/t065_test_gridintersect.py
View File

@ -128,7 +128,7 @@ def test_rect_grid_3d_point_outside():
botm = np.concatenate([np.ones(4), np.zeros(4)]).reshape(2, 2, 2)
gr = get_rect_grid(top=np.ones(4), botm=botm)
ix = GridIntersect(gr, method="structured")
result = ix.intersect_point(Point(25., 25., .5))
result = ix.intersect(Point(25., 25., .5))
assert len(result) == 0
return result
@ -142,7 +142,7 @@ def test_rect_grid_3d_point_inside():
botm = np.concatenate([np.ones(4), .5 * np.ones(4), np.zeros(4)]).reshape(3, 2, 2)
gr = get_rect_grid(top=np.ones(4), botm=botm)
ix = GridIntersect(gr, method="structured")
result = ix.intersect_point(Point(2., 2., .2))
result = ix.intersect(Point(2., 2., .2))
assert result.cellids[0] == (1, 1, 0)
return result
@ -156,7 +156,7 @@ def test_rect_grid_3d_point_above():
botm = np.concatenate([np.ones(4), np.zeros(4)]).reshape(2, 2, 2)
gr = get_rect_grid(top=np.ones(4), botm=botm)
ix = GridIntersect(gr, method="structured")
result = ix.intersect_point(Point(2., 2., 2))
result = ix.intersect(Point(2., 2., 2))
assert len(result) == 0
return result
@ -169,7 +169,7 @@ def test_rect_grid_point_outside():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_point(Point(25., 25.))
result = ix.intersect(Point(25., 25.))
assert len(result) == 0
return result
@ -182,7 +182,7 @@ def test_rect_grid_point_on_outer_boundary():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_point(Point(20., 10.))
result = ix.intersect(Point(20., 10.))
assert len(result) == 1
assert np.all(result.cellids[0] == (0, 1))
return result
@ -196,7 +196,7 @@ def test_rect_grid_point_on_inner_boundary():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_point(Point(10., 10.))
result = ix.intersect(Point(10., 10.))
assert len(result) == 1
assert np.all(result.cellids[0] == (0, 0))
return result
@ -210,7 +210,7 @@ def test_rect_grid_multipoint_in_one_cell():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_point(MultiPoint([Point(1., 1.), Point(2., 2.)]))
result = ix.intersect(MultiPoint([Point(1., 1.), Point(2., 2.)]))
assert len(result) == 1
assert result.cellids[0] == (1, 0)
return result
@ -224,7 +224,7 @@ def test_rect_grid_multipoint_in_multiple_cells():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_point(MultiPoint([Point(1., 1.), Point(12., 12.)]))
result = ix.intersect(MultiPoint([Point(1., 1.), Point(12., 12.)]))
assert len(result) == 2
assert result.cellids[0] == (1, 0)
assert result.cellids[1] == (0, 1)
@ -242,7 +242,7 @@ def test_rect_grid_point_outside_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_point(Point(25., 25.))
result = ix.intersect(Point(25., 25.))
assert len(result) == 0
return result
@ -255,7 +255,7 @@ def test_rect_grid_point_on_outer_boundary_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_point(Point(20., 10.))
result = ix.intersect(Point(20., 10.))
assert len(result) == 1
assert np.all(result.cellids[0] == (0, 1))
return result
@ -269,7 +269,7 @@ def test_rect_grid_point_on_inner_boundary_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_point(Point(10., 10.))
result = ix.intersect(Point(10., 10.))
assert len(result) == 1
assert np.all(result.cellids[0] == (0, 0))
return result
@ -283,16 +283,16 @@ def test_rect_vertex_grid_point_in_one_cell_shapely(rtree=True):
return
gr = get_rect_vertex_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_point(Point(4., 4.))
result = ix.intersect(Point(4., 4.))
assert len(result) == 1
assert result.cellids[0] == 0
result = ix.intersect_point(Point(4., 6.))
result = ix.intersect(Point(4., 6.))
assert len(result) == 1
assert result.cellids[0] == 0
result = ix.intersect_point(Point(6., 6.))
result = ix.intersect(Point(6., 6.))
assert len(result) == 1
assert result.cellids[0] == 0
result = ix.intersect_point(Point(6., 4.))
result = ix.intersect(Point(6., 4.))
assert len(result) == 1
assert result.cellids[0] == 0
return result
@ -306,7 +306,7 @@ def test_rect_grid_multipoint_in_one_cell_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_point(MultiPoint([Point(1., 1.), Point(2., 2.)]))
result = ix.intersect(MultiPoint([Point(1., 1.), Point(2., 2.)]))
assert len(result) == 1
assert result.cellids[0] == (1, 0)
return result
@ -320,7 +320,7 @@ def test_rect_grid_multipoint_in_multiple_cells_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_point(MultiPoint([Point(1., 1.), Point(12., 12.)]))
result = ix.intersect(MultiPoint([Point(1., 1.), Point(12., 12.)]))
assert len(result) == 2
assert result.cellids[0] == (0, 1)
assert result.cellids[1] == (1, 0)
@ -337,7 +337,7 @@ def test_tri_grid_point_outside(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_point(Point(25., 25.))
result = ix.intersect(Point(25., 25.))
assert len(result) == 0
return result
@ -352,7 +352,7 @@ def test_tri_grid_point_on_outer_boundary(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_point(Point(20., 10.))
result = ix.intersect(Point(20., 10.))
assert len(result) == 1
assert np.all(result.cellids[0] == 0)
return result
@ -368,7 +368,7 @@ def test_tri_grid_point_on_inner_boundary(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_point(Point(10., 10.))
result = ix.intersect(Point(10., 10.))
assert len(result) == 1
assert np.all(result.cellids[0] == 0)
return result
@ -384,7 +384,7 @@ def test_tri_grid_multipoint_in_one_cell(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_point(MultiPoint([Point(1., 1.), Point(2., 2.)]))
result = ix.intersect(MultiPoint([Point(1., 1.), Point(2., 2.)]))
assert len(result) == 1
assert result.cellids[0] == 1
return result
@ -400,7 +400,7 @@ def test_tri_grid_multipoint_in_multiple_cells(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_point(MultiPoint([Point(1., 1.), Point(12., 12.)]))
result = ix.intersect(MultiPoint([Point(1., 1.), Point(12., 12.)]))
assert len(result) == 2
assert result.cellids[0] == 0
assert result.cellids[1] == 1
@ -418,7 +418,7 @@ def test_rect_grid_linestring_outside():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_linestring(LineString([(25., 25.), (21., 5.)]))
result = ix.intersect(LineString([(25., 25.), (21., 5.)]))
assert len(result) == 0
return result
@ -431,7 +431,7 @@ def test_rect_grid_linestring_in_2cells():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_linestring(LineString([(5., 5.), (15., 5.)]))
result = ix.intersect(LineString([(5., 5.), (15., 5.)]))
assert len(result) == 2
assert result.lengths.sum() == 10.
assert result.cellids[0] == (1, 0)
@ -447,7 +447,7 @@ def test_rect_grid_linestring_on_outer_boundary():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_linestring(LineString([(15., 20.), (5., 20.)]))
result = ix.intersect(LineString([(15., 20.), (5., 20.)]))
assert len(result) == 2
assert result.lengths.sum() == 10.
assert result.cellids[1] == (0, 0)
@ -463,7 +463,7 @@ def test_rect_grid_linestring_on_inner_boundary():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_linestring(LineString([(5., 10.), (15., 10.)]))
result = ix.intersect(LineString([(5., 10.), (15., 10.)]))
assert len(result) == 2
assert result.lengths.sum() == 10.
assert result.cellids[0] == (0, 0)
@ -479,7 +479,7 @@ def test_rect_grid_multilinestring_in_one_cell():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_linestring(MultiLineString(
result = ix.intersect(MultiLineString(
[LineString([(1., 1), (9., 1.)]), LineString([(1., 9.), (9., 9.)])]))
assert len(result) == 1
assert result.lengths == 16.
@ -495,7 +495,7 @@ def test_rect_grid_linestring_in_and_out_of_cell():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_linestring(
result = ix.intersect(
LineString([(5., 9), (15., 5.), (5., 1.)]))
assert len(result) == 2
assert result.cellids[0] == (1, 0)
@ -512,7 +512,7 @@ def test_rect_grid_linestring_in_and_out_of_cell2():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_linestring(LineString(
result = ix.intersect(LineString(
[(5, 15), (5., 9), (15., 5.), (5., 1.)]))
assert len(result) == 3
# assert result.cellids[0] == (1, 0)
@ -532,7 +532,7 @@ def test_rect_grid_linestring_outside_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_linestring(LineString([(25., 25.), (21., 5.)]))
result = ix.intersect(LineString([(25., 25.), (21., 5.)]))
assert len(result) == 0
return result
@ -545,7 +545,7 @@ def test_rect_grid_linestring_in_2cells_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_linestring(LineString([(5., 5.), (15., 5.)]))
result = ix.intersect(LineString([(5., 5.), (15., 5.)]))
assert len(result) == 2
assert result.lengths.sum() == 10.
assert result.cellids[0] == (1, 0)
@ -561,7 +561,7 @@ def test_rect_grid_linestring_on_outer_boundary_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_linestring(LineString([(15., 20.), (5., 20.)]))
result = ix.intersect(LineString([(15., 20.), (5., 20.)]))
assert len(result) == 2
assert result.lengths.sum() == 10.
assert result.cellids[0] == (0, 0)
@ -577,7 +577,7 @@ def test_rect_grid_linestring_on_inner_boundary_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_linestring(LineString([(5., 10.), (15., 10.)]))
result = ix.intersect(LineString([(5., 10.), (15., 10.)]))
assert len(result) == 2
assert result.lengths.sum() == 10.
assert result.cellids[0] == (0, 0)
@ -593,7 +593,7 @@ def test_rect_grid_multilinestring_in_one_cell_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_linestring(MultiLineString(
result = ix.intersect(MultiLineString(
[LineString([(1., 1), (9., 1.)]), LineString([(1., 9.), (9., 9.)])]))
assert len(result) == 1
assert result.lengths == 16.
@ -609,7 +609,7 @@ def test_rect_grid_linestring_in_and_out_of_cell_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_linestring(
result = ix.intersect(
LineString([(5., 9), (15., 5.), (5., 1.)]))
assert len(result) == 2
assert result.cellids[0] == (1, 0)
@ -628,7 +628,7 @@ def test_tri_grid_linestring_outside(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_linestring(LineString([(25., 25.), (21., 5.)]))
result = ix.intersect(LineString([(25., 25.), (21., 5.)]))
assert len(result) == 0
return result
@ -643,7 +643,7 @@ def test_tri_grid_linestring_in_2cells(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_linestring(LineString([(5., 5.), (5., 15.)]))
result = ix.intersect(LineString([(5., 5.), (5., 15.)]))
assert len(result) == 2
assert result.lengths.sum() == 10.
assert result.cellids[0] == 1
@ -661,7 +661,7 @@ def test_tri_grid_linestring_on_outer_boundary(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_linestring(LineString([(15., 20.), (5., 20.)]))
result = ix.intersect(LineString([(15., 20.), (5., 20.)]))
assert len(result) == 2
assert result.lengths.sum() == 10.
assert result.cellids[0] == 2
@ -679,7 +679,7 @@ def test_tri_grid_linestring_on_inner_boundary(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_linestring(LineString([(5., 10.), (15., 10.)]))
result = ix.intersect(LineString([(5., 10.), (15., 10.)]))
assert len(result) == 2
assert result.lengths.sum() == 10.
assert result.cellids[0] == 0
@ -697,7 +697,7 @@ def test_tri_grid_multilinestring_in_one_cell(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_linestring(MultiLineString(
result = ix.intersect(MultiLineString(
[LineString([(1., 1), (9., 1.)]), LineString([(2., 2.), (9., 2.)])]))
assert len(result) == 1
assert result.lengths == 15.
@ -716,7 +716,7 @@ def test_rect_grid_polygon_outside():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_polygon(
result = ix.intersect(
Polygon([(21., 11.), (23., 17.), (25., 11.)]))
assert len(result) == 0
return result
@ -730,7 +730,7 @@ def test_rect_grid_polygon_in_2cells():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_polygon(
result = ix.intersect(
Polygon([(2.5, 5.0), (7.5, 5.0), (7.5, 15.), (2.5, 15.)]))
assert len(result) == 2
assert result.areas.sum() == 50.
@ -745,7 +745,7 @@ def test_rect_grid_polygon_on_outer_boundary():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_polygon(
result = ix.intersect(
Polygon([(20., 5.0), (25., 5.0), (25., 15.), (20., 15.)]))
assert len(result) == 0
return result
@ -759,7 +759,7 @@ def test_rect_grid_polygon_on_inner_boundary():
return
gr = get_rect_grid()
ix = GridIntersect(gr, method="structured")
result = ix.intersect_polygon(
result = ix.intersect(
Polygon([(5., 10.0), (15., 10.0), (15., 5.), (5., 5.)]))
assert len(result) == 2
assert result.areas.sum() == 50.
@ -777,7 +777,7 @@ def test_rect_grid_multipolygon_in_one_cell():
p1 = Polygon([(1., 1.), (8., 1.), (8., 3.), (1., 3.)])
p2 = Polygon([(1., 9.), (8., 9.), (8., 7.), (1., 7.)])
p = MultiPolygon([p1, p2])
result = ix.intersect_polygon(p)
result = ix.intersect(p)
assert len(result) == 1
assert result.areas.sum() == 28.
return result
@ -794,7 +794,7 @@ def test_rect_grid_multipolygon_in_multiple_cells():
p1 = Polygon([(1., 1.), (19., 1.), (19., 3.), (1., 3.)])
p2 = Polygon([(1., 9.), (19., 9.), (19., 7.), (1., 7.)])
p = MultiPolygon([p1, p2])
result = ix.intersect_polygon(p)
result = ix.intersect(p)
assert len(result) == 2
assert result.areas.sum() == 72.
return result
@ -810,7 +810,7 @@ def test_rect_grid_polygon_with_hole():
ix = GridIntersect(gr, method="structured")
p = Polygon([(5., 5.), (5., 15.), (25., 15.), (25., -5.),
(5., -5.)], holes=[[(9., -1), (9, 11), (21, 11), (21, -1)]])
result = ix.intersect_polygon(p)
result = ix.intersect(p)
assert len(result) == 3
assert result.areas.sum() == 104.
return result
@ -827,7 +827,7 @@ def test_rect_grid_polygon_outside_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_polygon(
result = ix.intersect(
Polygon([(21., 11.), (23., 17.), (25., 11.)]))
assert len(result) == 0
return result
@ -841,7 +841,7 @@ def test_rect_grid_polygon_in_2cells_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_polygon(
result = ix.intersect(
Polygon([(2.5, 5.0), (7.5, 5.0), (7.5, 15.), (2.5, 15.)]))
assert len(result) == 2
assert result.areas.sum() == 50.
@ -856,7 +856,7 @@ def test_rect_grid_polygon_on_outer_boundary_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_polygon(
result = ix.intersect(
Polygon([(20., 5.0), (25., 5.0), (25., 15.), (20., 15.)]))
assert len(result) == 0
return result
@ -870,7 +870,7 @@ def test_rect_grid_polygon_on_inner_boundary_shapely(rtree=True):
return
gr = get_rect_grid()
ix = GridIntersect(gr, method='vertex', rtree=rtree)
result = ix.intersect_polygon(
result = ix.intersect(
Polygon([(5., 10.0), (15., 10.0), (15., 5.), (5., 5.)]))
assert len(result) == 2
assert result.areas.sum() == 50.
@ -888,7 +888,7 @@ def test_rect_grid_multipolygon_in_one_cell_shapely(rtree=True):
p1 = Polygon([(1., 1.), (8., 1.), (8., 3.), (1., 3.)])
p2 = Polygon([(1., 9.), (8., 9.), (8., 7.), (1., 7.)])
p = MultiPolygon([p1, p2])
result = ix.intersect_polygon(p)
result = ix.intersect(p)
assert len(result) == 1
assert result.areas.sum() == 28.
return result
@ -905,7 +905,7 @@ def test_rect_grid_multipolygon_in_multiple_cells_shapely(rtree=True):
p1 = Polygon([(1., 1.), (19., 1.), (19., 3.), (1., 3.)])
p2 = Polygon([(1., 9.), (19., 9.), (19., 7.), (1., 7.)])
p = MultiPolygon([p1, p2])
result = ix.intersect_polygon(p)
result = ix.intersect(p)
assert len(result) == 2
assert result.areas.sum() == 72.
return result
@ -921,7 +921,7 @@ def test_rect_grid_polygon_with_hole_shapely(rtree=True):
ix = GridIntersect(gr, method='vertex', rtree=rtree)
p = Polygon([(5., 5.), (5., 15.), (25., 15.), (25., -5.),
(5., -5.)], holes=[[(9., -1), (9, 11), (21, 11), (21, -1)]])
result = ix.intersect_polygon(p)
result = ix.intersect(p)
assert len(result) == 3
assert result.areas.sum() == 104.
return result
@ -937,7 +937,7 @@ def test_rect_grid_polygon_in_edge_in_cell(rtree=True):
ix = GridIntersect(gr, method='vertex', rtree=rtree)
p = Polygon([(0., 5.), (3., 0.), (7., 0.),
(10., 5.), (10., -1.), (0., -1.)])
result = ix.intersect_polygon(p)
result = ix.intersect(p)
assert len(result) == 1
assert result.areas.sum() == 15.
return result
@ -953,7 +953,7 @@ def test_tri_grid_polygon_outside(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_polygon(
result = ix.intersect(
Polygon([(21., 11.), (23., 17.), (25., 11.)]))
assert len(result) == 0
return result
@ -969,7 +969,7 @@ def test_tri_grid_polygon_in_2cells(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_polygon(
result = ix.intersect(
Polygon([(2.5, 5.0), (5.0, 5.0), (5.0, 15.), (2.5, 15.)]))
assert len(result) == 2
assert result.areas.sum() == 25.
@ -986,7 +986,7 @@ def test_tri_grid_polygon_on_outer_boundary(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_polygon(
result = ix.intersect(
Polygon([(20., 5.0), (25., 5.0), (25., 15.), (20., 15.)]))
assert len(result) == 0
return result
@ -997,7 +997,7 @@ def test_tri_grid_polygon_on_inner_boundary(rtree=True):
if gr == -1:
return
ix = GridIntersect(gr, rtree=rtree)
result = ix.intersect_polygon(
result = ix.intersect(
Polygon([(5., 10.0), (15., 10.0), (15., 5.), (5., 5.)]))
assert len(result) == 4
assert result.areas.sum() == 50.
@ -1017,7 +1017,7 @@ def test_tri_grid_multipolygon_in_one_cell(rtree=True):
p1 = Polygon([(1., 1.), (8., 1.), (8., 3.), (3., 3.)])
p2 = Polygon([(5., 5.), (8., 5.), (8., 8.)])
p = MultiPolygon([p1, p2])
result = ix.intersect_polygon(p)
result = ix.intersect(p)
assert len(result) == 1
assert result.areas.sum() == 16.5
return result
@ -1036,7 +1036,7 @@ def test_tri_grid_multipolygon_in_multiple_cells(rtree=True):
p1 = Polygon([(1., 1.), (19., 1.), (19., 3.), (1., 3.)])
p2 = Polygon([(1., 9.), (19., 9.), (19., 7.), (1., 7.)])
p = MultiPolygon([p1, p2])
result = ix.intersect_polygon(p)
result = ix.intersect(p)
assert len(result) == 4
assert result.areas.sum() == 72.
return result
@ -1054,7 +1054,7 @@ def test_tri_grid_polygon_with_hole(rtree=True):
ix = GridIntersect(gr, rtree=rtree)
p = Polygon([(5., 5.), (5., 15.), (25., 15.), (25., -5.),
(5., -5.)], holes=[[(9., -1), (9, 11), (21, 11), (21, -1)]])
result = ix.intersect_polygon(p)
result = ix.intersect(p)
assert len(result) == 6
assert result.areas.sum() == 104.
return result
@ -1072,7 +1072,7 @@ def test_point_offset_rot_structured_grid():
sgr = get_rect_grid(angrot=45., xyoffset=10.)
p = Point(10., 10 + np.sqrt(200.))
ix = GridIntersect(sgr, method="structured")
result = ix.intersect_point(p)
result = ix.intersect(p)
# assert len(result) == 1.
return result
@ -1086,7 +1086,7 @@ def test_linestring_offset_rot_structured_grid():
sgr = get_rect_grid(angrot=45., xyoffset=10.)
ls = LineString([(5, 10. + np.sqrt(200.)), (15, 10. + np.sqrt(200.))])
ix = GridIntersect(sgr, method="structured")
result = ix.intersect_linestring(ls)
result = ix.intersect(ls)
# assert len(result) == 2.
return result
@ -1103,7 +1103,7 @@ def test_polygon_offset_rot_structured_grid():
(15, 10. + 1.5 * np.sqrt(200.)),
(5, 10. + 1.5 * np.sqrt(200.))])
ix = GridIntersect(sgr, method="structured")
result = ix.intersect_polygon(p)
result = ix.intersect(p)
# assert len(result) == 3.
return result
@ -1117,7 +1117,7 @@ def test_point_offset_rot_structured_grid_shapely(rtree=True):
sgr = get_rect_grid(angrot=45., xyoffset=10.)
p = Point(10., 10 + np.sqrt(200.))
ix = GridIntersect(sgr, method="vertex", rtree=rtree)
result = ix.intersect_point(p)
result = ix.intersect(p)
# assert len(result) == 1.
return result
@ -1131,7 +1131,7 @@ def test_linestring_offset_rot_structured_grid_shapely(rtree=True):
sgr = get_rect_grid(angrot=45., xyoffset=10.)
ls = LineString([(5, 10. + np.sqrt(200.)), (15, 10. + np.sqrt(200.))])
ix = GridIntersect(sgr, method="vertex", rtree=rtree)
result = ix.intersect_linestring(ls)
result = ix.intersect(ls)
# assert len(result) == 2.
return result
@ -1148,7 +1148,7 @@ def test_polygon_offset_rot_structured_grid_shapely(rtree=True):
(15, 10. + 1.5 * np.sqrt(200.)),
(5, 10. + 1.5 * np.sqrt(200.))])
ix = GridIntersect(sgr, method="vertex", rtree=rtree)
result = ix.intersect_polygon(p)
result = ix.intersect(p)
# assert len(result) == 3.
return result

462
autotest/t074_test_geospatial_util.py Normal file
View File

@ -0,0 +1,462 @@
polygon = {'type': 'Polygon',
'coordinates': (((-121.389308, 38.560816),
(-121.363391, 38.568835),
(-121.358641, 38.565972),
(-121.359327, 38.562767),
(-121.369932, 38.560575),
(-121.370609, 38.557232),
(-121.385435, 38.555018),
(-121.389308, 38.560816)),)}
poly_w_hole = {'type': 'Polygon',
'coordinates': (((-121.383097, 38.565764),
(-121.342866, 38.579086),
(-121.342739, 38.578995),
(-121.323309, 38.578953),
(-121.358295, 38.561163),
(-121.379047, 38.559053),
(-121.382318, 38.562934),
(-121.383097, 38.565764)),
((-121.367281, 38.567214),
(-121.362633, 38.562622),
(-121.345857, 38.570301),
(-121.352168, 38.572258),
(-121.367281, 38.567214)))}
multipolygon = {'type': 'MultiPolygon',
'coordinates': [[((-121.433775, 38.544254),
(-121.424263, 38.547474),
(-121.422917, 38.540376),
(-121.433775, 38.544254))],
[((-121.456113, 38.552220),
(-121.440092, 38.548303),
(-121.440053, 38.537820),
(-121.459991, 38.541350),
(-121.456113, 38.552220))]]}
point = {'type': 'Point',
'coordinates': (-121.358560, 38.567760)}
multipoint = {'type': 'MultiPoint',
'coordinates': ((-121.366489, 38.565485),
(-121.365405, 38.563835),
(-121.363352, 38.566422),
(-121.362895, 38.564504),
(-121.360556, 38.565530))}
linestring = {'type': 'LineString',
'coordinates': ((-121.360899, 38.563478),
(-121.358161, 38.566511),
(-121.355936, 38.564727),
(-121.354738, 38.567047),
(-121.356678, 38.568741),
(-121.361583, 38.568072),
(-121.363066, 38.565530),
(-121.364664, 38.567359))}
multilinestring = {'type': 'MultiLineString',
'coordinates': (((-121.370653, 38.566021),
(-121.368600, 38.563255),
(-121.364207, 38.563746),
(-121.364550, 38.561605)),
((-121.370710, 38.560713),
(-121.371338, 38.561873),
(-121.370881, 38.563122),
(-121.372478, 38.563345),
(-121.373448, 38.560802),
(-121.374361, 38.562363),
(-121.373733, 38.564906),
(-121.371052, 38.567091)))}
def test_import_geospatial_utils():
from flopy.utils.geospatial_utils import GeoSpatialUtil, \
GeoSpatialCollection
return
def test_polygon():
from flopy.utils.geospatial_utils import GeoSpatialUtil
from flopy.utils.geometry import Shape, Polygon
poly = Shape.from_geojson(polygon)
gi1 = poly.__geo_interface__
if not isinstance(poly, Polygon):
raise AssertionError()
gu = GeoSpatialUtil(poly)
shp = gu.shape
shply = gu.shapely
points = gu.points
geojson = gu.geojson
fp_geo = gu.flopy_geometry
geo_types = [shp, shply, points, geojson, fp_geo]
for geo in geo_types:
if geo is None: # if shapely or geojson is not installed
continue
t = GeoSpatialUtil(geo, 'polygon').flopy_geometry
gi2 = t.__geo_interface__
is_equal = gi1 == gi2
if not is_equal:
raise AssertionError("GeoSpatialUtil polygon conversion error")
def test_polygon_with_hole():
from flopy.utils.geospatial_utils import GeoSpatialUtil
from flopy.utils.geometry import Shape, Polygon
poly = Shape.from_geojson(poly_w_hole)
gi1 = poly.__geo_interface__
if not isinstance(poly, Polygon):
raise AssertionError()
gu = GeoSpatialUtil(poly)
shp = gu.shape
shply = gu.shapely
points = gu.points
geojson = gu.geojson
fp_geo = gu.flopy_geometry
geo_types = [shp, shply, points, geojson, fp_geo]
for geo in geo_types:
if geo is None: # if shapely or geojson is not installed
continue
t = GeoSpatialUtil(geo, 'polygon').flopy_geometry
gi2 = t.__geo_interface__
is_equal = gi1 == gi2
if not is_equal:
raise AssertionError("GeoSpatialUtil polygon conversion error")
def test_multipolygon():
from flopy.utils.geospatial_utils import GeoSpatialUtil
from flopy.utils.geometry import Shape, MultiPolygon
poly = Shape.from_geojson(multipolygon)
gi1 = poly.__geo_interface__
if not isinstance(poly, MultiPolygon):
raise AssertionError()
gu = GeoSpatialUtil(poly)
shp = gu.shape
shply = gu.shapely
points = gu.points
geojson = gu.geojson
fp_geo = gu.flopy_geometry
geo_types = [shp, shply, points, geojson, fp_geo]
for geo in geo_types:
if geo is None: # if shapely or geojson is not installed
continue
t = GeoSpatialUtil(geo, 'multipolygon').flopy_geometry
gi2 = t.__geo_interface__
is_equal = gi1 == gi2
if not is_equal:
raise AssertionError("GeoSpatialUtil multipolygon "
"conversion error")
def test_point():
from flopy.utils.geospatial_utils import GeoSpatialUtil
from flopy.utils.geometry import Shape, Point
pt = Shape.from_geojson(point)
gi1 = pt.__geo_interface__
if not isinstance(pt, Point):
raise AssertionError()
gu = GeoSpatialUtil(pt)
shp = gu.shape
shply = gu.shapely
points = gu.points
geojson = gu.geojson
fp_geo = gu.flopy_geometry
geo_types = [shp, shply, points, geojson, fp_geo]
for geo in geo_types:
if geo is None: # if shapely or geojson is not installed
continue
t = GeoSpatialUtil(geo, 'point').flopy_geometry
gi2 = t.__geo_interface__
is_equal = gi1 == gi2
if not is_equal:
raise AssertionError("GeoSpatialUtil point "
"conversion error")
def test_multipoint():
from flopy.utils.geospatial_utils import GeoSpatialUtil
from flopy.utils.geometry import Shape, MultiPoint
mpt = Shape.from_geojson(multipoint)
gi1 = mpt.__geo_interface__
if not isinstance(mpt, MultiPoint):
raise AssertionError()
gu = GeoSpatialUtil(mpt)
shp = gu.shape
shply = gu.shapely
points = gu.points
geojson = gu.geojson
fp_geo = gu.flopy_geometry
geo_types = [shp, shply, points, geojson, fp_geo]
for geo in geo_types:
if geo is None: # if shapely or geojson is not installed
continue
t = GeoSpatialUtil(geo, 'multipoint').flopy_geometry
gi2 = t.__geo_interface__
is_equal = gi1 == gi2
if not is_equal:
raise AssertionError("GeoSpatialUtil multipoint "
"conversion error")
def test_linestring():
from flopy.utils.geospatial_utils import GeoSpatialUtil
from flopy.utils.geometry import Shape, LineString
lstr = Shape.from_geojson(linestring)
gi1 = lstr.__geo_interface__
if not isinstance(lstr, LineString):
raise AssertionError()
gu = GeoSpatialUtil(lstr)
shp = gu.shape
shply = gu.shapely
points = gu.points
geojson = gu.geojson
fp_geo = gu.flopy_geometry
geo_types = [shp, shply, points, geojson, fp_geo]
for geo in geo_types:
if geo is None: # if shapely or geojson is not installed
continue
t = GeoSpatialUtil(geo, 'linestring').flopy_geometry
gi2 = t.__geo_interface__
is_equal = gi1 == gi2
if not is_equal:
raise AssertionError("GeoSpatialUtil linestring "
"conversion error")
def test_multilinestring():
from flopy.utils.geospatial_utils import GeoSpatialUtil
from flopy.utils.geometry import Shape, MultiLineString
mlstr = Shape.from_geojson(multilinestring)
gi1 = mlstr.__geo_interface__
if not isinstance(mlstr, MultiLineString):
raise AssertionError()
gu = GeoSpatialUtil(mlstr)
shp = gu.shape
shply = gu.shapely
points = gu.points
geojson = gu.geojson
fp_geo = gu.flopy_geometry
geo_types = [shp, shply, points, geojson, fp_geo]
for geo in geo_types:
if geo is None: # if shapely or geojson is not installed
continue
t = GeoSpatialUtil(geo, 'multilinestring').flopy_geometry
gi2 = t.__geo_interface__
is_equal = gi1 == gi2
if not is_equal:
raise AssertionError("GeoSpatialUtil multilinestring "
"conversion error")
def test_polygon_collection():
from flopy.utils.geospatial_utils import GeoSpatialCollection
from flopy.utils.geometry import Shape, Collection
col = [Shape.from_geojson(polygon),
Shape.from_geojson(poly_w_hole),
Shape.from_geojson(multipolygon)]
gi1 = [i.__geo_interface__ for i in col]
col = Collection(col)
gc1 = GeoSpatialCollection(col)
shapetype = gc1.shapetype
shp = gc1.shape
shply = gc1.shapely
points = gc1.points
geojson = gc1.geojson
fp_geo = gc1.flopy_geometry
collections = [shp, shply, points, geojson, fp_geo]
for col in collections:
if col is None: # if geojson or shapely is not installed
continue
gc2 = GeoSpatialCollection(col, shapetype)
gi2 = [i.flopy_geometry.__geo_interface__ for i in gc2]
for ix, gi in enumerate(gi2):
is_equal = gi == gi1[ix]
if not is_equal:
raise AssertionError("GeoSpatialCollection Polygon "
"conversion error")
def test_point_collection():
from flopy.utils.geospatial_utils import GeoSpatialCollection
from flopy.utils.geometry import Shape, Collection
col = [Shape.from_geojson(point),
Shape.from_geojson(multipoint)]
gi1 = [i.__geo_interface__ for i in col]
col = Collection(col)
gc1 = GeoSpatialCollection(col)
shapetype = gc1.shapetype
shp = gc1.shape
shply = gc1.shapely
points = gc1.points
geojson = gc1.geojson
fp_geo = gc1.flopy_geometry
collections = [shp, shply, points, geojson, fp_geo]
for col in collections:
if col is None: # if geojson or shapely is not installed
continue
gc2 = GeoSpatialCollection(col, shapetype)
gi2 = [i.flopy_geometry.__geo_interface__ for i in gc2]
for ix, gi in enumerate(gi2):
is_equal = gi == gi1[ix]
if not is_equal:
raise AssertionError("GeoSpatialCollection Point "
"conversion error")
def test_linestring_collection():
from flopy.utils.geospatial_utils import GeoSpatialCollection
from flopy.utils.geometry import Shape, Collection
col = [Shape.from_geojson(linestring),
Shape.from_geojson(multilinestring)]
gi1 = [i.__geo_interface__ for i in col]
col = Collection(col)
gc1 = GeoSpatialCollection(col)
shapetype = gc1.shapetype
shp = gc1.shape
shply = gc1.shapely
points = gc1.points
geojson = gc1.geojson
fp_geo = gc1.flopy_geometry
collections = [shp, shply, points, geojson, fp_geo]
for col in collections:
if col is None: # if geojson or shapely is not installed
continue
gc2 = GeoSpatialCollection(col, shapetype)
gi2 = [i.flopy_geometry.__geo_interface__ for i in gc2]
for ix, gi in enumerate(gi2):
is_equal = gi == gi1[ix]
if not is_equal:
raise AssertionError("GeoSpatialCollection Linestring "
"conversion error")
def test_mixed_collection():
from flopy.utils.geospatial_utils import GeoSpatialCollection
from flopy.utils.geometry import Shape, Collection
col = [Shape.from_geojson(polygon),
Shape.from_geojson(poly_w_hole),
Shape.from_geojson(multipolygon),
Shape.from_geojson(point),
Shape.from_geojson(multipoint),
Shape.from_geojson(linestring),
Shape.from_geojson(multilinestring)]
gi1 = [i.__geo_interface__ for i in col]
col = Collection(col)
gc1 = GeoSpatialCollection(col)
shapetype = gc1.shapetype
shp = gc1.shape
shply = gc1.shapely
points = gc1.points
geojson = gc1.geojson
fp_geo = gc1.flopy_geometry
collections = [shp, shply, points, geojson, fp_geo]
for col in collections:
if col is None: # if geojson or shapely is not installed
continue
gc2 = GeoSpatialCollection(col, shapetype)
gi2 = [i.flopy_geometry.__geo_interface__ for i in gc2]
for ix, gi in enumerate(gi2):
is_equal = gi == gi1[ix]
if not is_equal:
raise AssertionError("GeoSpatialCollection conversion error")
if __name__ == "__main__":
test_import_geospatial_utils()
test_polygon()
test_polygon_with_hole()
test_multipolygon()
test_point()
test_multipoint()
test_linestring()
test_multilinestring()
test_polygon_collection()
test_point_collection()
test_linestring_collection()
test_mixed_collection()

132
examples/Notebooks/flopy3_grid_intersection_demo.ipynb
View File

File diff suppressed because one or more lines are too long

45
flopy/export/shapefile_utils.py
View File

@ -90,7 +90,12 @@ def write_gridlines_shapefile(filename, mg):
def write_grid_shapefile(
filename, mg, array_dict, nan_val=np.nan, epsg=None, prj=None # -1.0e9,
filename,
mg,
array_dict,
nan_val=np.nan,
epsg=None,
prj=None, # -1.0e9,
):
"""
Method to write a shapefile of gridded input data
@ -357,7 +362,8 @@ def model_attributes_to_shapefile(
for ilay in range(a.model.modelgrid.nlay):
u2d = a[ilay]
name = "{}_{}".format(
shape_attr_name(u2d.name), ilay + 1
shape_attr_name(u2d.name),
ilay + 1,
)
arr = u2d.array
assert arr.shape == horz_shape
@ -463,7 +469,12 @@ def get_pyshp_field_info(dtypename):
def get_pyshp_field_dtypes(code):
"""Returns a numpy dtype for a pyshp field type."""
dtypes = {"N": np.int, "F": np.float, "L": np.bool, "C": np.object}
dtypes = {
"N": np.int,
"F": np.float,
"L": np.bool,
"C": np.object,
}
return dtypes.get(code, np.object)
@ -480,7 +491,7 @@ def shp2recarray(shpname):
recarray : np.recarray
"""
from ..utils.geometry import shape
from ..utils.geospatial_utils import GeoSpatialCollection
sf = import_shapefile(check_version=False)
@ -489,7 +500,7 @@ def shp2recarray(shpname):
(str(f[0]), get_pyshp_field_dtypes(f[1])) for f in sfobj.fields[1:]
]
geoms = [shape(s) for s in sfobj.iterShapes()]
geoms = GeoSpatialCollection(sfobj).flopy_geometry
records = [
tuple(r) + (geoms[i],) for i, r in enumerate(sfobj.iterRecords())
]
@ -510,14 +521,18 @@ def recarray2shp(
):
"""
Write a numpy record array to a shapefile, using a corresponding
list of geometries.
list of geometries. Method supports list of flopy geometry objects,
flopy Collection object, shapely Collection object, and geojson
Geometry Collection objects
Parameters
----------
recarray : np.recarray
Numpy record array with attribute information that will go in the
shapefile
geoms : list of flopy.utils.geometry objects
geoms : list of flopy.utils.geometry, shapely geometry collection,
flopy geometry collection, shapefile.Shapes,
list of shapefile.Shape objects, or geojson geometry collection
The number of geometries in geoms must equal the number of records in
recarray.
shpname : str
@ -536,6 +551,7 @@ def recarray2shp(
subsequent use. See flopy.reference for more details.
"""
from ..utils.geospatial_utils import GeoSpatialCollection
if len(recarray) != len(geoms):
raise IndexError(
@ -546,6 +562,9 @@ def recarray2shp(
raise Exception("Recarray is empty")
geomtype = None
geoms = GeoSpatialCollection(geoms).flopy_geometry
for g in geoms:
try:
geomtype = g.shapeType
@ -722,7 +741,10 @@ class CRS(object):
if self.wktstr is not None:
sp = [
p
for p in [self.standard_parallel_1, self.standard_parallel_2]
for p in [
self.standard_parallel_1,
self.standard_parallel_2,
]
if p is not None
]
sp = sp if len(sp) > 0 else None
@ -797,7 +819,12 @@ class CRS(object):
end = s[strt:].find("]") + strt
try:
return float(self.wktstr[strt:end].split(",")[1])
except (IndexError, TypeError, ValueError, AttributeError):
except (
IndexError,
TypeError,
ValueError,
AttributeError,
):
pass
def _getgcsparam(self, txt):

457
flopy/utils/geometry.py
View File

@ -4,7 +4,307 @@ Container objects for working with geometric information
import numpy as np
class Polygon:
class Shape(object):
"""
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 = (
tuple()
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, "
"and Point: Supplied shape type {}".format(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] + [i for i in 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__
@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
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(Collection, self).__init__(geometries)
def __repr__(self):
return "Shapes: {}".format(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 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
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
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(MultiPolygon, self).__init__(polygons)
def __repr__(self):
return "MultiPolygon: {}".format(list(self))
@property
def __geo_interface__(self):
return {
"type": "MultiPolygon",
"coordinates": [g.__geo_interface__["coordinates"] for g in self],
}
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:
----------
polygons : list
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(MultiLineString, self).__init__(linestrings)
def __repr__(self):
return "LineString: {}".format(list(self))
@property
def __geo_interface__(self):
return {
"type": "MultiLineString",
"coordinates": [g.__geo_interface__["coordinates"] for g in self],
}
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:
----------
polygons : list
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(MultiPoint, self).__init__(points)
def __repr__(self):
return "MultiPoint: {}".format(list(self))
@property
def __geo_interface__(self):
return {
"type": "MultiPoint",
"coordinates": [g.__geo_interface__["coordinates"] for g in self],
}
class Polygon(Shape):
type = "Polygon"
shapeType = 5 # pyshp
@ -48,11 +348,11 @@ class Polygon:
Multi-polygons not yet supported.
z information is only stored if it was entered.
"""
self.exterior = tuple(map(tuple, exterior))
self.interiors = (
tuple()
if interiors is None
else (map(tuple, i) for i in interiors)
super(Polygon, self).__init__(
self.type,
coordinates=None,
exterior=exterior,
interiors=interiors,
)
def __eq__(self, other):
@ -80,15 +380,6 @@ class Polygon:
xmax = np.max(self._exterior_x)
return xmin, ymin, xmax, ymax
@property
def geojson(self):
return {
"coordinates": tuple(
[self.exterior] + [i for i in self.interiors]
),
"type": self.type,
}
@property
def pyshp_parts(self):
from ..export.shapefile_utils import import_shapefile
@ -96,7 +387,7 @@ class Polygon:
# exterior ring must be clockwise (negative area)
# interiors rings must be counter-clockwise (positive area)
shapefile = import_shapefile(check_version=False)
shapefile = import_shapefile()
exterior = list(self.exterior)
if shapefile.signed_area(exterior) > 0:
@ -140,21 +431,22 @@ class Polygon:
import matplotlib.pyplot as plt
except ImportError:
print("This feature requires matplotlib.")
if ax is None:
fig, ax = plt.subplots()
else:
fig = ax.figure
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)
plt.show()
except:
print("could not plot polygon feature")
return ax
class LineString:
class LineString(Shape):
type = "LineString"
shapeType = 3
has_z = False
@ -192,9 +484,7 @@ class LineString:
z information is only stored if it was entered.
"""
self.coords = list(map(tuple, coordinates))
if len(self.coords[0]) == 3:
self.has_z = True
super(LineString, self).__init__(self.type, coordinates)
def __eq__(self, other):
if not isinstance(other, LineString):
@ -227,10 +517,6 @@ class LineString:
xmax = np.max(self.x)
return xmin, ymin, xmax, ymax
@property
def geojson(self):
return {"coordinates": tuple(self.coords), "type": self.type}
@property
def pyshp_parts(self):
return [self.coords]
@ -240,18 +526,19 @@ class LineString:
import matplotlib.pyplot as plt
except ImportError:
print("This feature requires matplotlib.")
if ax is None:
fig, ax = plt.subplots()
else:
fig = ax.figure
plt.plot(self.x, self.y, **kwargs)
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)
# plt.show()
return ax
class Point:
class Point(Shape):
type = "Point"
shapeType = 1
has_z = False
@ -296,12 +583,7 @@ class Point:
-----
z information is only stored if it was entered.
"""
while len(coordinates) == 1:
coordinates = coordinates[0]
self.coords = coordinates
if len(coordinates) == 3:
self.has_z = True
super(Point, self).__init__(self.type, coordinates)
def __eq__(self, other):
if not isinstance(other, Point):
@ -334,10 +616,6 @@ class Point:
xmax = np.max(self.x)
return xmin, ymin, xmax, ymax
@property
def geojson(self):
return {"coordinates": tuple(self.coords), "type": self.type}
@property
def pyshp_parts(self):
return self.coords
@ -347,15 +625,17 @@ class Point:
import matplotlib.pyplot as plt
except ImportError:
print("This feature requires matplotlib.")
if ax is None:
fig, ax = plt.subplots()
else:
fig = ax.figure
plt.scatter(self.x, self.y, **kwargs)
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
def rotate(x, y, xoff, yoff, angrot_radians):
"""
@ -383,7 +663,13 @@ def rotate(x, y, xoff, yoff, angrot_radians):
def transform(
x, y, xoff, yoff, angrot_radians, length_multiplier=1.0, inverse=False
x,
y,
xoff,
yoff,
angrot_radians,
length_multiplier=1.0,
inverse=False,
):
"""
Given x and y array-like values calculate the translation about an
@ -439,19 +725,31 @@ def shape(pyshp_shpobj):
>>> flopy_geom = shape(list(sfobj.iterShapes())[0])
"""
types = {5: Polygon, 3: LineString, 1: Point}
flopy_geometype = types[pyshp_shpobj.shapeType]
return flopy_geometype(pyshp_shpobj.points)
import warnings
warnings.warn(
"Method will be Deprecated, calling GeoSpatialUtil",
DeprecationWarning,
)
from .geospatial_utils import GeoSpatialUtil
return GeoSpatialUtil(pyshp_shpobj).flopy_geometry
def get_polygon_area(verts):
def get_polygon_area(geom):
"""
Calculate the area of a closed polygon
Parameters
----------
verts : numpy.ndarray
polygon vertices
geom : geospatial representation of polygon
accepted types:
vertices np.array([(x, y),....])
geojson.Polygon
shapely.Polygon
shapefile.Shape
Returns
-------
@ -459,6 +757,14 @@ def get_polygon_area(verts):
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):
@ -471,14 +777,19 @@ def get_polygon_area(verts):
return a
def get_polygon_centroid(verts):
def get_polygon_centroid(geom):
"""
Calculate the centroid of a closed polygon
Parameters
----------
verts : numpy.ndarray
polygon vertices
geom : geospatial representation of polygon
accepted types:
vertices np.array([(x, y),....])
geojson.Polygon
shapely.Polygon
shapefile.Shape
Returns
-------
@ -486,6 +797,14 @@ def get_polygon_centroid(verts):
(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
@ -502,16 +821,20 @@ def get_polygon_centroid(verts):
return cx, cy
def is_clockwise(x, y):
def is_clockwise(*geom):
"""
Determine if a ring is defined clockwise
Parameters
----------
x : numpy ndarray
The x-coordinates of the ring
y : numpy ndarray
The y-coordinate of the ring
*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
-------
@ -519,6 +842,14 @@ def is_clockwise(x, y):
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[-1])

435
flopy/utils/geospatial_utils.py Normal file
View File

@ -0,0 +1,435 @@
try:
import shapely
from shapely.geometry import (
MultiPolygon,
Polygon,
Point,
MultiPoint,
LineString,
MultiLineString,
)
except ImportError:
shapely = None
try:
import geojson
except ImportError:
geojson = None
import numpy as np
from flopy.utils.geometry import Shape, Collection
geojson_classes = {}
if geojson is not None:
geojson_classes = {
"polygon": geojson.Polygon,
"multipolygon": geojson.MultiPolygon,
"point": geojson.Point,
"multipoint": geojson.MultiPoint,
"linestring": geojson.LineString,
"multilinestring": geojson.MultiLineString,
}
shape_types = {
"multipolygon": "MultiPolygon",
"polygon": "Polygon",
"point": "Point",
"multipoint": "MultiPoint",
"linestring": "LineString",
"multilinestring": "MultiLineString",
}
class GeoSpatialUtil(object):
"""
Geospatial utils are a unifying method to provide conversion between
commonly used geospatial input types
Parameters
----------
obj : geospatial object
obj can accept any of the following objects:
shapefile.Shape object
flopy.utils.geometry objects
list of vertices
geojson geometry objects
shapely.geometry objects
shapetype : str
shapetype is required when a list of vertices is supplied for obj
"""
def __init__(self, obj, shapetype=None):
from ..export.shapefile_utils import import_shapefile
self.__shapefile = import_shapefile()
self.__obj = obj
self.__geo_interface = {}
self._geojson = None
self._shapely = None
self._shape = None
self._flopy_geometry = None
self._points = None
self.__shapetype = None
if shapetype is not None:
shapetype = shapetype.lower()
if isinstance(obj, self.__shapefile.Shape):
self.__geo_interface = self.__obj.__geo_interface__
elif isinstance(obj, (Shape, Collection)):
geo_interface = obj.__geo_interface__
if geo_interface["type"] == "GeometryCollection":
raise TypeError("GeometryCollections are not supported")
self.__geo_interface = geo_interface
elif isinstance(obj, (np.ndarray, list, tuple)):
if shapetype is None or shapetype not in shape_types:
err = "shapetype must be one of the following: " + " , ".join(
geojson_classes.keys()
)
raise AssertionError(err)
self.__geo_interface = {
"type": shape_types[shapetype],
"coordinates": list(obj),
}
if geojson is not None:
if isinstance(obj, geojson.Feature):
self.__geo_interface = {
"type": obj.geometry.type,
"coordinates": obj.geometry.coordinates,
}
elif isinstance(
obj,
(
geojson.Point,
geojson.MultiPoint,
geojson.Polygon,
geojson.MultiPolygon,
geojson.LineString,
geojson.MultiLineString,
),
):
self.__geo_interface = {
"type": obj.type,
"coordinates": obj.coordinates,
}
if shapely is not None:
if isinstance(
obj,
(
Point,
MultiPoint,
Polygon,
MultiPolygon,
LineString,
MultiLineString,
),
):
self.__geo_interface = obj.__geo_interface__
@property
def __geo_interface__(self):
"""
Geojson standard representation of a geometry
Returns
-------
dict
"""
return self.__geo_interface
@property
def shapetype(self):
"""
Shapetype string for a geometry
Returns
-------
str
"""
if self.__shapetype is None:
self.__shapetype = self.__geo_interface["type"]
return self.__shapetype
@property
def points(self):
"""
Returns a list of vertices to the user
Returns
-------
list
"""
if self._points is None:
self._points = self.__geo_interface["coordinates"]
return self._points
@property
def shapely(self):
"""
Returns a shapely.geometry object to the user
Returns
-------
shapely.geometry.<shape>
"""
if shapely is not None:
if self._shapely is None:
self._shapely = shapely.geometry.shape(self.__geo_interface)
return self._shapely
@property
def geojson(self):
"""
Returns a geojson object to the user
Returns
-------
geojson.<shape>
"""
if geojson is not None:
if self._geojson is None:
cls = geojson_classes[self.__geo_interface["type"].lower()]
self._geojson = cls(self.__geo_interface["coordinates"])
return self._geojson
@property
def shape(self):
"""
Returns a shapefile.Shape object to the user
Returns
-------
shapefile.shape
"""
if self._shape is None:
self._shape = self.__shapefile.Shape._from_geojson(
self.__geo_interface
)
return self._shape
@property
def flopy_geometry(self):
"""
Returns a flopy geometry object to the user
Returns
-------
flopy.utils.geometry.<Shape>
"""
if self._flopy_geometry is None:
self._flopy_geometry = Shape.from_geojson(self.__geo_interface)
return self._flopy_geometry
class GeoSpatialCollection(object):
"""
The GeoSpatialCollection class allows a user to convert between
Collection objects from common geospatial libraries.
Parameters
----------
obj : collection object
obj can accept the following types
str : shapefile name
shapefile.Reader object
list of [shapefile.Shape, shapefile.Shape,]
shapefile.Shapes object
flopy.utils.geometry.Collection object
list of [flopy.utils.geometry, ...] objects
geojson.GeometryCollection object
geojson.FeatureCollection object
shapely.GeometryCollection object
list of [[vertices], ...]
shapetype : list
optional list of shapetypes that is required when vertices are
supplied to the class as the obj parameter
"""
def __init__(self, obj, shapetype=None):
from ..export.shapefile_utils import import_shapefile
self.__shapefile = import_shapefile()
self.__obj = obj
self.__collection = []
self._geojson = None
self._shapely = None
self._shape = None
self._flopy_geometry = None
self._points = None
self.__shapetype = None
if isinstance(obj, str):
with self.__shapefile.Reader(obj) as r:
for shape in r.shapes():
self.__collection.append(GeoSpatialUtil(shape))
elif isinstance(obj, self.__shapefile.Reader):
for shape in obj.shapes():
self.__collection.append(GeoSpatialUtil(shape))
elif isinstance(obj, self.__shapefile.Shapes):
for shape in obj:
self.__collection.append(GeoSpatialUtil(shape))
elif isinstance(obj, Collection):
for shape in obj:
self.__collection.append(GeoSpatialUtil(shape))
elif isinstance(obj, (np.ndarray, list, tuple)):
if isinstance(obj[0], (Shape, Collection, self.__shapefile.Shape)):
for shape in obj:
self.__collection.append(GeoSpatialUtil(shape))
else:
if shapetype is None:
err = "a list of shapetypes must be provided"
raise AssertionError(err)
elif isinstance(shapetype, str):
shapetype = [shapetype] * len(obj)
for ix, geom in enumerate(obj):
self.__collection.append(
GeoSpatialUtil(geom, shapetype[ix])
)
if geojson is not None:
if isinstance(
obj,
(
geojson.GeometryCollection,
geojson.FeatureCollection,
geojson.MultiLineString,
geojson.MultiPoint,
geojson.MultiPolygon,
),
):
for geom in obj.geometries:
self.__collection.append(GeoSpatialUtil(geom))
if shapely is not None:
if isinstance(
obj,
(
shapely.geometry.collection.GeometryCollection,
MultiPoint,
MultiLineString,
MultiPolygon,
),
):
for geom in list(obj):
self.__collection.append(GeoSpatialUtil(geom))
def __iter__(self):
"""
Iterator method that allows the user to get a list of GeoSpatialUtil
objects from the GeoSpatialCollection object
Returns
-------
GeoSpatialUtil
"""
yield from self.__collection
@property
def shapetype(self):
"""
Returns a list of shapetypes to the user
Returns
-------
list of str
"""
if self.__shapetype is None:
self.__shapetype = [i.shapetype for i in self.__collection]
return self.__shapetype
@property
def points(self):
"""
Property returns a multidimensional list of vertices
Returns
-------
list of vertices
"""
if self._points is None:
self._points = [i.points for i in self.__collection]
return self._points
@property
def shapely(self):
"""
Property that returns a shapely.geometry.collection.GeometryCollection
object to the user
Returns
-------
shapely.geometry.collection.GeometryCollection object
"""
if shapely is not None:
if self._shapely is None:
self._shapely = shapely.geometry.collection.GeometryCollection(
[i.shapely for i in self.__collection]
)
return self._shapely
@property
def geojson(self):
"""
Property that returns a geojson.GeometryCollection object to the user
Returns
-------
geojson.GeometryCollection
"""
if geojson is not None:
if self._geojson is None:
self._geojson = geojson.GeometryCollection(
[i.geojson for i in self.__collection]
)
return self._geojson
@property
def shape(self):
"""
Property that returns a shapefile.Shapes object
Returns
-------
shapefile.Shapes object
"""
if self._shape is None:
self._shape = self.__shapefile.Shapes()
for geom in self.__collection:
self._shape.append(geom.shape)
return self._shape
@property
def flopy_geometry(self):
"""
Property that returns a flopy.util.geometry.Collection object
Returns
-------
flopy.util.geometry.Collectionnos object
"""
if self._flopy_geometry is None:
self._flopy_geometry = Collection(
[i.flopy_geometry for i in self.__collection]
)
return self._flopy_geometry

95
flopy/utils/gridgen.py
View File

@ -40,9 +40,16 @@ def features_to_shapefile(features, featuretype, filename):
Parameters
----------
features : list
List of point, line, or polygon features
point, line, or polygon features. Method accepts
feature can be:
a list of geometries
flopy.utils.geometry.Collection object
shapely.geometry.Collection object
geojson.GeometryCollection object
list of shapefile.Shape objects
shapefile.Shapes object
featuretype : str
Must be 'point', 'line', or 'polygon'
Must be 'point', 'line', 'linestring', or 'polygon'
filename : string
name of the shapefile to write
@ -51,30 +58,42 @@ def features_to_shapefile(features, featuretype, filename):
None
"""
from .geospatial_utils import GeoSpatialCollection
shapefile = import_shapefile(check_version=True)
if featuretype.lower() not in ["point", "line", "polygon"]:
if featuretype.lower() == "line":
featuretype = "LineString"
features = GeoSpatialCollection(features, featuretype).flopy_geometry
if featuretype.lower() not in [
"point",
"line",
"linestring",
"polygon",
]:
raise Exception("Unrecognized feature type: {}".format(featuretype))
if featuretype.lower() == "line":
if featuretype.lower() in ("line", "linestring"):
wr = shapefile.Writer(filename, shapeType=shapefile.POLYLINE)
wr.field("SHAPEID", "N", 20, 0)
for i, line in enumerate(features):
wr.line(line)
wr.line(line.__geo_interface__["coordinates"])
wr.record(i)
elif featuretype.lower() == "point":
wr = shapefile.Writer(filename, shapeType=shapefile.POINT)
wr.field("SHAPEID", "N", 20, 0)
for i, point in enumerate(features):
wr.point(point[0], point[1])
wr.point(*point.__geo_interface__["coordinates"])
wr.record(i)
elif featuretype.lower() == "polygon":
wr = shapefile.Writer(filename, shapeType=shapefile.POLYGON)
wr.field("SHAPEID", "N", 20, 0)
for i, polygon in enumerate(features):
wr.poly(polygon)
wr.poly(polygon.__geo_interface__["coordinates"])
wr.record(i)
wr.close()
@ -270,8 +289,14 @@ class Gridgen(object):
Parameters
----------
feature : str or list
feature can be either a string containing the name of a polygon
shapefile or it can be a list of polygons
feature can be:
a string containing the name of a polygon
a list of polygons
flopy.utils.geometry.Collection object of Polygons
shapely.geometry.Collection object of Polygons
geojson.GeometryCollection object of Polygons
list of shapefile.Shape objects
shapefile.Shapes object
layers : list
A list of layers (zero based) for which this active domain
applies.
@ -308,9 +333,16 @@ class Gridgen(object):
"""
Parameters
----------
features : str or list
features can be either a string containing the name of a shapefile
or it can be a list of points, lines, or polygons
features : str, list, or collection object
features can be
a string containing the name of a shapefile
a list of points, lines, or polygons
flopy.utils.geometry.Collection object
a list of flopy.utils.geometry objects
shapely.geometry.Collection object
geojson.GeometryCollection object
a list of shapefile.Shape objects
shapefile.Shapes object
featuretype : str
Must be either 'point', 'line', or 'polygon'
level : int
@ -385,7 +417,11 @@ class Gridgen(object):
qtgfname = os.path.join(self.model_ws, "quadtreegrid.dfn")
if os.path.isfile(qtgfname):
os.remove(qtgfname)
cmds = [self.exe_name, "quadtreebuilder", "_gridgen_build.dfn"]
cmds = [
self.exe_name,
"quadtreebuilder",
"_gridgen_build.dfn",
]
buff = subprocess.check_output(cmds, cwd=self.model_ws)
if verbose:
print(buff)
@ -473,7 +509,11 @@ class Gridgen(object):
), "Could not create export dfn file: {}".format(fname)
# Export shapefiles
cmds = [self.exe_name, "grid_to_shapefile_poly", "_gridgen_export.dfn"]
cmds = [
self.exe_name,
"grid_to_shapefile_poly",
"_gridgen_export.dfn",
]
buff = []
try:
buff = subprocess.check_output(cmds, cwd=self.model_ws)
@ -482,7 +522,10 @@ class Gridgen(object):
fn = os.path.join(self.model_ws, "qtgrid.shp")
assert os.path.isfile(fn)
except:
print("Error. Failed to export polygon shapefile of grid", buff)
print(
"Error. Failed to export polygon shapefile of grid",
buff,
)
cmds = [
self.exe_name,
@ -497,10 +540,17 @@ class Gridgen(object):
fn = os.path.join(self.model_ws, "qtgrid_pt.shp")
assert os.path.isfile(fn)
except:
print("Error. Failed to export polygon shapefile of grid", buff)
print(
"Error. Failed to export polygon shapefile of grid",
buff,
)
# Export the usg data
cmds = [self.exe_name, "grid_to_usgdata", "_gridgen_export.dfn"]
cmds = [
self.exe_name,
"grid_to_usgdata",
"_gridgen_export.dfn",
]
buff = []
try:
buff = subprocess.check_output(cmds, cwd=self.model_ws)
@ -523,7 +573,11 @@ class Gridgen(object):
except:
print("Error. Failed to export vtk file", buff)
cmds = [self.exe_name, "grid_to_vtk_sv", "_gridgen_export.dfn"]
cmds = [
self.exe_name,
"grid_to_vtk_sv",
"_gridgen_export.dfn",
]
buff = []
try:
buff = subprocess.check_output(cmds, cwd=self.model_ws)
@ -532,7 +586,10 @@ class Gridgen(object):
fn = os.path.join(self.model_ws, "qtg_sv.vtu")
assert os.path.isfile(fn)
except:
print("Error. Failed to export shared vertex vtk file", buff)
print(
"Error. Failed to export shared vertex vtk file",
buff,
)
return

75
flopy/utils/gridintersect.py
View File

@ -6,6 +6,7 @@ except ImportError:
plt = None
from .geometry import transform
from .geospatial_utils import GeoSpatialUtil
try:
from shapely.geometry import (
@ -136,15 +137,8 @@ class GridIntersect:
if self.rtree:
self.strtree = STRtree(self._get_gridshapes())
# set interesection methods
self.intersect_point = self._intersect_point_shapely
self.intersect_linestring = self._intersect_linestring_shapely
self.intersect_polygon = self._intersect_polygon_shapely
elif self.method == "structured" and mfgrid.grid_type == "structured":
self.intersect_point = self._intersect_point_structured
self.intersect_linestring = self._intersect_linestring_structured
self.intersect_polygon = self._intersect_polygon_structured
pass
else:
raise ValueError(
@ -155,6 +149,64 @@ class GridIntersect:
)
)
def intersect(self, shp, **kwargs):
"""
Method to intersect a shape with a model grid
Parameters
----------
shp : shapely.geometry, geojson object, shapefile.Shape,
or flopy geomerty object
sort_by_cellid : bool
Sort results by cellid
keepzerolengths : bool
boolean method to keep zero length intersections for
linestring intersection
Returns
-------
numpy.recarray
a record array containing information about the intersection
"""
gu = GeoSpatialUtil(shp)
shp = gu.shapely
sort_by_cellid = kwargs.pop("sort_by_cellid", True)
keepzerolengths = kwargs.pop("keepzerolengths", False)
if gu.shapetype in ("Point", "MultiPoint"):
if (
self.method == "structured"
and self.mfgrid.grid_type == "structured"
):
rec = self._intersect_point_structured(shp)
else:
rec = self._intersect_point_shapely(shp, sort_by_cellid)
elif gu.shapetype in ("LineString", "MultiLineString"):
if (
self.method == "structured"
and self.mfgrid.grid_type == "structured"
):
rec = self._intersect_linestring_structured(
shp, keepzerolengths
)
else:
rec = self._intersect_linestring_shapely(
shp, keepzerolengths, sort_by_cellid
)
elif gu.shapetype in ("Polygon", "MultiPolygon"):
if (
self.method == "structured"
and self.mfgrid.grid_type == "structured"
):
rec = self._intersect_polygon_structured(shp)
else:
rec = self._intersect_polygon_shapely(shp, sort_by_cellid)
else:
err = "Shapetype {} is not supported".format(gu.shapetype)
raise TypeError(err)
return rec
def _set_method_get_gridshapes(self):
"""internal method, set self._get_gridshapes to the certain method for
obtaining gridcells."""
@ -557,7 +609,8 @@ class GridIntersect:
Parameters
----------
shp : shapely.geometry
shp : shapely.geometry, geojson geometry, shapefile.shape,
or flopy geometry object
shape to intersect with the grid
Returns
@ -567,6 +620,8 @@ class GridIntersect:
the shape intersects with
"""
# query grid
shp = GeoSpatialUtil(shp).shapely
qresult = self.query_grid(shp)
# filter result further if possible (only strtree and filter methods)
qfiltered = self.filter_query_result(qresult, shp)
@ -881,7 +936,7 @@ class GridIntersect:
x0 = [x[0]]
y0 = [y[0]]
(i, j) = self.intersect_point(Point(x0[0], y0[0])).cellids[0]
(i, j) = self.intersect(Point(x0[0], y0[0])).cellids[0]
Xe, Ye = self.mfgrid.xyedges
xmin = Xe[j]
xmax = Xe[j + 1]

139
flopy/utils/rasters.py
View File

@ -15,11 +15,6 @@ try:
except ImportError:
scipy = None
try:
import shapely
except ImportError:
shapely = None
class Raster(object):
"""
@ -234,17 +229,21 @@ class Raster(object):
y = np.linspace(y1, y0, ylen)
self.__xcenters, self.__ycenters = np.meshgrid(x, y)
def sample_point(self, x, y, band):
def sample_point(self, *point, band=1):
"""
Method to get nearest raster value at a user provided
point
Parameters
----------
x : float
x coordinate
y : float
y coordinate
*point : point geometry representation
accepted data types:
x, y values : ex. sample_point(1, 3, band=1)
tuple of x, y: ex sample_point((1, 3), band=1)
shapely.geometry.Point
geojson.Point
flopy.geometry.Point
band : int
raster band to re-sample
@ -252,6 +251,15 @@ class Raster(object):
-------
value : float
"""
from .geospatial_utils import GeoSpatialUtil
if isinstance(point[0], (tuple, list, np.ndarray)):
point = point[0]
geom = GeoSpatialUtil(point, shapetype="Point")
x, y = geom.points
# 1: get grid.
rxc = self.xcenters
ryc = self.ycenters
@ -282,14 +290,14 @@ class Raster(object):
Parameters
----------
polygon : (shapely.geometry.Polygon or GeoJSON-like dict)
The values should be a GeoJSON-like dict or object
implements the Python geo interface protocal.
polygon : list, geojson, shapely.geometry, shapefile.Shape
sample_polygon method accepts any of these geometries:
Alternatively if the user supplies the vectors
of a polygon in the format [(x0, y0), ..., (xn, yn)]
a single shapely polygon will be created for
cropping the data
a list of (x, y) points, ex. [(x1, y1), ...]
geojson Polygon object
shapely Polygon object
shapefile Polygon shape
flopy Polygon shape
band : int
raster band to re-sample
@ -399,14 +407,14 @@ class Raster(object):
Parameters
----------
polygon : (shapely.geometry.Polygon or GeoJSON-like dict)
The values should be a GeoJSON-like dict or object
implements the Python geo interface protocal.
polygon : list, geojson, shapely.geometry, shapefile.Shape
crop method accepts any of these geometries:
Alternatively if the user supplies the vectors
of a polygon in the format [(x0, y0), ..., (xn, yn)]
a single shapely polygon will be created for
cropping the data
a list of (x, y) points, ex. [(x1, y1), ...]
geojson Polygon object
shapely Polygon object
shapefile Polygon shape
flopy Polygon shape
invert : bool
Default value is False. If invert is True then the
@ -519,14 +527,14 @@ class Raster(object):
Parameters
----------
polygon : (shapely.geometry.Polygon or GeoJSON-like dict)
The values should be a GeoJSON-like dict or object
implements the Python geo interface protocal.
polygon : list, geojson, shapely.geometry, shapefile.Shape
_sample_rio_dataset method accepts any of these geometries:
Alternatively if the user supplies the vectors
of a polygon in the format [(x0, y0), ..., (xn, yn)]
a single shapely polygon will be created for
cropping the data
a list of (x, y) points, ex. [(x1, y1), ...]
geojson Polygon object
shapely Polygon object
shapefile Polygon shape
flopy Polygon shape
invert : bool
Default value is False. If invert is True then the
@ -546,20 +554,13 @@ class Raster(object):
else:
from rasterio.mask import mask
if shapely is None:
msg = (
"Raster()._sample_rio_dataset(): error "
+ 'importing shapely - try "pip install shapely"'
)
raise ImportError(msg)
else:
from shapely import geometry
from .geospatial_utils import GeoSpatialUtil
if isinstance(polygon, list) or isinstance(polygon, np.ndarray):
shapes = [geometry.Polygon([[x, y] for x, y in polygon])]
if isinstance(polygon, (list, tuple, np.ndarray)):
polygon = [polygon]
else:
shapes = [polygon]
geom = GeoSpatialUtil(polygon, shapetype="Polygon")
shapes = [geom]
rstr_crp, rstr_crp_affine = mask(
self._dataset, shapes, crop=True, invert=invert
@ -586,14 +587,14 @@ class Raster(object):
Parameters
----------
polygon : (shapely.geometry.Polygon or GeoJSON-like dict)
The values should be a GeoJSON-like dict or object
implements the Python geo interface protocal.
polygon : list, geojson, shapely.geometry, shapefile.Shape
_intersection method accepts any of these geometries:
Alternatively if the user supplies the vectors
of a polygon in the format [(x0, y0), ..., (xn, yn)]
a single shapely polygon will be created for
cropping the data
a list of (x, y) points, ex. [(x1, y1), ...]
geojson Polygon object
shapely Polygon object
shapefile Polygon shape
flopy Polygon shape
invert : bool
Default value is False. If invert is True then the
@ -604,36 +605,14 @@ class Raster(object):
mask : np.ndarray (dtype = bool)
"""
if shapely is None:
msg = (
"Raster()._intersection(): error "
+ 'importing shapely try "pip install shapely"'
)
raise ImportError(msg)
else:
from shapely import geometry
from .geospatial_utils import GeoSpatialUtil
# step 1: check the data type in shapes
if isinstance(polygon, geometry.Polygon):
polygon = list(polygon.exterior.coords)
if isinstance(polygon, (list, tuple, np.ndarray)):
polygon = [polygon]
elif isinstance(polygon, dict):
# geojson, get coordinates=
if polygon["geometry"]["type"].lower() == "polygon":
polygon = [
[x, y] for x, y in polygon["geometry"]["coordinates"]
]
geom = GeoSpatialUtil(polygon, shapetype="Polygon")
else:
raise TypeError("Shape type must be a polygon")
elif isinstance(polygon, np.ndarray):
# numpy array, change to a list
polygon = list(polygon)
else:
# this is a list of coordinates
pass
polygon = geom.points[0]
# step 2: create a grid of centoids
xc = self.xcenters
@ -756,8 +735,8 @@ class Raster(object):
for band, arr in self.__arr_dict.items():
foo.write(arr, band)
@classmethod
def load(cls, raster):
@staticmethod
def load(raster):
"""
Static method to load a raster file
into the raster object
@ -783,7 +762,7 @@ class Raster(object):
bands = dataset.indexes
meta = dataset.meta
return cls(
return Raster(
array,
bands,
meta["crs"],

21
flopy/utils/triangle.py
View File

@ -4,6 +4,7 @@ import subprocess
from ..mbase import which
from ..utils.cvfdutil import centroid_of_polygon
from ..plot.plotutil import plot_cvfd
from ..utils.geospatial_utils import GeoSpatialUtil
class Triangle(object):
@ -59,15 +60,29 @@ class Triangle(object):
Parameters
----------
polygon : list
polygon is a list of (x, y) points
polygon : list, geojson, shapely.geometry, shapefile.Shape
add polygon method accepts any of these geometries:
a list of (x, y) points
geojson Polygon object
shapely Polygon object
shapefile Polygon shape
flopy.utils.geometry.Polygon object
Returns
-------
None
"""
self._polygons.append(polygon)
if isinstance(polygon, (list, tuple, np.ndarray)):
polygon = [polygon]
geom = GeoSpatialUtil(polygon, shapetype="Polygon")
polygon = geom.points
self._polygons.append(polygon[0])
if len(polygon) > 1:
for hole in polygon[1:]:
self.add_hole(hole)
return
def add_hole(self, hole):