0

I have two shapefiles one contains a polygon and on grid points inside the polygons. As shown below.

enter image description here enter image description here

I have classified grid points in N clusters. I want to place N rectangles of size 25 acres inside the polygon eg. say for 1 rectangle per 80 acres such that they do not overlap each other and one of them covers points from all clusters. I tried this method but only one rectangle of required size or none is being added.

Thanks!

TK7372
  • 49
  • 6
  • The question is a bit unclear. Can you add a picture of the expected solution? Also, can you make it clear what the N clusters are in your picture? And, of course, please show your code. – Daniel Junglas Aug 22 '22 at 14:47
  • The second picture shows example where grid points are classified in two clusters i.e. 0 and 1. Let's say we have a polygon of 80 acres then 1 rectangle of 25 acres in area should be place inside the main polygon. If polygon is of 160 acres there will 2 such rectangles. Out of all placed rectangle there should be at least one rectangle that covers some points from all present clusters. In above example the rectangle should cover few grid points from both clusters 0 and 1. – TK7372 Aug 22 '22 at 14:53

1 Answers1

0
  • using answer I previously provided, it can be used to divide a rectangle of points into .25 acre plots
  • clearly you need to consider your CRS to be able to define an area, hence have used UTM. Then number of square meters in an acre can be used to define size of .25 acres in meters
  • have synthesized your points as you have described them
import geopandas as gpd
import numpy as np
import shapely.geometry
import matplotlib.pyplot as plt
import math

fig, ax = plt.subplots(1, 2, figsize=(10, 6))

y = np.linspace(41.604, 41.611, 11)
x = np.linspace(-99.711, -99.707, 6)
Y, X = np.meshgrid(y, x)

# create the complete grid as per example
gdf = gpd.GeoDataFrame(
    geometry=gpd.points_from_xy(x=X.flatten(), y=Y.flatten()), crs="epsg:4386"
)

# cluster the subset of points
Y, X = np.meshgrid(y[3:6], x[0:3])
gdf = (
    gdf.sjoin(
        gpd.GeoDataFrame(
            geometry=gpd.points_from_xy(x=X.flatten(), y=Y.flatten()), crs="epsg:4386"
        ),
        how="left",
    )
    .rename(columns={"index_right": "cluster"})
    .assign(cluster=lambda d: (~d["cluster"].isna()).astype(int))
)

gdf.plot(column="cluster", legend=True, categorical=True, ax=ax[0])

# points to polygons
gdf_p = gdf.dissolve("cluster").envelope
gdf_p.reset_index().plot(column="cluster", legend=True, categorical=True, ax=ax[1])
gdf_p = gdf_p.to_crs(gdf_p.estimate_utm_crs())

# .25 acres expressed in meter dimensions
size_m = math.sqrt(4046.86 * 0.25)
size_m = (size_m, size_m)

for i, (c, g) in enumerate(gdf_p.iteritems()):
    rect(g, size=size_m, clip=False, include_poly=True).set_crs(gdf_p.crs).to_crs(
        gdf.crs
    ).exterior.plot(ax=ax[i], color="r")


enter image description here

Rob Raymond
  • 29,118
  • 3
  • 14
  • 30