Cómo crear cuando sea posible un rectángulo en polígonos con agujeros

Question

No he encontrado una pregunta similar en este sitio. ¿Hay alguna manera de crear automáticamente un polígono rectangular de 500 metros cuadrados dentro de otro polígono con anillos utilizando QGIS? (Tengo 28.000 polígonos con anillos)

Especifico mi pregunta: mi trabajo consiste en estudiar una zona urbanizada (véase A) .

En algunos casos, el terreno aún puede albergar construcciones. Para encontrar estas situaciones, hice un amortiguador alrededor de los marcos con un retroceso

(véase B ),

Corté estos topes de la tierra

(véase C )

Los objetos obtenidos (C) tienen una superficie total que no es relevante para tener en cuenta su edificabilidad (debido a los agujeros). Se trata de insertar, cuando sea posible, en cada objeto de C, un trozo de terreno rectangular con una superficie determinada (por ejemplo: 500 m² , cf rectángulo rojo).

¿Puedes darme una pista?

Answer 1

Cree puntos aleatorios en cada parcela, tamponelos, cree rectángulos a partir de los tampones utilizando st_expand, seleccione los rectángulos dentro de cada parcela que no intersecten con ningún edificio:

--Create n random points in each plot
create table pointsinplot as
with cte as (
select (st_dump(st_generatepoints(geom, 10000))).geom from
(select distinct p.geom from plot p
left join
building b
on st_within(st_centroid(b.geom), p.geom)) sub
)
select row_number() over() as rnum, cte.geom from cte
left join building b
on st_intersects(cte.geom, b.geom)
where b.geom is null --no points allowed inside the buildings

--Buffer each point, expand the buffers into a rectangle by random x an y factors and rotate by random radians
create table randomrectangles as
with cte as (
select st_expand(st_buffer(geom, 2),floor(random() * 10 + 1)::int,floor(random() * 10 + 1)::int) geom from public.pointsinplot
)
select row_number() over() as rnum, st_rotate(geom, random()*2*pi(), st_centroid(geom)) geom from cte

--Find the rectangles within plots, not intersecting the buildings
drop table if exists potential_sites;
create table potential_sites as
with cte as (
select plotid, row_number() over() as rnum, r.geom from randomrectangles r
join plot p
on st_within(r.geom, p.geom)
)
select cte.plotid, cte.rnum, cte.geom, st_area(cte.geom) from cte
left join building b
on st_intersects(cte.geom, b.geom)
where b.geom is null

--Find largest rectangle per plot
SELECT DISTINCT ON (plotid)
       plotid, geom
FROM   potential_sites
ORDER  BY plotid, st_area(geom) DESC;

Answer 2

Principios: qué puede esperar de esta solución

Esta solución es un enfoque heurístico. Devuelve las parcelas en las que seguro se puede colocar un rectángulo de 500 m2 así como otras parcelas "candidatas" con una alta probabilidad de que sea posible, pero que deben comprobarse manualmente. Véase la captura de pantalla:

Resultado: edificio inicial (gris) con búfer (wihte) y parcela restante (rojo). Los dos edificios posibles de 500 m² (sombreados en gris con contorno negro). El inferior está completamente dentro de la parcela, por lo que se reconoce automáticamente como solución; el superior no está completamente, pero casi dentro de la parcela (sólo un 3% de su superficie fuera de la parcela), por lo que es un candidato a tener en cuenta y, de hecho, desplazando un poco el polígono, quedaría completamente dentro de la parcela:

La solución paso a paso

1. Obtener la parcela que queda al recortar los edificios amortiguados (utilizando Difference herramienta).

2. Corta estrechos "pasillos" de conexión entre las formas principales mediante un amortiguamiento con un valor negativo de, digamos, -5 m, y aplica a continuación el valor inverso, positivo, para el amortiguamiento: 5 m.

Izquierda: edificio 8gris), búfer (amarillo), parcela restante (rojo); centro: parcela restante (rojo), búfer negativo (azul); derecha: parcela restante (delineada en rojo), búfer negativo (delineado en azul), las dos parcelas resultantes que utilizaremos a partir de ahora (naranja sólido):

3. Convertir de Multipart a partes individuales y añadir un id único (aquí: fid ) al resultado.

4. Utilice Menu Processing / Toolbox / Minimum bounding geometry y crear el Minimum oriented rectangle basado en el fid para obtener una caja separada para cada característica.

5. Escala los polígonos del paso 4 con un factor de escala para que tengan el tamaño que desees (500 metros cuadrados, en tu caso). El factor de escala es el square root of (500 divided by the area of the polygon) (pseudocódigo).

   Las parcelas de arriba (naranja), el rectángulo orientado mínimo del paso 4 (azul) y estos rectángulos reducidos a 500 m2 (rojo) - aquí utilizando el generador de Geometría:

   La expresión 1 (ver abajo) devuelve el Minimum oriented rectangle con una superficie de 500 m2 y el mismo centroide.

6. Ahora calcule la intersección del polígono escalado (paso 5) con el polígono de la parcela del paso 3. Si el polígono a escala de 500 m2 está completamente dentro de la parcela, la intersección es la misma que la salida del paso 5. Si está no completamente dentro de ella, algunas partes están cortadas. Calcule ahora el área de esta intersección con la expresión 2 de abajo (una extensión de la expresión 1) creando un nuevo campo de atributo con la calculadora Field.

7. Ahora seleccione todas las características en las que el campo calculado en el paso 6 esté cerca del valor de 500 m2. Lo ideal sería buscar sólo los valores que tengan exactamente 500 m². Pero deberíamos restar un cierto umbral de tolerancia para tener en cuenta el hecho de que, a veces, desplazar el polígono escalado dentro de la parcela podría dar lugar a una intersección mayor. Por tanto, en función de sus necesidades y del tipo de datos que tenga, seleccione todos los valores que se encuentren en un determinado intervalo (por ejemplo, de 450 a 500).

Expresiones

Expresión 1

 with_variable(
     'scale',
     sqrt (500/$area),
     make_rectangle_3points(
         project (
             centroid ($geometry),
             @scale*length (make_line (point_n ($geometry, 1), centroid ($geometry))),
             azimuth (point_n ($geometry, 1), centroid ($geometry))
         ),
         project (
             centroid ($geometry),
             @scale*length (make_line (point_n ($geometry, 2), centroid ($geometry))),
             azimuth (point_n ($geometry, 2), centroid ($geometry))
         ),
         project (
             centroid ($geometry),
             @scale*length (make_line (point_n ($geometry, 2), centroid ($geometry))),
             azimuth (point_n ($geometry, 3), centroid ($geometry))
         )
     )
 )

Expresión 2

 area (
     intersection (
         with_variable(
             'scale',
             sqrt (5000/area($geometry)),
             make_rectangle_3points(
                 project (
                     centroid ($geometry),
                     @scale*length (make_line (point_n ($geometry, 1), centroid ($geometry))),
                     azimuth (point_n ($geometry, 1), centroid ($geometry))
                 ),
                 project (
                     centroid ($geometry),
                     @scale*length (make_line (point_n ($geometry, 2), centroid ($geometry))),
                     azimuth (point_n ($geometry, 2), centroid ($geometry))
                 ),
                 project (
                     centroid ($geometry),
                     @scale*length (make_line (point_n ($geometry, 2), centroid ($geometry))),
                     azimuth (point_n ($geometry, 3), centroid ($geometry))
                 )
             )
         ),
         overlay_nearest ('parcel',$geometry)[0]
     )
 )