在角度不垂直的多边形上制作网格

问题描述

您好，我可以轻松地在矩形多边形上创建网格。另一方面，如何在旋转的多边形上制作网格？

目标和多边形按照相同的对齐方式切割成网格示例。

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+lenght),(x,y+lenght)]) )
grid = gpd.GeoDataFrame({'geometry':polygons},crs = gdf.crs.to_string())
grid['CASE'] = "FILLED"
print(grid)

预先感谢您的想法。

经过长时间搜索和修改我的代码后，我找到了一个有效的解决方案。谢谢大家的帮助。

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,Linestring


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


print(essias)


d = {'geometry': [polygon([(32250,crs="epsg:4326")

#xmin,ymax = essias.total_bounds

#print(xmin,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,ymax=gdf1.total_bounds
lenght = 2500
wide = 3500

cols = list(range(int(np.floor(xmin)),lenght))
rows.reverse()
polygons = []

for x in cols:
    for y in rows:
        polygons.append( polygon([(x,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,"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")

解决方法

您可以使用 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())

geopandas python shapely