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))
4
  • 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 4, 2020 at 8: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. Sep 4, 2020 at 8:13
  • You would be fine if some polygons don't have any point?
    – til_b
    Sep 4, 2020 at 9:07
  • No I need points inside each polygon @til_b Sep 4, 2020 at 9:42

2 Answers 2

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
2
  • Filtering the warning is not a solution, because it will return an error in the next release. Sep 4, 2020 at 12:23
  • @gene I tried the above methods as suggested by you, it is hardly give me any points(like 5-10) put together. Sep 4, 2020 at 12:45
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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.