9

I have a shapefile and I would like the exploded view with GeoPandas in several polygons, I cannot use ArcGIS and keep the polygons that have a given.

import csv 
import glob
import time
import sys
import datetime
import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Polygon

gdf = gpd.read_file("dossier/essais.shp")

time.sleep(10)
xmin,ymin,xmax,ymax = gdf.total_bounds

lenght = 2065
wide = 3069

cols = list(range(int(np.floor(xmin)), int(np.ceil(xmax)), wide))
rows = list(range(int(np.floor(ymin)), int(np.ceil(ymax)), lenght))
rows.reverse()

polygons = []
for x in cols:
    for y in rows:
        polygons.append( Polygon([(x,y), (x+wide, y), (x+wide, y+lenght), (x, y+lenght)]) )

grid = gpd.GeoDataFrame({'geometry':polygons}, crs = gdf.crs.to_string())
grid.to_file("dossier/grid.shp")

#df4.head(1000).to_file("dossier/"+file)
gdf.to_file("dossier/grille.shp")

enter image description here

4

1 Answer 1

10
gdf = gpd.read_file("dossier/essais.shp")
gdf_merged_geom = gdf.unary_union

#
# other lines
#

grid = gpd.GeoDataFrame({'geometry':polygons}, crs=gdf.crs.to_string())

# add new CASE column
grid["CASE"] = "FILLED"

for index, polygon in grid.iterrows():
    if polygon.geometry.disjoint(gdf_merged_geom):
        grid.loc[index, "CASE"] = "EMPTY"

print(grid.head())

# OUTPUT
#           geometry      CASE
# 0 POLYGON ((282...    EMPTY
# 1 POLYGON ((282...    FILLED
# 2 POLYGON ((282...    FILLED
# 3 POLYGON ((282...    FILLED
# 4 POLYGON ((282...    FILLED

Sample data:

enter image description here

Result:

custom_draw(grid[grid["CASE"]=="EMPTY"].geometry)
# OR 
# grid[grid["CASE"]=="EMPTY"].plot()

enter image description here

2
  • 3
    Please use unary_union instead of cascaded_union as recommended in the Shapely manual (The unary union function can operate on different geometry types, not only polygons as is the case for the older cascaded union.)
    – gene
    Jan 23, 2021 at 12:18
  • @gene Done. Thank you. I misunderstood the note in the shapely manual. Jan 23, 2021 at 15:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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