Rasterio: Numpy for Geospatial Data

Cześć, tu Alex z DEV STORIES DOT EU. Rasterio: Numpy dla danych geoprzestrzennych, odcinek 1 z 5. Chcesz wyciągnąć wartości pikseli ze zdjęcia satelitarnego, ale standardowe biblioteki do obrazów nie radzą sobie z formatami geoprzestrzennymi, a oficjalne bindingi przypominają pisanie kodu w C wewnątrz Pythona. Właśnie przez ten zgrzyt powstało Rasterio. Rasterio to biblioteka, która czyta i zapisuje dane geoprzestrzenne w formie siatki, takie jak pliki GeoTIFF. Bardzo możliwe, że słyszałeś o GDAL, bibliotece geoprzestrzennej będącej standardem w branży. Powszechnym mitem jest to, że Rasterio zastępuje GDAL. Wcale tak nie jest. Rasterio tak naprawdę używa GDAL pod spodem. Różnica tkwi w interfejsie. Jeśli kiedykolwiek korzystałeś z natywnych bindingów GDAL w Pythonie, wiesz, że wystawiają one toporne C API, pełne ręcznego zarządzania zasobami i składni, która wydaje się zupełnie obca w Pythonie. Rasterio zostawia to w tyle. Zapewnia nowoczesny interfejs w podejściu Python-first. Traktuje datasety rastrowe jak normalne pliki, a ich piksele obsługuje jak standardowe Numpy arrays. Załóżmy, że masz na dysku plik GeoTIFF ze zdjęciami Landsat. Otwierasz ten plik za pomocą funkcji rasterio kropka open. Ponieważ Rasterio trzyma się standardowych wzorców Pythona, robisz to używając with statement. To tworzy context manager. Podobnie jak przy otwieraniu zwykłego pliku tekstowego, blok with gwarantuje, że dataset zostanie czysto zamknięty, a pamięć zwolniona w momencie zakończenia działania twojego kodu. Nigdy nie musisz ręcznie niszczyć obiektów ani martwić się o memory leaks z niezamkniętych plików. Wewnątrz tego bloku with, funkcja open daje ci obiekt dataset reader. Przed odczytaniem jakichkolwiek właściwych pikseli, możesz zbadać plik. Możesz zapytać dataset o jego właściwość count, która mówi ci o całkowitej liczbie pasm, czyli warstw, na zdjęciu. Dla zdjęcia Landsat, ten count może wynosić coś koło siedmiu lub jedenastu, reprezentując różne długości fal światła. Możesz odczytać właściwości width i height, aby uzyskać wymiary przestrzenne w pikselach. Możesz również sprawdzić właściwość dtypes. To daje ci data type pikseli dla każdego pasma, na przykład unsigned 16-bit integers lub 32-bit floats. Wszystkie te metadane są dostępne natychmiast, bez ładowania ciężkich danych obrazu do pamięci systemowej. Oto kluczowa kwestia. Kiedy jesteś gotowy, by spojrzeć na właściwe piksele, wywołujesz metodę read na obiekcie dataset. Możesz przekazać konkretny indeks, aby załadować tylko jedną warstwę. Zwróć uwagę na ten szczegół: pasma w Rasterio są one-indexed, co oznacza, że prosisz o pasmo pierwsze, a nie zerowe. Jeśli wywołasz read i przekażesz jedynkę, metoda zwróci te piksele jako standardowy, dwuwymiarowy Numpy array. Jeśli wywołasz read bez żadnych argumentów, odczyta wszystko, zwracając trójwymiarowy Numpy array zawierający wszystkie pasma. Nie ma żadnego specjalnego, własnościowego pixel object. Dostajesz po prostu array liczb. Kiedy już masz ten Numpy array, wracasz na znane terytorium. Możesz go slice'ować, odpalać operacje statystyczne, albo przekazać bezpośrednio do modeli machine learning, używając dokładnie tego samego kodu, co dla każdej innej macierzy. Rasterio działa jak most. Obsługuje skomplikowane formaty plików i metadane na dysku, i przekazuje ci czysty obiekt matematyczny w pamięci. Jeśli podoba ci się podcast i chcesz wesprzeć to, co robimy, możesz wyszukać DevStoriesEU na Patreon — to ogromna pomoc. Prawdziwa siła Rasterio nie polega na tym, że wymyśla na nowo przetwarzanie geoprzestrzenne, ale na tym, że czyni je nudnym; zamieniając skomplikowane zdjęcia satelitarne w standardowe Numpy arrays, pozwala ci przestać walczyć z formatami GIS i zacząć robić prawdziwe data science. Dzięki za wysłuchanie, miłego kodowania wszystkim!

Cześć, tu Alex z DEV STORIES DOT EU. Rasterio: Numpy dla danych geoprzestrzennych, odcinek 2 z 5. Podstawowy Numpy array to po prostu płaska siatka liczb. Jeśli to zwizualizujesz, na monitorze pojawi się po prostu kolorowy prostokąt. Co tak naprawdę zmienia te abstrakcyjne liczby w precyzyjną mapę Płaskowyżu Kolorado? Odpowiedzią jest transformacja afiniczna. Ludzie często zakładają, że dane geoprzestrzenne to po prostu obraz, którego rogi są luźno przypięte do współrzędnych mapy. Wcale tak nie jest. Transformacja to ścisłe powiązanie matematyczne, które oblicza dokładne położenie każdego pojedynczego piksela w twoim arrayu w świecie rzeczywistym. Zanim zmapujesz piksele na Ziemię, musisz wiedzieć, jak ją mierzysz. To jest Coordinate Reference System. Odczytujesz to z datasetu, odwołując się do dataset dot crs. To property informuje Rasterio, czy twoje współrzędne przestrzenne są mierzone w stopniach długości i szerokości geograficznej, czy w metrach względem określonego równika i południka zerowego. Gdy system odniesienia jest już ustawiony, potrzebujesz logiki mapowania. Znajdziesz ją, sprawdzając dataset dot transform. Zwraca to macierz transformacji afinicznej. Brzmi to jak ciężka algebra liniowa, ale w praktyce to zestaw sześciu liczb, które opisują fizyczny ślad twojego grida. Macierz definiuje dokładne współrzędne przestrzenne x i y lewego górnego rogu lewego górnego piksela. Definiuje fizyczną szerokość jednego piksela, na przykład dokładnie trzydzieści metrów. Definiuje fizyczną wysokość jednego piksela, która zazwyczaj jest ujemna, ponieważ wiersze obrazu liczone są w dół, podczas gdy współrzędne mapy rosną w górę. Macierz przechowuje również wartości rotacji na wypadek, gdyby satelita był pochylony w momencie robienia zdjęcia. Oto kluczowa sprawa. Nigdy nie musisz mnożyć tych wartości macierzy ręcznie. Rasterio daje ci bezpośrednie metody, żeby skakać tam i z powrotem między przestrzenią pikseli a przestrzenią współrzędnych. Załóżmy, że analizujesz obraz i znajdujesz anomalię w konkretnej lokalizacji, powiedzmy wiersz pięćsetny i kolumna tysięczna. Aby dowiedzieć się, gdzie to jest na prawdziwej planecie, przekazujesz te dwa integery do metody dataset dot xy. Rasterio przepuszcza wiersz i kolumnę przez macierz afiniczną i zwraca dokładne współrzędne przestrzenne x i y dla środkowego punktu tego konkretnego piksela. Możesz też uruchomić tę maszynę w drugą stronę. Powiedzmy, że masz dokładne współrzędne przestrzenne dla fizycznej lokalizacji. Może to kanion położony dokładnie sto kilometrów na wschód i pięćdziesiąt kilometrów na południe od originu twojego obrazu. Musisz wyciągnąć dane pikseli dla tego kanionu. Bierzesz swoje przestrzenne wartości x i y i przekazujesz je do metody dataset dot index. Rasterio odwraca macierz afiniczną, przetwarza twoje współrzędne i zwraca dokładne integery wiersza i kolumny. Następnie używasz tych integerów, żeby zrobić slice na twoim Numpy arrayu i odczytać dokładny pomiar zarejestrowany przez satelitę. Transformacja afiniczna izoluje twój array surowych danych od geometrii przestrzennej, gwarantując, że niezależnie od tego, jak zrobisz crop lub odczytasz swoje dane, nigdy nie stracisz swojej precyzyjnej pozycji na kuli ziemskiej. To wszystko na dziś. Dzięki za wysłuchanie — idź zbudować coś fajnego.

Cześć, tu Alex z DEV STORIES DOT EU. Rasterio: numpy dla danych geoprzestrzennych, odcinek 3 z 5. Jeśli spróbujesz wczytać dziesięciogigabajtowy plik GeoTIFF bezpośrednio do pamięci laptopa, twój skrypt prawie na pewno się wywali. Możesz potrzebować tylko małego wycinka pikseli do przetestowania algorytmu, ale system próbuje połknąć cały kontynent naraz. Rozwiązaniem tego problemu z pamięcią jest Windowed Processing. Windowed Processing pozwala ci czytać lub zapisywać konkretne, prostokątne wycinki pliku rastrowego, podczas gdy reszta pliku leży bezpiecznie na dysku twardym. Opiera się to całkowicie na indeksach. Nie zajmujemy się dzisiaj współrzędnymi geograficznymi. Operujemy wyłącznie na wierszach i kolumnach pikseli. Aby wczytać taki wycinek, musisz powiedzieć Rasterio, które piksele ma pobrać. Robisz to, importując obiekt Window z modułu rasterio windows. Standardowym sposobem na zdefiniowanie obiektu Window jest podanie czterech liczb: column offset, row offset, width oraz height. Zwróć na to uwagę. Kolejność argumentów ma znaczenie. Najpierw podajesz column offset, a potem row offset. Column offset to twoja początkowa pozycja w poziomie, mierzona od lewej krawędzi. Row offset to początkowa pozycja w pionie, mierzona w dół od górnej krawędzi. Po offsetach określasz, ile pikseli szerokości i wysokości ma mieć twoje zaznaczenie. Wyobraź sobie ogromny, gigabajtowy plik testowy. Chcesz wyciągnąć z niego mały kafelek o szerokości 256 i wysokości 512 pikseli. Chcesz, żeby ten kafelek zaczynał się tysiąc pikseli od lewej i dwa tysiące pikseli od góry. Tworzysz instancję obiektu Window, przekazując tysiąc jako column offset, dwa tysiące jako row offset, 256 jako width i 512 jako height. Następnie otwierasz swój dataset za pomocą standardowej funkcji open. Kiedy wywołujesz metodę read, nie zostawiasz pustych argumentów. Puste argumenty mówią Rasterio, żeby załadował cały plik do pamięci. Zamiast tego, przypisujesz swój nowy obiekt Window do keyword argumentu window wewnątrz metody read. Rasterio tłumaczy to okno na byte offsets na dysku, pobiera tylko ten konkretny chunk 256 na 512 pikseli i zwraca go jako numpy array. Twoje zużycie pamięci ledwie drgnie. Jeśli na co dzień pracujesz z numpy, przekazywanie offsetów i wymiarów może wydawać się trochę nienaturalne. Prawdopodobnie jesteś przyzwyczajony do definiowania slice'ów za pomocą indeksów start i stop dla wierszy i kolumn. Rasterio wspiera ten workflow za pomocą alternatywnej metody o nazwie Window dot from slices. Metoda ta przyjmuje dwa argumenty. Pierwszy z nich to tupla definiująca indeksy start i stop dla wiersza. Drugi to tupla definiująca indeksy start i stop dla kolumny. Zauważ, że kolejność to najpierw wiersze, a potem kolumny. Odzwierciedla to dokładnie sposób, w jaki robisz slice'y na dwuwymiarowym numpy array. Pod spodem ta metoda oblicza za ciebie offsety i wymiary, po czym zwraca standardowy obiekt Window, który możesz przekazać do metody read dokładnie tak, jak wcześniej. Wymuszając windowed reads w swoich data pipelines, na stałe uniezależniasz fizyczny rozmiar plików wejściowych od limitów sprzętowych twojej maszyny. To wszystko w tym odcinku. Dzięki za wysłuchanie i koduj dalej!

Cześć, tu Alex z DEV STORIES DOT EU. Rasterio: Numpy dla danych geoprzestrzennych, odcinek 4 z 5. Właśnie pobrałeś ogromny obraz satelitarny obejmujący cały stan, ale interesują cię tylko piksele wewnątrz granic jednej, konkretnej farmy. Ładowanie całej siatki marnuje pamięć, a odczyt na podstawie offsetów pikseli jest bezużyteczny, gdy twoja granica to nieregularny polygon. Rozwiązaniem jest maskowanie wektorami. W poprzednim odcinku przyjrzeliśmy się windowed reading. To podejście wymagało dokładnych indeksów pikseli, na przykład zapytania o pięćdziesiąty wiersz i setną kolumnę. Maskowanie działa zupełnie inaczej. Maskowanie opiera się na rzeczywistych współrzędnych geograficznych. Zamiast liczyć piksele, podajesz geometrię wektorową, na przykład polygon zdefiniowany przez rzeczywistą długość i szerokość geograficzną. Rasterio samo określa, które piksele wpadają do tego kształtu i zajmuje się złożoną translacją między ciągłą przestrzenią współrzędnych wektorowych a dyskretną siatką rastrową. Głównym narzędziem do tego workflow jest funkcja mask, znajdująca się w module rasterio dot mask. Przekazujesz do niej dwa główne argumenty. Po pierwsze, otwarty dataset rastrowy. Po drugie, iterable kształtów geometrycznych. Te kształty muszą być sformatowane jako słowniki zgodne ze specyfikacją GeoJSON. Zazwyczaj wczytujesz shapefile zawierający polygon twojej farmy, wyciągasz słownik geometrii dla tego konkretnego feature'a, wrzucasz go do standardowej listy w Pythonie i przekazujesz do funkcji mask. I tu jest kluczowa sprawa. Funkcja mask ma dwa różne zachowania sterowane pojedynczym parametrem typu boolean o nazwie crop. Musisz zrozumieć różnicę między zwykłym maskowaniem datasetu a jego faktycznym cropowaniem. Jeśli odpalisz funkcję z parametrem crop ustawionym na false, Rasterio sprawdzi polygon twojej farmy i ustawi każdy piksel poza tą granicą na wartość nodata. Domyślną wartość nodata odczytuje bezpośrednio z oryginalnych metadanych rastra. Wyjściowy array, który otrzymasz, ma dokładnie takie same wymiary jak twój ogromny, oryginalny obraz całego stanu. Piksele wewnątrz farmy zachowują swoje rzeczywiste wartości, ale wszystko na zewnątrz jest wyzerowane. Całkowity zasięg przestrzenny pliku nie ulega zmianie. Jeśli ustawisz crop na true, zachowanie się zmienia. Rasterio nadal ustawia piksele poza twoim polygonem na nodata. Jednak potem odrzuca te ogromne obszary pustej przestrzeni. Oblicza prostokątny bounding box twojego nieregularnego polygonu farmy i pomniejsza zwracany array tak, aby idealnie wpasował się w ten konkretny box. Twój gigantyczny array całego stanu zostaje fizycznie zredukowany do małej, łatwej w zarządzaniu siatki, zawierającej tylko farmę i niewielki padding z wartości nodata wokół jej nieregularnych krawędzi. Ponieważ cropowanie fizycznie zmienia wymiary siatki, zmienia również geograficzny punkt początkowy. Lewy górny róg nie jest już lewym górnym rogiem stanu. Jest teraz lewym górnym rogiem bounding boxa farmy. Dlatego funkcja mask zwraca tuplę zawierającą dwa elementy. Pierwszy element to twój nowo wycropowany numpy array. Drugi element to zupełnie nowy obiekt affine transform. Ten zaktualizowany transform przechowuje nowe współrzędne origin. Jeśli chcesz zapisać ten nowo wycropowany array do nowego pliku, musisz przekazać ten zaktualizowany transform do swojego write profile. Jeśli przez pomyłkę użyjesz oryginalnego transformu dla całego stanu, obraz twojej małej farmy zostanie rozciągnięty z powrotem na cały stan. Prawdziwą zaletą funkcji mask jest to, że automatycznie ogarnia za ciebie całe to przestrzenne przeliczanie. Mam nadzieję, że to było przydatne. Dzięki za wysłuchanie i miłego dnia!

Cześć, tu Alex z DEV STORIES DOT EU. Rasterio: Numpy dla danych geoprzestrzennych, odcinek 5 z 5. Czasami masz piękne zdjęcia satelitarne o wysokiej rozdzielczości, ale twój model klimatyczny akceptuje tylko zgrubne dane wejściowe o niskiej rozdzielczości. Nie możesz po prostu zrobić slice'a arraya, ponieważ wynikowe piksele nadal muszą reprezentować dokładnie ten sam fizyczny obszar, tylko w siatce o mniejszej gęstości. Aby zrobić to poprawnie, potrzebujesz resamplingu i zapisu. Nie chodzi tu o reprojekcję. Układ współrzędnych pozostaje dokładnie taki sam. Zmieniamy tylko rozdzielczość. Załóżmy, że musisz zdownsamplować raster o połowę. W Rasterio resampling wykonujesz dokładnie w momencie wczytywania danych do pamięci. Wywołując metodę read na swoim datassecie, przekazujesz parametr o nazwie out_shape. Zamiast wczytywać całego arraya, Rasterio wczytuje go i natychmiast skaluje, aby dopasować do żądanego shape'a. W przypadku rastra jednopasmowego, twój out_shape byłby tuplą określającą jedno pasmo, a następnie połowę pierwotnej wysokości i połowę pierwotnej szerokości. Musisz również powiedzieć Rasterio, jak obliczyć wartości tych nowych, większych pikseli. Robisz to, przekazując parametr resampling. Rasterio udostępnia w tym celu enuma Resampling. Jeśli masz do czynienia z danymi ciągłymi, takimi jak temperatura, możesz przekazać Resampling dot bilinear, który oblicza średnią ważoną otaczających pikseli. Jeśli masz do czynienia z danymi kategorialnymi, takimi jak pokrycie terenu, użyłbyś Resampling dot nearest, aby zachować dokładnie te same, oryginalne wartości klas. Oto kluczowa sprawa. Masz teraz nowego, mniejszego arraya Numpy, ale nadal reprezentuje on dokładnie ten sam obszar geograficzny na Ziemi. Ponieważ array ma mniej pikseli, każdy pojedynczy piksel pokrywa większą przestrzeń fizyczną. Oznacza to, że twój oryginalny transform Affine jest teraz całkowicie błędny. Jeśli zapiszesz nowego arraya na dysku przy użyciu starego transformu, twój obraz na mapie skurczy się, pokrywając tylko jedną czwartą pierwotnego obszaru. Musisz przeskalować transform. Robisz to, biorąc oryginalny transform datasetu i mnożąc go przez macierz skalującą Affine. Współczynniki skalowania znajdujesz, dzieląc oryginalną szerokość przez twoją nową szerokość i oryginalną wysokość przez twoją nową wysokość. Ponieważ zmniejszyliśmy wymiary o połowę, oba te współczynniki wynoszą dokładnie dwa. Kiedy pomnożysz oryginalny transform przez te współczynniki skalowania, otrzymasz nowy transform, w którym wymiary pikseli są podwojone, ale punkt początkowy w lewym górnym rogu pozostaje idealnie zakotwiczony. Teraz masz swojego zresamplowanego arraya i przeskalowany transform. Ostatnim krokiem jest zapisanie poprawnego pliku GeoTIFF na dysku. Aby to zrobić, potrzebujesz metadanych. Najbezpieczniejszym podejściem jest skopiowanie właściwości profile z twojego źródłowego datasetu. To da ci słownik zawierający driver pliku, układ współrzędnych, typ danych i wartość nodata. Następnie aktualizujesz ten słownik profile o twoją nową wysokość, nową szerokość i nowo obliczony transform. Mając gotowy zaktualizowany profile, wywołujesz rasterio dot open. Podajesz nową nazwę pliku, ustawiasz mode na zapis i rozpakowujesz swój słownik profile jako keyword arguments. Tworzy to pusty dataset na dysku, sformatowany dokładnie według twoich specyfikacji. Na koniec wywołujesz metodę write na tym nowym obiekcie datasetu i przekazujesz swojego zdownsamplowanego arraya Numpy. Array jest zapisywany do pliku, a twój nowy GeoTIFF o niższej rozdzielczości jest gotowy dla modelu klimatycznego. Pamiętaj, modyfikowanie danych rastrowych to zawsze dwuczęściowy kontrakt: jeśli zmienisz shape arraya Numpy, musisz jawnie zmienić transform Affine, aby go dopasować, w przeciwnym razie twój ślad przestrzenny będzie zepsuty. To już koniec naszej serii. Zachęcam cię do przeczytania oficjalnej dokumentacji Rasterio, aby poznać inne algorytmy resamplingu i spróbować zbudować te pipeline'y w praktyce. Jeśli chcesz zaproponować tematy do przyszłej serii, odwiedź devstories dot eu. Chciałbym poświęcić chwilę, aby podziękować ci za słuchanie — to bardzo nam pomaga. Miłego dnia!