CODEWITHSUNDEEP
pythonshapefilegeopandas
Samira Kumar
Samira Kumar

Reputation: 521

How to split a shape file vertically in python?

I've a shapefile and I'd like to split the shape into "n" vertical strips. Is there a pure python method to make this and save the output as a shape file for each vertical strip? I've the shapefile of USA and would like to split the polygon into few vertical strips.

import geopandas as gpd

gdf = gpd.read_file('cb_2017_us_nation_20m/cb_2017_us_nation_20m.shp')

AFFGEOID    GEOID   NAME    geometry
0100000US   US  United States   (POLYGON ((-136.563223 58.035052, -136.538708 ...

img

Based on the above image (sorry for poor quality), I'd like to split the polygon into "n" sub-polygon shapes.

I tried few GIS software, but couldn't come up with any solution. Any pure python method to achieve this?

Any help would be appreciated.

Upvotes: 1

Views: 1023

Answers (2)

Bera
Bera

Reputation: 1949

import geopandas as gpd
from shapely.ops import unary_union, polygonize
from shapely.geometry import LineString, Polygon
import numpy as np
import os

df = gpd.read_file(r"C:\Users\bera\Desktop\gistest\usa.shp")
n_parts = 7

xmin, ymin, xmax, ymax = df.total_bounds #Find the bound coordinates
x_step = (xmax-xmin)/n_parts
x_coords = np.arange(xmin, xmax+x_step, x_step) #List of x coordinates where to vertical split

line_list = [LineString([(x, ymin), (x, ymax)]) for x in x_coords] #Vertical lines
df_line = gpd.GeoDataFrame(geometry=line_list, crs=df.crs)

multiline = unary_union(df_line.geometry) #Create one multiline of all split lines

#Create a convex hull around the multiline and union it with the multiline, 
#   then polygonize to create n_parts number of polygons
polygon_list = list(polygonize(unary_union([multiline, multiline.convex_hull.boundary])))
df_split = gpd.GeoDataFrame(geometry=polygon_list, crs=df.crs)

ax = df_split.plot(color="red", alpha=0.2, edgecolor="black", zorder=1)
df.plot(ax=ax, color="blue", alpha=0.3, zorder=2)

enter image description here

df_split["split_id"] = range(df_split.shape[0])

intersection = gpd.overlay(df1=df, df2=df_split, how="intersection")
intersection.plot(column="split_id", cmap="Dark2")

enter image description here

out_folder = r"C:\Users\bera\Desktop\gistest\Split_us"
for split_id, subframe in intersection.groupby("split_id"):
    out_file = os.path.join(out_folder, f"Split_{split_id}.shp")
    subframe.to_file(out_file)

Upvotes: 0

Bruno Carballo
Bruno Carballo

Reputation: 1196

  1. Use a projection in meters (as opposed to a lat,lng projection)
  2. Create a geodata frame where each row is a “vertical rectangle”
  3. Get the intersection of each row in your rectangles geodata frame with the original shape file

To make the rectangles use shapely:

from shapely.geometry import Polygon

Then each rectangle is going to be made of four points which you can get from the total extent of your shapefile and dividing the total extent in “x” by the number of rectangles that you want (you’ll want to use a for loop to build each rectangle)

Upvotes: 0

Related Questions