Is your feature request related to a problem? Please describe.
From discussion in #2646, shapely supports operations on numpy arrays with geometries, but these are not really supported by GeoSpatialUtil.
Describe the solution you'd like
import shapely
from flopy.utils.geospatial_utils import GeoSpatialUtil
x = np.ones(100)
y = np.ones(100)
pts = shapely.points(x, y)
gu = GeoSpatialUtil(points)
gu.shapely # np.array of points
Currently the above code raises AssertionError: shapetype must be one of the following:
This will allow users to pass arrays of geometries to GridIntersect operations.