Jak określić sąsiadujące identyfikatory kafelków w QGIS?

Podczas ostatniego szkolenia zapytano mnie, czy QGIS może automatycznie obliczyć numery następnych / poprzednich i powyżej / poniżej stron książek z mapami utworzonych za pomocą generatora atlasu. Udało mi się wypracować dość rozsądne wyrażenie etykiety dla regularnej siatki, jeśli znasz szerokość i wysokość siatki.

Ale potem zaczęliśmy myśleć o realistycznych przykładach, w których nie chcemy rysować stron, które nie zawierają naszej dzielnicy zainteresowań, takich jak ten w moim rodzinnym hrabstwie:

Więc tego popołudnia bawiłem się skryptem Pythona, aby obliczyć 4 sąsiadów, którymi byłem zainteresowany dla każdej komórki siatki i dodałem te wartości do mojej siatki (jest to w dużej mierze oparte na tutorialu Ujavala Gandhiego ):

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()

To działa dobrze.

Ale szczerze mówiąc, stworzenie punktu testowego na północ, a następnie przetestowanie wszystkich możliwych sąsiadów wydaje się błędne. Jednak po popołudniowym zniszczeniu mózgu nie mogę wymyślić lepszego sposobu na określenie, kto jest północnym sąsiadem konkretnej komórki sieci?

Idealnie chciałbym, aby coś prostego umieścić w polu tekstowym kompozytora wydruku, ale podejrzewam, że to zbyt wiele, o co mogę prosić.

Jeśli nie dopasowujesz każdego zasięgu strony (z warstwy indeksu) dokładnie do kompozytora, ale zamiast tego masz nakładające się granice z sąsiednimi stronami (jak pokazano na drugim zrzucie ekranu), możesz użyć etykiet z warstwy indeksu, z wadą byłyby wewnątrz granicy mapy.

Jeśli nie ma żadnego nakładania się, możesz powtórzyć technikę, którą z powodzeniem stosowałem w przeszłości (przypadkowo w E & W Sussex!) W MapInfo, gdzie napisałem mały skrypt, który wygenerował zestaw czterech punktów dla każdej funkcji indeksu , przesunąć do sąsiednich elementów, z atrybutami zarówno numeru arkusza, jak i kierunku przesunięcia. Warstwę punktową wykorzystano następnie do ponownego wygenerowania etykiet, przy czym kierunek przesunięcia umożliwia dostosowanie orientacji etykiet w celu uzyskania lepszego efektu.

Nie próbowałem tego, ale możesz uniknąć generowania osobnej warstwy danych w QGIS dzięki nowej funkcji generowania stylów generatora geometrii, dzięki czemu byłoby bardziej eleganckie i dynamiczne rozwiązanie, którego nie można było uzyskać w MapInfo!

W rzeczywistości wykonałeś już większość pracy wymaganej do określenia kafelków, które chcesz wydrukować za pomocą atlasu. Chodzi o to, jak dopasować wszystko razem, aby wyświetlać tylko potrzebne identyfikatory kafelków. Aby zademonstrować mój pomysł, użyję w tym przykładzie obrazu DEM i pliku wektorowego siatki, jak widać poniżej:

Najpierw musimy pokazać etykietę każdej siatki.

W widoku układu użyłem siatki jako warstwy pokrycia w atlasie, utworzyłem dwie mapy: mapę okna głównego widoku i mapę indeksu, która pokazuje tylko siatkę, jak widać poniżej:

Następnie wykonałem następujące czynności:

1. Skorygowałem skalę mapy indeksu, aby pokazać cały zasięg siatki, a następnie poprawiłem skalę
2. Poprawiłem zasięg widoku, aby zapobiec panoramowaniu mapy podczas korzystania Preview atlas, oraz
3. Włączyłem opcję, Overviewaby zobaczyć zasięg i lokalizację głównej mapy widoku, jak widać poniżej:

W przypadku mapy okna głównego widoku ustawiłem skalę w zakresie każdego bloku siatki, aby mieć pewność, że skala nie zostanie zmieniona, jeśli cokolwiek się stanie, jak widać poniżej;

Korzystając z mapy indeksu, możesz łatwo zobaczyć identyfikator i lokalizację każdego kafelka w odniesieniu do innego kafelka, nawet po wyłączeniu siatki z głównego okna mapy widoku. Na przykład poniższa mapa ma identyfikator kafelka = 14 i można zobaczyć otaczające identyfikatory kafelków.

Aktualizacja :

Zaktualizuję moją odpowiedź, ponieważ zdałem sobie sprawę, że chcesz pokazać indeks numerów stron otaczających układów, a nie identyfikatory otaczających układów.

Aby ułatwić zrozumienie tego procesu, zaktualizuję numery identyfikacyjne na mapie indeksu, aby pokazać numer strony układu, jak pokazano poniżej:

Ponieważ identyfikatory, które mam, zaczynają się od 0 (zero), identyfikator pierwszej siatki pokazany na mapie indeksu zacznie się od 3. Dlatego chcę zmienić numer strony, aby zaczynał się od 1, odejmując 2 od numeru identyfikacyjnego w Atlasie: Page number: ID -2, użyję bieżącego numeru strony jako odwołania w wyrażeniu, aby utworzyć etykiety dla bieżącej strony, poprzedniej strony, następnej strony, strony w górę i poniżej strony w następujący sposób:

• Bieżąca strona ma to wyrażenie w polu tekstowym etykiety: Current Page Number: [%@atlas_pagename%]

• Wyrażenie poprzedniej strony: [%if((@atlas_pagename = 1), Null, '↑ Page Number: ' || (@atlas_pagename - 1))%]ponieważ przed 1 nie ma żadnych stron

• Wyrażenie następnej strony: [%if( (@atlas_pagename = 25), Null, '↓ Page Number: ' || (@atlas_pagename + 1))%]ponieważ po 25 nie ma stron

• Wyrażenie w górę strony: [%if((@atlas_pagename <= 6),NULL,'↑ Page Number: ' || (@atlas_pagename -6))%]ponieważ nie ma stron przed 6 w górnym kierunku

• Wyrażenie poniżej strony: [%if((@atlas_pagename >= 20), Null, '↓ Page Number: ' || (@atlas_pagename + 6))%]ponieważ po 20 nie ma stron w dolnym kierunku

Niektóre wyniki wyjściowe:

To rozwiązanie będzie działać dla prostokątnych siatek i jest automatyczne (powinno działać w każdym scenariuszu bez ręcznego dostosowywania).

Załóżmy, że masz siatkę z numerami stron. Możesz uruchomić mój skrypt Przetwarzanie, wybierając warstwę siatki i pole numeru strony jako parametry. Skrypt tworzy cztery pola ( right, left, above, below) w warstwie siatki i oblicza odpowiedni identyfikator strony sąsiada dla każdej komórki siatki. Następnie możesz użyć wyrażeń (np. [% if( "left" is not NULL, 'to page' || "left", "" ) %]), Aby wyświetlić etykiety stron sąsiadów.

Wystarczy dodać moje repozytorium ( https://github.com/gacarrillor/QGIS-Resources.git ) z wtyczki QGIS Resource Sharing i zainstalować skrypt:

Jak to działa

Skrypt określa relację (prawą, lewą, górną lub dolną) poprzez porównanie współrzędnych obwiedni zarówno z bieżącej komórki siatki, jak i każdej przecinającej się komórki. Okazuje się, że dla każdej relacji brakuje jednej ze współrzędnych.

Jeśli relacja jest above, brakująca współrzędna jest yMin, tzn. Wszystkie pozostałe 3 współrzędne z bieżącego obwiedni komórki siatki będą obecne w obwiedni powyższej komórki. Pamiętaj, że QGIS ograniczające pola są zdefiniowane w następującej kolejności: [xMin, yMin, xMax, yMax].

Na przykład liczbowy weźmy prostokąty o bokach długości 1. Powiedzmy, że obwiednia bieżącej komórki jest zdefiniowana jako bbox1=[0,0,1,1]. Obwiednia powyższej komórki zostałaby zdefiniowana jako bbox2=[0,1,1,2]. Współrzędne X z bbox1 są obecne w bbox2, podczas gdy brak jest bbox1 yMinwe współrzędnych Y bbox2.

W ten sposób możemy zdefiniować nasze 4 relacje (o: present, #: missing):

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

Jak widać, brakujący indeks daje nam wszystkie potrzebne informacje.