Generate points inside multipolygons using GeoPandas

Question

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))

Answer 1

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

Answer 2

You may also try to buffer the df to avoid self-intersection in multi-polygons layers:

df = gpd.read_file(multipolygonshp)

# buffering to avoid self-intersection in multipolygons layers
df = geom.buffer(0.1)

# Use GeoSeries.total_bounds if you would like the limits of the entire series
x_min, y_min, x_max, y_max = df.total_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))

# You can directly retrieve the points that intersect the polygon
gdf_points = gdf_points[gdf_points.within(gdf_points.unary_union)]

gdf_points

0       POINT (318400.364 5064051.518)
1       POINT (264133.118 5067828.107)
2       POINT (307008.396 5090454.061)
3       POINT (287648.343 5091872.906)
4       POINT (275455.794 5068725.441)
                     ...              
2995    POINT (265071.293 5082035.520)
2996    POINT (300584.749 5085800.227)
2997    POINT (316275.872 5093292.001)
2998    POINT (308961.236 5093976.585)
2999    POINT (283761.834 5062090.378)
Length: 3000, dtype: geometry