Geographic Information Systems Asked on December 27, 2020
I have a shapefile which has multipolygon geometries. I want to generate random points inside them using GeoPandas. I tried with the below code but it’s not giving geodataframe of points.
df = gpd.read_file(multipolygonshp)
geom= df['geometry']
final_geom = geom.unary_union
x_min, y_min, x_max, y_max = final_geom.bounds
n = 3000
x = np.random.uniform(x_min, x_max, n)
y = np.random.uniform(y_min, y_max, n)
gdf_points = gpd.GeoSeries(gpd.points_from_xy(x, y))
gdf_points = gpd.GeoSeries(gdf_points.within(final_geom.unary_union))
<ipython-input-125-ed976bd8ef24>:23: FutureWarning: You are passing non-geometry data to the GeoSeries constructor. Currently,
it falls back to returning a pandas Series. But in the future, we will start
to raise a TypeError instead.
gdf_points = gpd.GeoSeries(gdf_points.within(final_geom))
FutureWarning is a simple warnings statement that you can eliminate with
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
The result of gdf_points.within
is a Boolean, as you say (shapely binary predicates)
gdf_points = gpd.GeoSeries(gdf_points.within(final_geom))
gdf_points.head(5)
0 False
1 False
2 True
3 False
4 False
dtype: bool
And it is a Pandas serie (not a GeoSerie -> no geometry)
type(gdf_points)
pandas.core.series.Series
There are many solutions in GIS SE and I use here Geopandas equivalent to select by location
gdf_points = gpd.GeoSeries(gpd.points_from_xy(x, y))
subset = gdf_points[gdf_points.within(final_geom)]
subset.head(5)
2 POINT (0.49117 -0.09817)
5 POINT (-0.51586 -0.68984)
8 POINT (-0.42499 0.01249)
10 POINT (-0.21973 -0.35859)
11 POINT (-0.77146 0.17700)
dtype: geometry
type(subset)
geopandas.geoseries.GeoSeries
Look also Check if a point falls within a multipolygon with Python with GeoDataFrames, for example
gdf_points = gpd.GeoSeries(gpd.points_from_xy(x, y))
points = gpd.GeoDataFrame({'geometry':gdf_points})
poly = gpd.GeoDataFrame({'geometry':final_geom})
within_points = gpd.sjoin(points, poly, op = 'within')
within_points.head(5)
geometry index_right
2 POINT (0.49117 -0.09817) 0
14 POINT (0.66375 -0.09176) 0
17 POINT (0.55135 -0.36559) 0
39 POINT (0.69521 -0.37174) 0
59 POINT (0.60061 -0.29019) 0
Answered by gene on December 27, 2020
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP