4

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))
Vince
  • 20,017
  • 15
  • 45
  • 64
hashalluring
  • 119
  • 12
  • 1
    So the problem is you get a geoseries when you want a dataframe? How many multipolygons do you want to generate random points for? How important is it to get a certain number of points? – BERA Sep 04 '20 at 08:04
  • No I'm just getting the boolen values not the points in a geodataframe(gdf_points). I have 30+ polygon geometries inside the shapefile and I just want few points inside each geometry not necessary to have certain number of points. – hashalluring Sep 04 '20 at 08:13
  • You would be fine if some polygons don't have any point? – til_b Sep 04 '20 at 09:07
  • No I need points inside each polygon @til_b – hashalluring Sep 04 '20 at 09:42

2 Answers2

2

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
Vince
  • 20,017
  • 15
  • 45
  • 64
gene
  • 54,868
  • 3
  • 110
  • 187
1

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
vbferreira
  • 11
  • 1