Como determinar os IDs de bloco vizinhos no QGIS?

Em um curso de treinamento recente, me perguntaram se o QGIS poderia calcular automaticamente os números de página seguinte / anterior e acima / abaixo de um livro de mapas criado usando o gerador de atlas. Consegui elaborar uma expressão de rótulo razoavelmente razoável para uma grade regular, se você souber a largura e a altura da grade.

Mas então começamos a pensar em exemplos realistas em que não queremos desenhar páginas que não contenham nosso distrito de interesse, como este em meu país natal:

Então, nesta tarde, eu joguei um script python para descobrir os 4 vizinhos nos quais eu estava interessado em cada célula da grade e adicionei esses valores à minha grade (isso é fortemente baseado no 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()

Isso funciona muito bem.

Mas, para ser sincero, criar um ponto de teste para o norte e testar todos os vizinhos possíveis parece errado. No entanto, depois de uma tarde destruindo meu cérebro, não consigo pensar em uma maneira melhor de determinar o que é o vizinho do norte de uma célula da grade em particular?

Idealmente, eu gostaria de algo simples o suficiente para colocar em uma caixa de texto do compositor de impressão, mas suspeito que seja pedir demais.

Se você não estiver ajustando cada extensão de página (da camada de índice) exatamente ao compositor, mas tiver bordas sobrepostas com páginas adjacentes (como mostrado na segunda captura de tela), poderá usar rótulos da camada de índice, com a desvantagem que eles estariam dentro da borda do mapa.

Se não houver sobreposição, você poderá replicar uma técnica que usei com sucesso no passado (coincidentemente na E&W Sussex!) No MapInfo, onde escrevi um pequeno script que gerava um conjunto de quatro pontos para cada recurso de índice , deslocamento para os recursos adjacentes, com atributos do número da folha e a direção do deslocamento. A camada de pontos foi então usada para gerar rótulos novamente, com a direção do deslocamento, permitindo que a orientação dos rótulos fosse ajustada para um efeito melhor.

Eu não tentei isso, mas você pode evitar gerar uma camada de dados separada no QGIS através do uso da nova funcionalidade de estilo de gerador de geometria, isso criaria uma solução mais elegante e dinâmica que não era possível no MapInfo!

Na verdade, você já realizou a maior parte do trabalho necessário para determinar os blocos que deseja imprimir usando o atlas. Mas o ponto é como ajustar tudo para mostrar apenas os IDs de bloco que você precisa. Para demonstrar minha ideia, usarei neste exemplo uma imagem DEM e um arquivo vetorial de grade, como você pode ver abaixo:

Primeiro, precisamos mostrar o rótulo de cada grade.

Na visualização de layout, usei grade como camada de cobertura no atlas, criei dois mapas: o mapa da janela principal e um mapa de índice que mostra apenas a grade, como você pode ver abaixo:

Então eu fiz o seguinte:

1. Ajustei a escala do mapa de índice para mostrar toda a extensão da grade e depois fixei a escala
2. Corrigi a extensão da visualização para impedir que o mapa se movesse ao usar Preview atlase
3. Eu ativei o Overviewpara ver a extensão e a localização do mapa principal, como você pode ver abaixo:

Para o mapa da janela principal, eu fixei a escala na extensão de cada bloco da grade, para garantir que a escala não será alterada se acontecer alguma coisa, como você pode ver abaixo;

Usando um mapa de índice, é possível ver facilmente o ID e o local de cada bloco com referência a outro bloco, mesmo quando você desativa a grade na janela principal do mapa da visualização. Por exemplo, o mapa a seguir tem um ID de bloco = 14 e você pode ver os IDs de bloco ao redor.

Atualização :

Atualizarei minha resposta porque percebi que você queria mostrar o índice de número de página dos layouts circundantes, não os IDs dos layouts circundantes.

Para facilitar o entendimento do processo, atualizarei os números de identificação no mapa de índice para mostrar o número da página de layout, conforme mostrado abaixo:

Como os IDs iniciados em 0 (Zero), o ID da primeira grade mostrada no mapa de índice começará em 3. Portanto, quero alterar o número da página para começar em 1 subtraindo 2 do número de ID no Atlas: Page number: ID -2, usarei o número da página atual como referência na expressão para criar rótulos para a página atual, página anterior, próxima página, página acima e página abaixo, da seguinte maneira:

• A página atual tem esta expressão na caixa de texto do rótulo: Current Page Number: [%@atlas_pagename%]

• Expressão de página anterior: [%if((@atlas_pagename = 1), Null, '↑ Page Number: ' || (@atlas_pagename - 1))%]como não há páginas antes de 1

• Expressão da próxima página: [%if( (@atlas_pagename = 25), Null, '↓ Page Number: ' || (@atlas_pagename + 1))%]como não há páginas após 25

• Up Page expression: [%if((@atlas_pagename <= 6),NULL,'↑ Page Number: ' || (@atlas_pagename -6))%]como não existem páginas antes das 6 na direção superior

• Abaixo Expressão da página: [%if((@atlas_pagename >= 20), Null, '↓ Page Number: ' || (@atlas_pagename + 6))%]como não há páginas após 20 na direção inferior

Alguns resultados de saída:

Esta solução funcionará para grades retangulares e é automática (deve funcionar em qualquer cenário sem ajustar nada manualmente).

Vamos supor que você tenha uma grade com números de página. Você pode executar meu script de processamento selecionando a camada de grade e seu campo de número de página como parâmetros. O script cria quatro campos ( right, left, above, below) na camada de grade e calcula o ID da página vizinha correspondente para cada célula da grade. Então você pode usar suas expressões (por exemplo, [% if( "left" is not NULL, 'to page' || "left", "" ) %]) para mostrar rótulos de páginas vizinhas.

Basta adicionar meu repositório ( https://github.com/gacarrillor/QGIS-Resources.git ) no plug-in QGIS Resource Sharing e instalar o script:

Como funciona

O script determina a relação (direita, esquerda, acima ou abaixo) comparando as coordenadas das caixas delimitadoras da célula atual da grade e de cada célula que se cruza. Acontece que, para cada relação, falta uma das coordenadas.

Se a relação for above, a coordenada ausente será yMin, ou seja, todas as outras 3 coordenadas da caixa delimitadora da célula da grade atual estarão presentes na caixa delimitadora da célula acima. Lembre-se que caixas delimitadoras QGIS são definidos na seguinte ordem: [xMin, yMin, xMax, yMax].

Para um exemplo numérico, vamos usar retângulos com lados de comprimento 1. Digamos que a caixa delimitadora da célula atual esteja definida como bbox1=[0,0,1,1]. A caixa delimitadora da célula acima seria definida como bbox2=[0,1,1,2]. As coordenadas X do bbox1 estão presentes no bbox2, enquanto as do bbox1 yMinestão ausentes nas coordenadas Y do bbox2.

Podemos definir nossas 4 relações dessa maneira (o: presente, #: ausente):

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

Como você pode ver, o índice ausente nos fornece todas as informações que precisamos.