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


11

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:

wprowadź opis zdjęcia tutaj

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.

wprowadź opis zdjęcia tutaj

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ć.


co jeśli nie ma sąsiada po jednej stronie. Czy chcesz wartość zamkniętej komórki w jednym kierunku, czy pozostawiłbyś pustkę?
radouxju

W tym przypadku cieszę się z pustej wartości, mogę łatwo ustawić etykietę, aby wyświetlała się tylko wtedy, gdy nie jest pusta lub pusta.
Ian Turton

Odpowiedzi:


3

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!


Naprawdę powinienem pomyśleć o użyciu etykiet innych wielokątów! :-) Po szybkim eksperymencie z Generatorem Geometrii mogę narysować ramkę ograniczającą, ale trudniej jest zbudować siatkę
Ian Turton

Myślałam wzdłuż linii generowania punktów etykiety przesuniętych do sąsiednich wielokątów, a nie siatek. Inną opcją byłoby rozszerzenie MBR elementu indeksu na sąsiednie elementy, aby umożliwić rysowanie etykiet.
Andy Harfoot,

Właśnie się pobawiłem i wygląda na to, że geometria wygenerowana przez styl generatora geometrii nie jest oznaczana, więc nie jest to bardziej eleganckie rozwiązanie, na jakie liczyłem.
Andy Harfoot,

8

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:

wprowadź opis zdjęcia tutaj

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:

wprowadź opis zdjęcia tutaj

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:

wprowadź opis zdjęcia tutaj

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;

wprowadź opis zdjęcia tutaj

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.

wprowadź opis zdjęcia tutaj

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:

wprowadź opis zdjęcia tutaj

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:

wprowadź opis zdjęcia tutaj

  • 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:

wprowadź opis zdjęcia tutaj

wprowadź opis zdjęcia tutaj

wprowadź opis zdjęcia tutaj

wprowadź opis zdjęcia tutaj


Chociaż jest to przydatne, nie odpowiada na jego pytanie.
Victor,

@Victor Dziękujemy za komentarz, zaktualizowałem swoją odpowiedź.
ahmadhanb,

działa to w twoim przykładzie (i jego), ponieważ boki mapy klawiszy / siatki są regularne. Gdyby nie były proste, nie działałoby, ponieważ liczba do dodania lub odjęcia (6 w twoim przykładzie) będzie się różnić w zależności od strony atlasu, na której jesteś.
Victor,

2
Zgadzam się z Tobą. Jeśli siatka nie jest regularna, proces będzie bardziej skomplikowany. Ale ponieważ chce zastosować go na regularnej siatce, metoda zastosowana w moim sugerowanym rozwiązaniu zadziała.
ahmadhanb,

tylko zauważając fakt, na wypadek, gdybyś miał inny dobry pomysł! Zwłaszcza, że ​​moja siatka nie jest regularna!
Victor,

2

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: wprowadź opis zdjęcia tutaj

wprowadź opis zdjęcia tutaj

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.

Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.