Comment déterminer les ID de tuile voisins dans QGIS?

11

Lors d'un récent cours de formation, on m'a demandé si QGIS pouvait calculer automatiquement les numéros de page suivant / précédent et supérieur / inférieur pour un atlas créé à l'aide du générateur d'atlas. J'ai réussi à trouver une expression d'étiquette assez raisonnable pour une grille régulière si vous connaissez la largeur et la hauteur de la grille.

Mais nous avons ensuite commencé à penser à des exemples réalistes où nous ne voulons pas dessiner des pages qui ne contiennent pas notre district d'intérêt, comme celui de mon comté d'origine:

entrez la description de l'image ici

Cet après-midi, j'ai donc joué à un script python pour déterminer les 4 voisins qui m'intéressaient pour chaque cellule de la grille et j'ai ajouté ces valeurs à ma grille (cela est fortement basé sur le tutoriel d' 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()

Cela fonctionne très bien.

entrez la description de l'image ici

Mais pour être honnête, le tout créant un point de test vers le Nord puis testant tous les voisins possibles semble faux. Cependant, après un après-midi à me casser la tête, je ne peux pas penser à une meilleure façon de déterminer quel est le voisin du nord d'une cellule de grille particulière?

Idéalement, j'aimerais quelque chose d'assez simple à mettre dans une zone de texte de compositeur d'impression, mais je soupçonne que c'est trop demander.

Ian Turton
la source
et s'il n'y a pas de voisin d'un côté. Voulez-vous la valeur de la cellule close dans une direction ou voulez-vous laisser un vide?
radouxju
Je suis heureux pour une valeur nulle dans ce cas, je peux facilement définir l'étiquette pour qu'elle s'affiche uniquement lorsqu'elle n'est pas nulle ou vide.
Ian Turton

Réponses:

3

Si vous n'ajustez pas chaque étendue de page (à partir de la couche d'index) exactement dans le compositeur, mais que vous avez plutôt des bordures qui se chevauchent avec des pages adjacentes (comme indiqué dans votre deuxième capture d'écran), vous pouvez utiliser des étiquettes de la couche d'index, avec l'inconvénient que ils seraient à l'intérieur de la frontière de la carte.

S'il n'y a pas de chevauchement, vous pouvez reproduire une technique que j'ai utilisée avec succès dans le passé (par coïncidence dans E & W Sussex!) Dans MapInfo, où j'ai écrit un petit script qui a généré un ensemble de quatre points pour chaque entité d'index , décalé dans les entités adjacentes, avec des attributs à la fois le numéro de feuille et la direction du décalage. La couche ponctuelle a ensuite été utilisée pour générer à nouveau des étiquettes, la direction du décalage permettant d'ajuster l'orientation des étiquettes pour un effet plus agréable.

Je n'ai pas essayé cela, mais vous pourriez éviter de générer une couche de données distincte dans QGIS en utilisant la nouvelle fonctionnalité de style du générateur de géométrie, cela ferait une solution plus élégante et dynamique qui n'était pas réalisable dans MapInfo!

Andy Harfoot
la source
J'aurais vraiment dû penser à utiliser simplement les étiquettes des autres polygones! :-) Après une rapide expérience avec le générateur de géométrie, je peux dessiner un cadre de délimitation, mais il est plus difficile de construire une grille
Ian Turton
Je pensais à générer des points d'étiquette décalés dans les polygones adjacents, plutôt que des grilles. Une autre option serait d'étendre le MBR de l'entité d'indexation aux entités adjacentes pour permettre de dessiner des étiquettes.
Andy Harfoot
Je viens de jouer et il semble que la géométrie générée par le style du générateur de géométrie ne soit pas étiquetée, ce n'est donc pas la solution la plus élégante que j'espérais.
Andy Harfoot
8

En fait, vous avez déjà effectué la plupart des travaux nécessaires pour déterminer les tuiles que vous souhaitez imprimer à l'aide d'atlas. Mais le fait est de savoir comment tout régler ensemble pour n'afficher que les ID de mosaïque dont vous avez besoin. Pour démontrer mon idée, je vais utiliser dans cet exemple une image DEM et un fichier vectoriel de grille, comme vous pouvez le voir ci-dessous:

entrez la description de l'image ici

Nous devons d'abord montrer l'étiquette de chaque grille.

Dans la vue de mise en page, j'ai utilisé la grille comme couche de couverture dans l'atlas, j'ai créé deux cartes: la carte de la fenêtre de vue principale et une carte d'index qui montre uniquement la grille, comme vous pouvez le voir ci-dessous:

entrez la description de l'image ici

Ensuite, j'ai fait ce qui suit:

  1. J'ai ajusté l'échelle de la carte d'index pour afficher toute l'étendue de la grille, puis j'ai fixé l'échelle
  2. J'ai corrigé l'étendue de la vue pour empêcher la carte de se déplacer lors de l'utilisation Preview atlas, et
  3. J'ai activé le Overviewpour voir l'étendue et l'emplacement de la carte de la vue principale, comme vous pouvez le voir ci-dessous:

entrez la description de l'image ici

Pour la carte de la fenêtre de vue principale, j'ai fixé l'échelle à l'étendue de chaque bloc de grille, pour être sûr que l'échelle ne sera pas modifiée si quelque chose se produit, comme vous pouvez le voir ci-dessous;

entrez la description de l'image ici

À l'aide d'une carte d'index, vous pouvez facilement voir l'ID et l'emplacement de chaque tuile en référence à une autre tuile, même lorsque vous désactivez la grille à partir de la fenêtre de carte de la vue principale. Par exemple, la carte suivante a un ID de tuile = 14 et vous pouvez voir les ID de tuile environnants.

entrez la description de l'image ici

Mise à jour :

Je mettrai à jour ma réponse car j'ai réalisé que vous vouliez afficher l'index des numéros de page des mises en page environnantes et non les ID des mises en page environnantes.

Pour faciliter la compréhension du processus, je mettrai à jour les numéros d'identification dans la carte d'index pour afficher le numéro de page de mise en page, comme indiqué ci-dessous:

entrez la description de l'image ici

Étant donné que les ID que je commence à partir de 0 (zéro), l'ID de la première grille affichée sur la carte d'index commencera à partir de 3. Par conséquent, je veux changer le numéro de page pour commencer à 1 en soustrayant 2 du numéro ID dans Atlas: Page number: ID -2, je vais utiliser le numéro de page actuel comme référence dans l'expression pour créer des étiquettes pour la page actuelle, la page précédente, la page suivante, la page précédente et la page inférieure, comme suit:

entrez la description de l'image ici

  • La page actuelle a cette expression dans la zone de texte de l'étiquette: Current Page Number: [%@atlas_pagename%]

  • Expression de la page précédente: [%if((@atlas_pagename = 1), Null, '↑ Page Number: ' || (@atlas_pagename - 1))%]puisqu'il n'y a pas de pages avant 1

  • Expression de la page suivante: [%if( (@atlas_pagename = 25), Null, '↓ Page Number: ' || (@atlas_pagename + 1))%]puisqu'il n'y a pas de pages après 25

  • Expression de la [%if((@atlas_pagename <= 6),NULL,'↑ Page Number: ' || (@atlas_pagename -6))%]page précédente : puisqu'il n'y a pas de pages avant 6 dans la direction supérieure

  • Au-dessous de l'expression de la page: [%if((@atlas_pagename >= 20), Null, '↓ Page Number: ' || (@atlas_pagename + 6))%]puisqu'il n'y a pas de pages après 20 dans la direction inférieure

Quelques résultats de sortie:

entrez la description de l'image ici

entrez la description de l'image ici

entrez la description de l'image ici

entrez la description de l'image ici

ahmadhanb
la source
Bien qu'utile, cela ne répond pas à sa question.
Victor
@Victor Merci pour votre commentaire, j'ai mis à jour ma réponse.
ahmadhanb
cela fonctionne dans votre exemple (et le sien), car les côtés de la carte / grille sont réguliers. S'ils n'étaient pas droits, cela ne fonctionnerait pas, car le nombre à ajouter ou à soustraire (6 dans votre exemple) varierait en fonction de la page d'atlas sur laquelle vous vous trouvez.
Victor
2
Je suis d'accord avec toi. Si la grille n'est pas régulière, le processus sera plus compliqué. Mais comme il veut l'appliquer sur une grille régulière, la méthode appliquée dans ma solution suggérée fonctionnera.
ahmadhanb
juste noter le fait, au cas où vous auriez une autre bonne idée! D'autant plus que ma grille n'est pas régulière!
Victor
2

Cette solution fonctionnera pour les grilles rectangulaires et elle est automatique (devrait fonctionner pour n'importe quel scénario sans rien ajuster manuellement).

Supposons que vous ayez une grille avec des numéros de page. Vous pouvez exécuter mon script de traitement en sélectionnant la couche de grille et son champ de numéro de page comme paramètres. Le script crée quatre champs ( right, left, above, below) dans la couche de grille et calcule l'ID de page voisin correspondant pour chaque cellule de grille. Ensuite, vous pouvez utiliser vos expressions (par exemple, [% if( "left" is not NULL, 'to page' || "left", "" ) %]) pour afficher les étiquettes des pages voisines.

Ajoutez simplement mon référentiel ( https://github.com/gacarrillor/QGIS-Resources.git ) à partir du plugin QGIS Resource Sharing et installez le script: entrez la description de l'image ici

entrez la description de l'image ici

Comment ça fonctionne

Le script détermine la relation (droite, gauche, au-dessus ou en dessous) en comparant les coordonnées des boîtes englobantes à la fois de la cellule de grille actuelle et de chaque cellule qui se croisent. Il s'avère que pour chaque relation, une des coordonnées est manquante.

Si la relation est above, la coordonnée manquante est yMin, c'est-à-dire que toutes les 3 autres coordonnées du cadre de délimitation de la cellule de grille actuelle seront présentes dans le cadre de délimitation de la cellule ci-dessus. Rappelez - vous que les boîtes englobantes QGIS sont définies dans cet ordre: [xMin, yMin, xMax, yMax].

Pour un exemple numérique, prenons des rectangles avec des côtés de longueur 1. Supposons que le cadre de sélection de la cellule actuelle soit défini comme bbox1=[0,0,1,1]. La boîte englobante de la cellule ci-dessus serait définie comme bbox2=[0,1,1,2]. Les coordonnées X de bbox1 sont présentes dans bbox2, alors que bbox1 yMinmanque dans les coordonnées Y de bbox2.

Nous pouvons définir nos 4 relations de cette façon (o: présent, #: manquant):

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

Comme vous pouvez le voir, l'index manquant nous donne toutes les informations dont nous avons besoin.

Germán Carrillo
la source