OPIS Polskiej Siatki Binarnej PolBiG

Celem utworzenia opisanej w artykule (link poniżej) binarnej siatki PolBiG (patrz: szczegółowy widok etykiet pól 64×64 i siatki 8×8 km, na tle granic województw i gmin polskich) było umożliwienie jej współpracy z programami rodziny GIS, precyzyjne i proste określanie powierzchni pól w dowolnie ustalonym zagęszczeniu dla pełnego zasięgu siatki, naturalny sposób kodowania (etykietowania) pól binarnych, łatwość przenoszenia się pomiędzy układami (powiązanym z Ziemią i płaskim), wreszcie ułatwione implementowanie siatki w przyszłych użytkowych aplikacjach „różniczkujących” powierzchniowo.

Artykuł na temat siatki - Polska Siatka Binarna - PolBiG
Patrz też: ATPOL - polska siatka geobotaniczna
----------------------------

Siatka PolBiG o wymiarach 1024×1024 km jest rozwinięta na płaszczyźnie stycznej do kuli (o promieniu R) w punkcie dowiązania (Xo, Yo)=(512, 512), a zastosowane odwzorowanie oparte jest o wzory matematyczne wyprowadzone z modelu azymutalnej transformacji równopowierzchniowej Lamberta (Lambert Azimuthal Equal Area Projection). Przekształcenie odwrotne f -1 określa z definicji wartości współrzędnych geograficznych w układzie odniesienia WGS 84 przenosząc kwadratową siatkę określoną w prostokątnym układzie współrzędnych (odpowiednikami kwadratów są pola określane tak po transformacji na elipsoidę WGS 84):
$$\begin{array}{l}\\Dla\;przekształcenia\;(x,\;y)=f\;(\lambda,\;\varphi):\\x=Xo+R\times p\times\cos\;\varphi\times\sin(\lambda-\beta)\\y=Yo+R\times p\times(\mathrm{co}s\;\alpha\times\sin\;\varphi-\sin\;\alpha\times\cos\;\varphi\times\cos(\lambda-\beta))\\\\Dla\;przekształcenia\;odwrotnego\;(\lambda,\;\varphi)=f^{-1}(x,\;y):\\\varphi=arc\;\sin(\cos\;c\times\sin\;\alpha+(y-Yo)\times\sin\;c\times\frac{\cos\;\alpha}g)\\\lambda=\beta+arc\;tg((x-Xo)\times\frac{\sin\;c}{g\;\times\;\cos\;c\;\times\;\cos\;\alpha\;-\;(y-Yo)\;\times\;\sin\;c\;\times\;\sin\;\alpha})\\\\gdzie:\\\boldsymbol p=\sqrt{2/(1+\sin\;\alpha\times\sin\;\varphi+\cos\;\varphi\times\cos\;\alpha\times\cos(\lambda-\beta))}\\\boldsymbol g=\sqrt{(x-Xo)^2+(y-Yo)^2}\;\;;\;\;\boldsymbol c=2\times arc\;\sin(\frac g{2\times R})\\\end{array}$$
Dla zainteresowanych udostępniam prosty skoroszyt MS Excel gdzie w dwóch kolejnych arkuszach przepisałem formalnie powyższe wzory transformacyjne (przeliczenia wewnątrz funkcji trygonometrycznych wykonujemy stosując miarę łukową kąta, w przypadku funkcji odwrotnych powracamy do miary stopniowej). Łatwo więc sprawdzić zgodność transformacyjną w obie strony, zgodność punktu dowiązania (Xo,Yo) = (512,512) dla współrzędnych geograficznych α = 52°N i β = 19°E a także wyróżnioną własność południka β = 19°E, któremu zawsze odpowiada wartość współrzędnej kartezjańskiej x = 512. Skoroszyt pozwala też zamieniać współrzędne geograficzne na współrzędne kartezjańskie siatki i odwrotnie w ilościach hurtowych, wystarczy wstawić je do odpowiedniego arkusza i przeliczać na zasadzie kopiuj - wklej.

ETYKIETOWANIE SIATKI PolBiG

Binarna siatka PolBiG pozwala na proste wprowadzenie etykietowania "oczek" siatki w dowolnym jej zagęszczeniu. Duże kwadraty siatki PolBiG (64×64 km) można oznaczyć małymi literami w kolejności wiersz-kolumna (podobnie jak w zapisie macierzowym) jedynie dla lepszej orientacji w położeniu geograficznym dużych pól binarnych. Ogólny schemat etykietowania opisany jest w cytowanej powyżej pracy - opiera się na numerowaniu cyframi 1,2,3,4 kolejnych ćwiartek kwadratów "zagęszczanych" w ten binarny sposób. Numeracja dopisywana jest w kolejnym podziale do etykiety kwadratu siatki po jej prawej stronie. Duże kwadraty binarne 64×64 km po trzykrotnym podziale dzielą się już na 64 "oczka" o wymiarach 8×8 km. Także kwadraty 8×8 km po następnym trzykrotnym podziale utworzą 64 mniejszych "oczek" o wymiarze 1×1 km. Dla obu tych przypadków zapis etykiet jest identyczny. Należy tylko pamiętać o dodaniu do etykiety bardziej gęstego podziału odpowiadającą mu etykietę podziału poprzedniego. Poza zapisem literowym dla "oczek" 8×8 km będą to 3 cyfry, dla 1×1 km będzie ich w sumie 6.

Przykład etykietowania kwadratów 8×8 km i 1×1 km:
  • Podglądamy mapę Polski z siatką PolBiG i opisem etykiet (wersja literowa lub cyfrowa)
  • Pobieramy uniwersalny szablon etykiet dla trzykrotnego podziału binarnego kwadratów (np z dużego kwadratu binarnego na kwadraty 8×8 km, lub z kwadratów 8×8 km na kw. 1×1 km)
  • Do etykiety odczytanej z mapy dopisujemy z jej prawej strony dalszy ciąg etykiety odczytany z szablonu. Przykładowo, Hel umiejscowiony jest w kwadracie #1244 (lub inaczej "dh"), a dokładniej w podkwadracie 243 (np 18,805217°E i 54,603947°N). Można więc oznaczyć kwadrat 8×8 km dla Helu jako #1244243 lub bardziej skrótowo dh243.
  • Aby uzyskać etykiety 1×1 km postępujemy analogicznie. Wystarczy zobaczyć wygląd siatki kilometrowej PolBiG na tle ośmiokilometrowej: popatrz ..., większa część miasteczka Hel znajduje się w podkwadracie 342 co łatwo odczytać z udostępnionego szablonu. Mamy więc w końcu dla Helu etykietę 1×1 km: #1244243342 lub krócej: dh243342.

KALKULATOR PolBiG:
Kalkulator binarny PolBiG
Kalkulator można pobrać w formacie XLSM (Excel) zawierającym makra (podprogramy w języku VBA), stąd konieczność akceptacji ich uruchomienia po pobraniu. Kalkulator przelicza wartości współrzędnych geograficznych na odpowiadające im etykiety (także w wersji literowo-cyfrowej) po wyborze wielkości kwadratu siatki. Współrzędne geograficzne muszą określać miejsca zawierające się w maksymalnej rozpiętości siatki PolBiG. W innym przypadku wynik obliczeń będzie zapewne błędny. Z drugiej strony - odwrotnie, kalkulator podaje współrzędne geograficzne wierzchołków i środka kwadratu wyliczając je w oparciu o symbol wpisanej etykiety. Przystosowany jest do obliczeń dla siatki binarnej: Pobierz kalkulator

SIATKI WZORCOWE w formacie KMZ:
Wzorcowa siatka (1024 × 1024 km) tworzona z kwadratów 8×8 km, w formacie KMZ (PolBiG) wraz z etykietami w wersji literowo-numerycznej i numerycznej 8×8 km do pobrania: Literowo-numeryczna... Numeryczna z #
Wzorcowa siatka (512 × 512 km: 1) tworzona z kwadratów 4×4 km, w formacie KMZ (PolBiG) wraz z etykietami w wersji literowo-numerycznej i numerycznej 4×4 km do pobrania: Literowo-numeryczna... Numeryczna z # plus odpowiadająca poligonowa 2×2 km.
Wzorcowa siatka (512 × 512 km: 2) tworzona z kwadratów 4×4 km, w formacie KMZ (PolBiG) wraz z etykietami w wersji literowo-numerycznej i numerycznej 4×4 km do pobrania: Literowo-numeryczna... Numeryczna z # plus odpowiadająca poligonowa 2×2 km.
Wzorcowa siatka (512 × 512 km: 3) tworzona z kwadratów 4×4 km, w formacie KMZ (PolBiG) wraz z etykietami w wersji literowo-numerycznej i numerycznej 4×4 km do pobrania: Literowo-numeryczna... Numeryczna z # plus odpowiadająca poligonowa 2×2 km.
Wzorcowa siatka (512 × 512 km: 4) tworzona z kwadratów 4×4 km, w formacie KMZ (PolBiG) wraz z etykietami w wersji literowo-numerycznej i numerycznej 4×4 km do pobrania: Literowo-numeryczna... Numeryczna z # plus odpowiadająca poligonowa 2×2 km.
  • UWAGA - Siatki tworzone były z linii korygowanych odpowiednio co 8 i co 4 km. Siatka 1024 × 1024 km o oczkach 4×4 podzielona została na 4 ćwiartki zgodnie z ustaloną konwencją aby ułatwić wczytywanie do G.E. Do tych ostatnich siatek dołączono odpowiadające im siatki 2×2 km w wersji poligonowej (kolor czerwony). Nie mają one co prawda etykiet, wystarczy jednak dołożyć cyfrę 1, 2, 3 lub 4 wg. ustalonej już konwencji do odpowiadającej etykiety 4×4 km aby zaetykietować odpowiednio kwadrat 2×2 km.
  • UWAGA - Można wczytać siatki 4×4 km i 2×2 km do G.E., w pierwszej kolejności wczytując siatkę 2×2 km. Ułatwi to rozróżnienie kolorystyczne konkretnych "oczek" i ich etykietowanie. Należy jednak uzbroić się w cierpliwość przy "zamykaniu" G.E. Może to trwać nawet kilka minut. Nie jest jasne co jest tego powodem, same siatki wczytują się stosunkowo szybko.
  • UWAGA - Wszelkie siatki o innych wielkościach i gęstościach (od boku o długości od 64 km do 1/1024 km), wraz z etykietami (zaetykietowano niepuste "oczka") w wersji numerycznej pobierzemy dynamicznie z pomocą udostępnionego tu programu. W siatkach wzorcowych wszystkie "oczka" siatki 8×8 km i 4×4 km zostały zaetykietowane.

SIATKI WZORCOWE w formacie Shapefile (spakowane zip):
Wzorcowa siatka (1024 × 1024 km) tworzona z kwadratów 8×8 km, w formacie Shapefile (PolBiG) wraz z etykietami w wersji literowo-numerycznej i numerycznej 8×8 km do pobrania: Literowo-numeryczna... Numeryczna z #
Wzorcowa siatka (1024 × 1024 km) tworzona z kwadratów 4×4 km, w formacie Shapefile (PolBiG) wraz z etykietami w wersji literowo-numerycznej i numerycznej 4×4 km do pobrania: Literowo-numeryczna... Numeryczna z #
Wzorcowa siatka (1024 × 1024 km) tworzona z kwadratów 2×2 km, w formacie Shapefile (PolBiG) wraz z etykietami w wersji literowo-numerycznej i numerycznej 4×4 km do pobrania: Literowo-numeryczna... Numeryczna z #
  • UWAGA - Ostatnia siatka ma co prawda etykiety dla kwadratów 4×4 km, wystarczy jednak dołożyć cyfrę 1, 2, 3 lub 4 wg. ustalonej już konwencji do odpowiadającej etykiety 4×4 km aby zaetykietować odpowiednio kwadrat 2×2 km.
  • UWAGA - Wszystkie siatki (wersja poligonowa) i etykiety w formacie SHP w oparciu o powyżej udostępnione pliki w formacie KMZ opracował Piotr Kobierski.
PRZYKŁADOWE SIATKI PolBiG (poligonowe, generowane z QGIS):

Binarna (1024 × 1024 km) - kwadrat 8×8 km
Binarna (1024 × 1024 km) - kwadrat 4×4 km
Binarna (część Małopolski) - kwadrat 1×1 km
Binarna (Kraków-Borek Fałęcki) - kwadrat 0,25×0,25 km

Dziesiętna (1000 × 1000 km) - kwadrat 10×10 km
Dziesiętna (1000 × 1000 km) - kwadrat 5×5 km

Jak wygenerować siatkę PolBiG z pomocą oprogramowania QGIS ?
  1. Instalujemy QGIS w aktualnie dostępnej wersji (opis dotyczy w. 3.2 dla Windows).
  2. Przechodzimy do Ustawień i dalej do Układu współrzędnych użytkownika.
  3. Tworzymy definicję układu współrzędnych przyciskając + w kolorze zielonym.
  4. Wprowadzamy nazwę układu np. Polbig1024  , aby wygenerować siatkę o maksymalnej wielkości 1024 km.
  5. Do okna Parametry wstawiamy następujący tekst (definicję układu zgodną z Proj4- https://proj.org):
    +proj=laea +lat_1=52 +lat_2=52 +lat_0=52 +lon_0=19 +a=6371000 +b=6371000 +ellps=sphere +x_0=512000 +y_0=512000

    lub też krótszym alternatywnym zapisie, zgodnym z Proj4 i z definicją "własnego" ukł. wsp. w pakiecie QGIS:
    +proj=laea +lat_0=52 +lon_0=19 +R=6371000 +x_0=512000 +y_0=512000
  6. Potwierdzamy naszą definicję utworzonego układu współrzędnych przyciskiem OK i klikamy w EPSG:4326, czyli bieżący układ współrzędnych WGS 84 - przycisk znajduje się w prawym dolnym narożu okna programu QGIS.
  7. W uruchomionym oknie, w "Dostępnych układach współrzędnych", na samym dole (trzeba użyć beleczki) rozwijamy  > "Układy współrzędnych użytkownika" i wybieramy nasz zdefiniowany układ przyciskiem Zastosuj a następnie OK.
  8. Teraz możemy w tym zdefiniowanym przez nas układzie wygenerować siatkę. Przechodzimy do menu głównego, wybieramy Wektor , dalej Narzędzia badawcze i Utwórz siatkę.
  9. W powstałym oknie (Utwórz Siatkę) wybieramy kolejno: Prostokąt (poligon), zasięg siatki (xmin, xmax, ymin, ymax). W naszym przypadku, kiedy chcemy wygenerować pełną siatkę PolBiG, będzie to wpis: 0,1024000,0,1024000 bo dane wpisujemy w metrach. Można oczywiście podać inne wartości, nie zmieni to w żaden sposób samej siatki a jedynie jej wielkość. Ważne jest jednak by podać taki zakres, aby wymiar "oczek siatki" mieścił się bez reszty w podanym zakresie.
  10. Dobieramy teraz wymiar "oczek" czyli Odległość w poziomie i Odległość w pionie. Zgodnie z powyższym mogą to być przykładowe wielkości 16 km, 8 km, 4 km, 2 km, 1 km, 1/2 km, 1/4 km, 1/8 km  itd. Te wielkości podajemy w metrach, przykładowo wpisujemy 8000 dla obu pól. Uwaga, mały wymiar "oczek" generuje później duże pliki, bo każde "oczko (w ramach zadeklarowanego zakresu) zostanie opisane. Wtedy zakres należy ustawić na mniejszy.
  11. Resztę pozostawiamy bez zmian, sprawdzając, by zaznaczone było okienko: Wczytaj plik po zakończeniu.
  12. Wciskamy przycisk Uruchomienia w tle i dalej OK. Powinna ukazać się wygenerowana siatka w idealnej formie.
  13. W pasku po lewej stronie głównego okna w tzw Warstwach ukaże się wiersz zatytułowany Siatka. Możemy teraz prawym przyciskiem myszki uruchomić menu kontekstowe Siatki i po rozwinięciu wybrać Eksportuj. Dalej wybieramy Save feature as (Zapisz warstwę jako, ... u mnie ukazuje się w języku angielskim).
  14. W powstałym oknie dobieramy Format ... przykładowo Keyhole Markup Language [KML] (pozwala na zobaczenie siatki w popularnym programie G. Earth), lub też format ESRI Shapefile (należy wtedy wskazać pusty folder i wprowadzić własną nazwę pliku - w folderze wygeneruje się zestaw plików powiązanych tą nazwą). W Układzie współrzędnych, poniżej, należy wybrać taki układ współrzędnych jakim jesteśmy zainteresowani, a więc domyślny - EPSG: 4326 - WGS 84. To zagwarantuje, że siatka "położona na mapie" rysowanej w tym układzie współrzędnych będzie wyglądała prawidłowo.
  15. Poniżej "odznaczamy" opcję Dodaj zapisany plik do mapy, jest to ważne o tyle, że w ten sposób nadal pozostaniemy w naszym własnym "kartezjańskim" układzie, mimo, że sam plik (pliki) wygenerują się już w ogólnoświatowym standardzie WGS 84. Będziemy więc mogli nadal generować inne siatki binarne pozostając nadal w układzie współrzędnych prostokątnych.
  16. Klikamy OK , po chwili u góry głównego okna programu powinien ukazać się komunikat w zielonym kolorze z inf. o zapisanym pliku.

Osoby, które chciałyby generować inne siatki, np. 10×10 km oczywiście mogą to robić podobnie, jednak należy użyć innej definicji własnego układu współrzędnych. Przykładowo, jeżeli cała siatka ma mieć 1000×1000 km (nieco mniejsza od binarnej), to dla takiej wielkości siatki naturalne będą oczka o boku 10 km, 5 km, 2 km czy 1 km. Wtedy własny układ definiujemy tak:
+proj=laea +lat_1=52 +lat_2=52 +lat_0=52 +lon_0=19 +a=6371000 +b=6371000 +ellps=sphere +x_0=500000 +y_0=500000
lub:  +proj=laea +lat_0=52 +lon_0=19 +R=6371000 +x_0=500000 +y_0=500000

Jedyne co uległo zmianie to wartość liczbowa x_0 i y_0 (siatka ma mniejszy zasięg liczony z centralnego punktu), cała reszta parametrów pozostała ta sama, co gwarantuje wygenerowanie odpowiedniej siatki (to samo umiejscowienie geograficzne punktu dowiązania i ta sama wartość parametru R).
Warto zauważyć, że siatki PolBiG: binarna i dziesiętna, są identyczne dla kilometrowych i dwukilometrowych  "oczek", jednak sposób ich etykietowania byłby inny.

Marek Verey, Kraków, 29 czerwiec 2020