make a grid on a polygon having an angle not perpendicular

Question

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.

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

Can you show a picture or explain what you want to obtain ? Do you want your grid to be clipped by your source polygon ? Or do you want the cells in your grid to be rotated like the source polygon ? — mgc, Feb 14 '21 at 18:05
Thanks for your help I modify my question according to your comments, but the goal is to create a grid in the polygon according to its angle — patrickblancseau, Feb 14 '21 at 18:25
I looked at your method and it rotates the whole grid and therefore I end up with more polygon than with the size of a polygon having undergone a rotation — patrickblancseau, Feb 15 '21 at 08:00

Matthew Borish · Answer 1 · 2021-02-15T06:45:32.407

0

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

edited Feb 15 '21 at 06:45

answered Feb 15 '21 at 03:42

Matthew Borish

3,016
2
13
25

Thank you for your help, I will try but in the opposite direction indeed my original polygon has already undergone "minimum_rotated_rectangle" I would have preferred to start from the obtained polygon, but I tried – patrickblancseau Feb 15 '21 at 06:17
You can adjust the first parameter for .rotate() as needed. I edited my answer to illustrate this. – Matthew Borish Feb 15 '21 at 06:47
its not focused, when I try to reset the angle of the original polygon, it does not happen, it is done separately instead of a block – patrickblancseau Feb 15 '21 at 11:55
1

I'm still not sure exactly what your question is, but glad you found a solution. Would love an up vote if .rotate() helped. Thanks! – Matthew Borish Feb 15 '21 at 18:03

make a grid on a polygon having an angle not perpendicular

1 Answers1