Puntos en polígonos utilizando el algoritmo “Ray Casting” en PyQGIS

Mirando en el Blog de Joel Lawhead encontré el algoritmo “Ray Casting”, ya escrito en Python, que permite evaluar el número de puntos en polígonos. Mi interés es adaptarlo para que funcione en el entorno de programación de PyQGIS y para ello le hice las modificaciones que propongo a continuación:

# Determine if a point is inside a given polygon or not
# Polygon is a list of (x,y) pairs. This function
# returns True or False.  The algorithm is called
# the "Ray Casting Method".

def point_in_poly(point,poly):
    x = point.x()
    y = point.y()

    n = len(poly)
    inside = False

    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x,p1y = p2x,p2y

    return inside

## Test
layer = iface.activeLayer()

iter = layer.getFeatures()

features = []

for feature in iter:
    geom = feature.geometry().asPolygon()
    features.append(geom)

points = [QgsPoint(411795.923604, 4448694.05657), 
          QgsPoint(386974.02657, 4447841.8026),
          QgsPoint(411245.830677, 4442551.55703),
          QgsPoint(412760.409476, 4454668.18741),
          QgsPoint(436387.838734, 4453456.52438),
          QgsPoint(447595.721844, 4447095.29342)]

## Call the function with the points and the polygon
count = [0]*(layer.featureCount())

for point in points:
    i = 0
    for feat in features:
        if point_in_poly(point, feat[0]) == True:
            count[i] += 1
        i += 1

print count

La lista de QgsPoint que se tiene dentro del código corresponde a los que se observan en la imagen siguiente:

point_polygon2

Visualmente se detecta que en uno de los features se tienen 3 puntos y en el otro 2. La ejecución del código anterior en la Python Console de QGIS permite obtener la lista [3, 2]; tal como era de esperar. Para probar un caso algo más complejo vamos a utilizar una grid de 12 filas por 25 columnas, es decir, 300 features; como en la imagen a continuación:

point_polygon1

A la lista de salida en la Python Console de la imagen le di forma de matriz de 12 x 25 para que se pudiera observar mejor que el conteo de los puntos en las celdas de la grid también fue el esperado.

Aunque pudiese parecer lo contrario, programar este código como plugin es aún más sencillo que en la Python Console porque se pueden usar filtros específicos en las QgsMapLayerComboBox que seleccionarán los vectoriales tipo punto y polígono que facilitan la programación. Por otra parte, con las bondades del “Data Provider”, introduciríamos la lista que contiene el conteo de puntos como un campo adicional en el vectorial tipo polígono.

En nuestro curso de PyQGIS se contemplan con detalle el desarrollo de éste y otros plugins mediante el empleo de ‘Plugin Builder’.

PRÓXIMA SESIÓN EL 20 DE OCTUBRE

Esta entrada fue publicada en PyQGIS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

Una respuesta a Puntos en polígonos utilizando el algoritmo “Ray Casting” en PyQGIS

  1. Pingback: Newsletter Octubre 2015 | El Blog de José Guerrero

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