Dividir un polígono irregular por la mitad usando las coordenadas de sus vértices en un script de Python

En un artículo anterior se determinó el área de un polígono irregular usando las coordenadas de sus vértices en un script de Python y, en este, se pretende reutilizar ese código para dividir el mismo polígono en dos partes iguales partiendo de su vértice A; tal como se encuentra en la siguiente imagen:

Intuitivamente, se percibe que el punto de corte podría ser algún punto situado sobre el lado que está entre los vértices C y D. En principio, se podría resolver el problema de manera analítica planteando un sistema de dos ecuaciones con dos incógnitas pero sucesivas divisiones para un mayor número de partes podrían resultar en una gran complejidad por lo que se prefirió un método numérico de convergencia.

Por tanto, a partir de los vértices C y D se determina su punto medio y se calcula el área para la porción superior resultante. Si ésta es mayor que el área mitad se sigue subiendo por puntos medios hasta que el área se hace menor que el área mitad para luego proceder a bajar de manera similar. Esto resulta en un pivoteo a través del punto para el cual el área quedaría dividida en dos partes iguales y el proceso se detiene en el punto de convergencia para una tolerancia que se estableció en una discrepacia menor que 1×10-6 en términos absolutos.

El código propuesto (py_divide_area.py) es el siguiente:

#!/usr/bin/env python
# -*- coding: utf-8

#Función para la estimación del área
def area(x,y,n):
	#Algoritmo para la determinación del área
	suma = x[0]*y[n-1] - x[n-1]*y[0]

	for i in range(n-1):
		suma += x[i+1]*y[i] - x[i]*y[i+1]

	return suma/2

###################################################

from os import system

system("clear")

file = open("coor.dat", "r")

coordenadas = file.read()

coordenadas = coordenadas.split("\n")

n = len(coordenadas) - 1

x = range(n)
y = range(n)

print "coordenadas x,y (archivo coor.dat)"

for i in range(n):
	x[i] = float((coordenadas[i].split())[0])
	y[i] = float((coordenadas[i].split())[1])
	print x[i], y[i]

a=area(x,y,n)

#Imprime el área determinada
print "\narea = %.2f m2\n" % a

area1_2 = a/2

n=4
x=x[:4]
y=y[:4]

iterador=0

x1 = x[2]
x2 = x[3]
y1 = y[2]
y2 = y[3]

control = 1 

while(control > 1e-6):

	iterador += 1

	#Cálculo punto medio de dos puntos

	x_med = (x1 + x2)/2
	y_med = (y1 + y2)/2

	x[3] = x_med
	y[3] = y_med

	b=area(x,y,n)

	delta = area1_2 - b	

	if delta < 0:

		x_old = x2
		y_old = y2

		x2 = x_med
		y2 = y_med

	if delta > 0:

		x1 = x_med 
		y1 = y_med

		x2 = x_old
		y2 = y_old

	control = abs(delta)

print "iteración %d: área mitad: %.5f m2 área final iterada: %.5f m2" % (iterador, area1_2, b)
print "x1=%.5f  y1=%.5f  x2=%.5f  y2=%.5f" % (x1,y1,x2,y2)

y produce, al ejecutarlo en cónsola de bash, esta salida:

coordenadas x,y (archivo coor.dat)
0.0 591.66
125.66 847.62
716.3 694.06
523.58 0.0
517.54 202.97

area = 272599.90 m2

iteración 54: área mitad: 136299.95195 m2 área final iterada: 136299.95195 m2
x1=675.07596  y1=545.59611  x2=675.07596  y2=545.59611

En la última línea se encuentran las coordenadas del punto de convergencia que resultó ser:

(xc, yc) = (675.07596, 545.59611).

Para probar si el programa trabajaba bien convertí rápidamente ambas semi áreas en shapefiles de puntos con este sencillo programa en python:

Script en Python, con librería ogr, para crear un Point shapefile y añadirle puntos

y luego, las cargué en QGIS donde las convertí en polígonos con el plugin Point2One y en modo de edición le añadí una columna de área en metros cuadrados. En la imagen siguiente se puede corroborar que el programa logró su cometido con error menor de 0,01 m2.

Pequeños ajustes permitirían dividir en áreas que correspondan a 1/3, 1/4, etc y usar otros puntos como partida de la división ya que las opciones son infinitas.

Esta entrada fue publicada en Código Python, QGIS, SIG, Software Libre. Guarda el enlace permanente.

8 respuestas a Dividir un polígono irregular por la mitad usando las coordenadas de sus vértices en un script de Python

  1. Pingback: Analizar geometría de un shapefile mediante el módulo OGR de Python |

  2. Alfonso dijo:

    Hola
    Actualmente en mi trabajo de tesis intento hallar áreas pequeñas de diferentes zonas de una finca, ya he revisado tutoriales en youtube pero ninguno me ha funcionado.

    intento hallar el área de un poligono que se por prueba de campo de debe dar alrededor de 3000 m2… pero de ninguna forma logro un resultado similar.

    Solo el polígono general me da 0.000008 supongo yo son km2, se acerca a la realidad pero no logro obtener un resultado plausible en los polígonos pequeños.

    agradecería mucho ayuda para este aspecto, ¿que debo configurar? ¿que me hace falta?

    gracias.

    • Alfonso dijo:

      Olvidaba comentar el programa que estoy utilizando es Quantum Gis 2.0.

      Gracias.

      • Los valores son típicos de que parece que estás intentando calcular áreas en un proyecto definido en coordenadas geográficas. Tienes que reproyectar en metros para que los parámetros de geometría tengan sentido.

  3. Manuel dijo:

    Buenas:

    Yo tengo un problema en el que llevo inmerso unos días y no consigo solucionar.
    Tengo que dividir una serie de polígonos en diferentes subpolígonos, pero lo tengo que realizar de forma automática ya que son muchísimos. Es decir, tengo que eliminar 216.000 polígonos por su mitad. Suelo trabajar con formato .shp de ARCGIS. Ruego me echéis una mano.

    • Yo no manejo ArcGis desde hace un buen tiempo pero sé que hay una extensión que se llama Parcel Editor que hace lo que tu quieres. Para automatizarlo vas a tener que aprender arcgis.scripting con Python.

  4. Estimado, necesito que me ayuden con un problema que no puedo resolver. Tengo un polígono y si no me equivoco es “multiparte” al cual lo quiero dividir con la herramienta “Dividir objetos espaciales” y no me permite realizarlo. Lo hice tanto seleccionando el polígono como no. Es un polígono muy complejo, pero no encuentro la forma de dividirlo. Alguna idea?. Mil gracias

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s