make a grid on a polygon having an angle not perpendicular

Viewed 486

Hello I have no difficulty creating a grid on a rectangular polygon. On the other hand how can I make a grid on a polygon undergoing a rotation?

The goal and that the polygon is cut into a grid following the same alignment example.

import geopandas as gpd
import numpy as np

from shapely.geometry import Point, Polygon, LineString
d = {'geometry': [Polygon([(32250, 175889), (33913, 180757), (29909, 182124), (28246, 177257), (32250, 175889)])]}

gdf = gpd.GeoDataFrame(d, crs="EPSG:4326")

#creating a grid on a rectangular
#on the other hand with an angle???????????????

xmin,ymin,xmax,ymax=gdf.total_bounds
lenght = 2500
wide = 3500

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 = []
polygons_minimum = []
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['CASE'] = "FILLED"
print(grid)

Thank you in advance for your ideas.

enter image description here

enter image description here

I found a solution that works after a long search and modification of my code. Thank you for your help everyone.

import csv 
import glob
import pandas as pd
import numpy as np
import time
import datetime
import sys, math
import geopandas as gpd
from shapely.geometry import Polygon
from shapely.geometry import Point
import shapely
from geopandas import GeoSeries
from shapely.geometry import Point, Polygon, LineString


essias=Polygon([[32250, 175889], [33913, 180757], [29909, 182124], [28246, 177257], [32250, 175889]])


print(essias)


d = {'geometry': [Polygon([(32250, 175889), (33913, 180757), (29909, 182124), (28246, 177257), (32250, 175889)])]}

gdf = gpd.GeoDataFrame(d, crs="EPSG:4326")

#xmin,ymin,xmax,ymax = essias.total_bounds

#print(xmin,ymin,xmax,ymax)


gdf_loc1 = gdf.iloc[0].geometry        
l = gdf_loc1.boundary
coords = [c for c in l.coords]
time.sleep(1)
segments = [shapely.geometry.LineString([a, b]) for a, b in zip(coords,coords[1:])]
longest_segment = max(segments, key=lambda x: x.length)
p1, p2 = [c for c in longest_segment.coords]
anglest = math.degrees(math.atan2(p2[1]-p1[1], p2[0]-p1[0])) #https://stackoverflow.com/questions/42258637/how-to-know-the-angle-between-two-points

edge_length=(Point(coords[0]).distance(Point(coords[1])),Point(coords[1]).distance(Point(coords[2])))
length = max(edge_length)
width = min(edge_length)

gdf1 = gdf.rotate(-anglest, origin=gdf.centroid.item())
gdf_loc = gdf1.iloc[0]        
df4_merged_geom = gdf1.cascaded_union
l1 = gdf1.boundary
xmin,ymin,xmax,ymax=gdf1.total_bounds
lenght = 2500
wide = 3500

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['CASE'] = "FILLED"

############################################
for index, polygon in grid.iterrows():
    if polygon.geometry.disjoint(df4_merged_geom):
        grid.loc[index, "CASE"] = "EMPTY"
grid_plein=grid[grid['CASE'].str.contains("FILLED")]
print("grid['CASE'].str.contains(FILLED)",grid['CASE'].str.contains("FILLED"))
grid_plein=grid_plein.reset_index(drop=True)

for index, polygon in grid.iterrows():
    if polygon.geometry.disjoint(df4_merged_geom):
        grid.loc[index, "CASE"] = "EMPTY"
grid_plein=grid[grid['CASE'].str.contains("FILLED")]
############################################


print("grid['CASE'].str.contains(FILLED)",grid['CASE'].str.contains("FILLED"))
grid_plein=grid_plein.reset_index(drop=True)

gdf_apres = grid_plein.rotate(anglest, origin=gdf.centroid.item())
gdf_loc_apres = gdf_apres.iloc[1] 
gdf_apres.to_file("dossier_VDR/grid_comple.shp")
1 Answers

You can use gpd .rotate().

# dissolve initial polys to obtain centroid for origin param
grid_d = grid.dissolve(by='CASE')

# make a copy of grid before rotation
grid_r = grid.copy(deep=True)

#rotate polys using grid_d centroid as origin
grid_r['geometry'] = grid_r['geometry'].rotate(-30, origin=grid_d.centroid.item())

enter image description here

Related