TransWikia.com

Generate points inside multipolygons using GeoPandas

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

One Answer

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

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP