Python w analizach geoprzestrzennych

Python jest językiem ogólnego przeznaczenia a nie środowiskiem interaktywnej analizy danych. Wiele ułatwień obecnych w środowisku R nie jest dostępnych ze względu na inną naturę jęzka.Podobnie jak w innych współczesnych narzędziach programistycznych dostęp do danych geoprzestrzennych odbywa się przy pomocy bibliotek GDAL/OGR, które będą omawiane w czasie całego kursu.

Import podstawowych bibliotek

Do realizacji zadania potrzebujemy podstawowego zestawu bibliotek:

  • numpy - do przetwarzania danych numerycznych (generuje warning)
  • osgeo gdal i ogr - do wczytania i podstawowego przetwarzania danych geoprzestrzennych
  • os - do realizacji zadań systemowych
  • matplotlib do wizualizacji wyników

UWAGA - w podstawowej wersji GDAL/OGR działają tylko z Python 2.

In [1]:
import matplotlib.pyplot as plt
import numpy as np
from osgeo import gdal
from osgeo import ogr
import os

/usr/lib/python2.7/dist-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

Ścieżki dostępu

Ścieżki dostępu wymagają dostosowania do lokalizacji danych, dotyczy to zmiennej "path", gdzie powinna zostać podana aktualna ścieżka dostępu do katalogu zawierającego ćwiczenie

In [2]:
# sprawdzić scieżki dostępu
path = os.path.join(os.path.expanduser("~"),"Documents")
shapefile = os.path.join(path,"zadanie0","sample.shp")
rastfile = os.path.join(path,"zadanie0","green.tif")
print shapefile
print rastfile

/home/jarekj/Documents/zadanie0/sample.shp
/home/jarekj/Documents/zadanie0/green.tif

wczytanie danych rastrowych

Wczytadnie danych rastrowych wymaga wykonania trzech poleceń:

  • otwarcia dostępu do danych rastrowych (dane nie są wczytywane
  • wskazania warstwy (pasma) w danych rastrowych, które będzie przetwarzane, domyślnie jest to pasmo 1 (nie ma zerowego)
  • pobranie danych do transformacji innych obiektów. transformacje te zawierają odpowiednio:

| easting | ew-resolution | 0 | northing | 0 | ns-resolution |

gdal.Open nie wczytuje danych do pamięci, a jedynie tworzy połączenie ze zbiorem danych.

In [3]:
raster = gdal.Open(rastfile)
band = raster.GetRasterBand(1)
gtr = raster.GetGeoTransform() # paramtery transformacji rastra
print gtr

(345000.0, 250.0, 0.0, 520000.0, 0.0, -250.0)

Wczytanie danych wektorowych

Ta czynność jest nieco bardziej skomplikowana (dokładniej: warper na funkcję języka C jest napisany mniej inteligentnie niż w przypadku gdal). W pierwszej kolejności tworzymy obiekt driver (drv) tu stosując rozwiązanie polegające na wykorzystaniu nazwy sterownika dostępu do danych

W następnym kroku otwieramy dostęp do danych i pobieramy połączenie do warstwy danych. Dane nie są wczytywane ale twożony jest tzw. iterator, a więc obiekt posiadający wyłącznie metodę .next(), która odczytuje jednorazowo aktualny element listy, a następnie zamienia go następnym. Iterator nie może być w prosty sposób przewinięty do początku, za każdym razem trzeba na nowo odnawiać połączenie (czyli wywoływać ponownie poniższe dwa polecenia). Ostatnie polecenie tworzy nam pustą listę, do której będziemy zapisywać wartości odczytywane z rastra.

In [4]:
drv = ogr.GetDriverByName("ESRI Shapefile")
vector = drv.Open(shapefile,0)
layer = vector.GetLayer()
lista=[]

Pętla odczytująca dane

Dane odczytywane są punkt po punkcie w pętli for:

  • pobierane są współrzędne x i y danego punktu
  • wspórzędne są przeliczane na index komórki w rastrze poprzez zastosowanie geotransformacji (więcej później)
  • polecenie ReadAsArray pobiera fragment danych zaczynający się w pixelu o współrzędnych px i py oraz o wielkości 1 i 1. Zwraca dwuwymiarową tablicę.
  • każdy element jest wydobywany z dwuwymiarowej tablicy i dodawany do pustej listy

na końcu lista zostaje wyświetlona.

In [5]:
for feature in layer:
    geom = feature.GetGeometryRef()
    x,y  = geom.GetX(), geom.GetY()
    px = int((x - gtr[0]) / gtr[1]) #x pixel
    py = int((y - gtr[3]) / gtr[5]) #y pixel
    val=band.ReadAsArray(px,py,1,1)
    lista.append(val[0][0])  

print lista

[0.15695313, 0.47809625, 0.24658442, 0.55048352, 0.66184962, 0.41637826, 0.38210374, 0.55072433, 0.47989923, 0.33881006, 0.63541627, 0.32926854, 0.3315393, 0.60909092, 0.17488512, 0.33083308, 0.47297901, 0.57804716, 0.6551103, 0.56595707, 0.70661676, 0.67248076, 0.64350539, 0.70449948, 0.67657459, 0.24534434, 0.38752335, 0.58717322, 0.03556066, 0.47259578, 0.38739946, 0.43412685, 0.59262776, 0.28067645, 0.22196788, 0.49889642, 0.45273513, 0.36439869, 0.31357798, 0.47129419, 0.32363704, 0.58946741, 0.71876633, 0.36591971, 0.32014713, 0.38484544, 0.21629527, 0.16487837, 0.13388792, 0.2803137, 0.60606372, 0.44894746, 0.41901439, 0.22451718, 0.5031144, 0.32262817, 0.61455905, 0.38783967, 0.63581491, 0.17688844, 0.55593538, 0.52711761, 0.622199, 0.19539903, 0.41544437, 0.64590353, 0.64931899, 0.45413813, 0.51736796, 0.48016798, nan, 0.71468055]

Tworzenie wykresu

Ostatnim krokiem jest przygotowanie i wyświetlenie wykresu, do czego zostanie wykorzystana biblioteka matplotlib. W pierwszym kroku lista jest zamieniana na numpy array, a następnie wyświetlany jest boxplot. Boxplot nie toleruje wartości pustych oraz nieskończoności. Z tego powodu należy przefiltrować obiekt array, aby usunąć z niego wyszystkie elementy, które nie są skończoną wartością liczbową.

In [6]:
arr = np.array(lista)
plt.figure()
plt.boxplot(arr[np.isfinite(arr)], 0)
plt.show()