¿Cómo determinar los identificadores de mosaico vecinos en QGIS?

En un curso de capacitación reciente me preguntaron si QGIS podía calcular automáticamente los números de página siguiente / anterior y superior / inferior para un libro de mapas creado usando el generador de atlas. Logré elaborar una expresión de etiqueta bastante razonable para una cuadrícula regular si conoces el ancho y la altura de la cuadrícula.

Pero luego comenzamos a pensar en ejemplos realistas en los que no queremos dibujar páginas que no contengan nuestro distrito de interés, como este de mi condado de origen:

Así que esta tarde jugué en un script de Python para resolver los 4 vecinos que me interesaban para cada celda de la cuadrícula y agregué esos valores a mi cuadrícula (esto se basa en gran medida en el tutorial de Ujaval Gandhi ):

for f in feature_dict.values():
    print 'Working on %s' % f[_NAME_FIELD]
    geom = f.geometry()
    # Find all features that intersect the bounding box of the current feature.
    # We use spatial index to find the features intersecting the bounding box
    # of the current feature. This will narrow down the features that we need
    # to check neighboring features.
    intersecting_ids = index.intersects(geom.boundingBox())
    # Initalize neighbors list and sum
    neighbors = []
    neighbors_sum = 0
    for intersecting_id in intersecting_ids:
        # Look up the feature from the dictionary
        intersecting_f = feature_dict[intersecting_id]
        int_geom = intersecting_f.geometry()
        centroid = geom.centroid()
        height = geom.boundingBox().height()
        width = geom.boundingBox().width()
        # For our purpose we consider a feature as 'neighbor' if it touches or
        # intersects a feature. We use the 'disjoint' predicate to satisfy
        # these conditions. So if a feature is not disjoint, it is a neighbor.
        if (f != intersecting_f and
            not int_geom.disjoint(geom)):
            above_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x(),
               centroid.asPoint().y()+height))
            below_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x(),
               centroid.asPoint().y()-height))
            left_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x()-width,
               centroid.asPoint().y()))
            right_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x()+width,
               centroid.asPoint().y()))
            above = int_geom.contains(above_point)   
            below = int_geom.contains(below_point)   
            left = int_geom.contains(left_point)
            right = int_geom.contains(right_point)
            if above:
                print "setting %d as above %d"%(intersecting_f['id'],f['id'])
                f['above']=intersecting_f['id']

            if below:
                print "setting %d as below %d"%(intersecting_f['id'],f['id'])
                f['below']=intersecting_f['id']

            if left:
                print "setting %d as left of %d"%(intersecting_f['id'],f['id'])
                f['left']=intersecting_f['id']

            if right:
                print "setting %d as right of %d"%(intersecting_f['id'],f['id'])
                f['right']=intersecting_f['id']

    # Update the layer with new attribute values.
    layer.updateFeature(f)

layer.commitChanges()

Esto funciona bien

Pero para ser honesto, la creación de un punto de prueba hacia el norte y luego probar todos los posibles vecinos parece estar mal. Sin embargo, después de una tarde de sacudir mi cerebro, no puedo pensar en una mejor manera de determinar cuál es el vecino del norte de una celda de red en particular.

Idealmente, me gustaría algo lo suficientemente simple como para ponerlo en un cuadro de texto del compositor de impresión, pero sospecho que es demasiado pedir.

Si no está ajustando cada extensión de página (desde la capa de índice) exactamente en el compositor, sino que tiene bordes superpuestos con páginas adyacentes (como se muestra en su segunda captura de pantalla), entonces podría usar etiquetas de la capa de índice, con la desventaja de que estarían dentro del borde del mapa.

Si no hay superposición, podría replicar una técnica que utilicé con éxito en el pasado (¡coincidentemente en E&W Sussex!) En MapInfo, donde escribí un pequeño script que generó un conjunto de cuatro puntos para cada entidad de índice , desplazamiento en las entidades adyacentes, con atributos tanto del número de hoja como de la dirección del desplazamiento. La capa de puntos se usó para generar etiquetas nuevamente, con la dirección de desplazamiento permitiendo que la orientación de las etiquetas se ajuste para un mejor efecto.

No lo he intentado, pero es posible que pueda evitar generar una capa de datos separada en QGIS mediante el uso de la nueva funcionalidad de diseño del generador de geometría, ¡esto sería una solución más elegante y dinámica que no se podía lograr en MapInfo!

En realidad, usted ya hizo la mayor parte del trabajo requerido para determinar los mosaicos que desea imprimir usando atlas. Pero el punto es cómo ajustar todo junto para mostrar solo las ID de mosaico que necesita. Para demostrar mi idea, usaré en este ejemplo una imagen DEM y un archivo vectorial de cuadrícula, como puede ver a continuación:

Primero necesitamos mostrar la etiqueta de cada cuadrícula.

En la vista de diseño, utilicé la cuadrícula como capa de cobertura en el atlas, creé dos mapas: el mapa de la ventana de vista principal y un mapa de índice que muestra solo la cuadrícula, como puede ver a continuación:

Luego hice lo siguiente:

1. Ajusté la escala del mapa de índice para mostrar la extensión de la cuadrícula completa y luego arreglé la escala
2. Arreglé la extensión de la vista para evitar que el mapa se desplace al usar Preview atlas, y
3. Permití Overviewver el alcance y la ubicación del mapa de la vista principal, como puede ver a continuación:

Para el mapa de la ventana de vista principal, arreglé la escala en la medida de cada bloque de cuadrícula, para asegurarme de que la escala no cambiará si sucede algo, como puede ver a continuación;

Con un mapa de índice, puede ver fácilmente la ID y la ubicación de cada mosaico con referencia a otro mosaico, incluso cuando apaga la cuadrícula desde la ventana del mapa de vista principal. Por ejemplo, el siguiente mapa tiene una ID de mosaico = 14, y puede ver las ID de mosaico circundantes.

Actualización :

Actualizaré mi respuesta porque me di cuenta de que deseaba mostrar el índice del número de página de los diseños circundantes, no los ID de los diseños circundantes.

Para facilitar la comprensión del proceso, actualizaré los números de identificación en el mapa de índice para mostrar el número de página de diseño, como se muestra a continuación:

Dado que las ID que tengo comienzan desde 0 (Cero), la ID de la primera cuadrícula que se muestra en el mapa de índice comenzará desde 3. Por lo tanto, quiero cambiar el número de página para comenzar desde 1 restando 2 del número de ID en Atlas: Page number: ID -2, luego usaré el número de página actual como referencia en la expresión para crear etiquetas para la página actual, la página anterior, la página siguiente, la página superior y la página inferior, de la siguiente manera:

• La página actual tiene esta expresión en el cuadro de texto de la etiqueta: Current Page Number: [%@atlas_pagename%]

• Expresión de página anterior: [%if((@atlas_pagename = 1), Null, '↑ Page Number: ' || (@atlas_pagename - 1))%]dado que no hay páginas anteriores a 1

• Expresión de página siguiente: [%if( (@atlas_pagename = 25), Null, '↓ Page Number: ' || (@atlas_pagename + 1))%]ya que no hay páginas después de 25

• Subir expresión de página: [%if((@atlas_pagename <= 6),NULL,'↑ Page Number: ' || (@atlas_pagename -6))%]ya que no hay páginas anteriores a 6 en la dirección superior

• Debajo de la expresión de la página: [%if((@atlas_pagename >= 20), Null, '↓ Page Number: ' || (@atlas_pagename + 6))%]ya que no hay páginas después de 20 en la dirección inferior

Algunos resultados de salida:

Esta solución funcionará para cuadrículas rectangulares y es automática (debería funcionar para cualquier escenario sin ajustar nada manualmente).

Supongamos que tiene una cuadrícula con números de página. Puede ejecutar mi secuencia de comandos de procesamiento seleccionando la capa de cuadrícula y su campo de número de página como parámetros. El script crea cuatro campos ( right, left, above, below) en la capa de cuadrícula y calcula la identificación de la página vecina correspondiente para cada celda de la cuadrícula. Luego puede usar sus expresiones (por ejemplo, [% if( "left" is not NULL, 'to page' || "left", "" ) %]) para mostrar las etiquetas de las páginas vecinas.

Simplemente agregue mi repositorio ( https://github.com/gacarrillor/QGIS-Resources.git ) desde el complemento QGIS Resource Sharing e instale el script:

Cómo funciona

La secuencia de comandos determina la relación (derecha, izquierda, arriba o abajo) comparando las coordenadas de los cuadros delimitadores tanto de la celda de la cuadrícula actual como de cada celda de intersección. Resulta que para cada relación, falta una de las coordenadas.

Si la relación es above, la coordenada que falta es yMin, es decir, todas las otras 3 coordenadas del cuadro delimitador de la celda de la cuadrícula actual estarán presentes en el cuadro delimitador de la celda anterior. Recuerde que los cuadros delimitadores de QGIS se definen en este orden: [xMin, yMin, xMax, yMax].

Para un ejemplo numérico, tomemos rectángulos con lados de longitud 1. Digamos que el cuadro delimitador de la celda actual se define como bbox1=[0,0,1,1]. El cuadro delimitador de la celda anterior se definiría como bbox2=[0,1,1,2]. Las coordenadas X de bbox1 están presentes en bbox2, mientras que las bbox1 yMinfaltan en las coordenadas Y de bbox2.

Podemos definir nuestras 4 relaciones de esta manera (o: presente, #: faltante):

right: [#,o,o,o]
above: [o,#,o,o]
left:  [o,o,#,o]
below: [o,o,o,#]

Como puede ver, el índice faltante nos brinda toda la información que necesitamos.