In [1]:
# setup
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
import warnings
warnings.filterwarnings('ignore')
In [2]:
import matplotlib.pyplot as plt
from matplotlib import colors
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
import numpy as np
import os

Geoprocessing

Przykład 1: Prezentacja danych rastrowych w naturalnych przedziałach

Celem zadania jest podział warstwy danych rastrowych na przedziały (percentyle), wyświetlenie percentyli w postaci wykresów pudełkowych oraz mapy kategoryzowanej gdzie wartości zostaną przypisane do percentyli. Pierwszy krok to zmiana katalogu roboczego oraz otwarcei warstwy rastrowej. Warstwa rastrowa posiada nodata o określonej wartości (minimalna wartość dla danego typu danych).

In [3]:
os.chdir('/home/jarekj/Dropbox/dydaktyka/programowanie2/data_raster')
In [4]:
raster = gdal.Open('pop_low.tif')
band = raster.GetRasterBand(1)
data = band.ReadAsArray()
nodata = band.GetNoDataValue()
nodata
Out[4]:
-1.7e+308

W celu wyświetlenia mapy z uwzględnieniem nodata można zamaskować wartości symbolizujące nodata, przy pomocy funkcji np.ma(sk) oraz wyświetlenie obrazu z dodaniem palety kolorów.

In [5]:
masked_data = np.ma.masked_equal(data,nodata)
plt.imshow(masked_data,cmap="rainbow")
plt.colorbar()
Out[5]:
<matplotlib.image.AxesImage at 0x7fe2278441d0>
Out[5]:
<matplotlib.colorbar.Colorbar at 0x7fe227862c88>

Do wykonania wykresów słupkowych użyjemy jedynie podzbioru danych. W pierwszej kolejności musimy pobrać tę część danych, która posiada wartości (nie null).

In [6]:
subset = masked_data[~masked_data.mask] # następuje spłaszczenie

Uzyjemy podzbioru 10% i 4 przedziałów.

In [7]:
sampling = 0.1 #
nclass = 4
breaks = np.arange(0,1,step=1./nclass) #2
  1. Liczb próbek to 10% całości danych bez wartości null, ponieważ musi to być liczba całkowita, wartość musi być zaokrąglona i zamieniona na liczbę całkowitą
  2. Następnie wybieramy ze zbioru indeksy próbek
  3. Wyliczamy kwantyle
  4. Następnie dokonujemy kategoryzacji danych na podstawie przedziałów
  5. Dzielimy cały zbiór na listę wartości dla poszczególnych kategorii, na podstawie listy unikalnych kategorii. Lista jest potrzebna do zbudowania serii wykresów pudełkowych
In [8]:
num_of_samples = int(round(len(subset)*sampling)) #1 musi być liczba całkowita
samples = np.random.choice(subset,num_of_samples,replace=False) #2
quants = np.quantile(subset,breaks) #3
quants
categories = np.digitize(samples,quants) #4
samplist = [samples[categories==cat] for cat in np.unique(categories)] #5
Out[8]:
array([ 0.        ,  0.35285249, 15.68345237, 62.47363497])

Wyświetlanie wykresów pudełkowych pozwala zdefiniować paletę kolorów. Niestety modyfikacja poszczególnych pudełek wymaga zastosowania specyficznej pętli.

In [9]:
cols = ["#FFA118","#FF111F","#AFA118","#9AABCD"]
boxes = plt.boxplot(samplist,vert=True,patch_artist=True)
for patch, color in zip(boxes['boxes'],cols):
    patch.set_facecolor(color)

Ostatnim krokiem jest przygotowanie mapy. Przy pomocy funkcji digitize() ponownie kategoryzujemy dane, a następnie maskujemy wartości zerowe, które pelnią fukcję wartości NULL. Wyświetlanie kolorów zgodnych z wykresem, musimy przygotować colormap(). Taką mapę można następnie wyświetlić w formie obrazu.

In [10]:
categories = np.digitize(masked_data,quants)
masked_categories = np.ma.masked_equal(categories,0)  

cmap = colors.ListedColormap(cols)
plt.imshow(masked_categories,cmap=cmap)
plt.colorbar()
Out[10]:
<matplotlib.image.AxesImage at 0x7fe17540f9e8>
Out[10]:
<matplotlib.colorbar.Colorbar at 0x7fe1753ae4a8>

Zadanie:

Kwantyle z próby użyte do podziału mapy i wykonania wykresów odbiegają znacząco od kwantyli z populacji, co powoduje że wykresy nie są do końca wiarygodne. Jak poprawić mechanizm próbkowania aby otrzymać bardziej wiarygodne przedziały (stratyfikacja). Wykorzystać informację o rozkładzie i pojęcie wartości odstających. Zaproponować inne przedziały oparte o liczności danych.

Przykład 2: Zamiana danych z formatu rastrowego w wektorowy na podstawie progowania wartości

Procesy konwersji jednego formatu na drugi, jest standardową procedurą geo-przetwarzania. W ramach procedury zmodyfikujemy raster zmiennoprzecinkowy, który zostanie zamieniony na wartości binarne, na podstawie jednej wartości.

In [11]:
raster = gdal.Open('pop_low.tif')
band = raster.GetRasterBand(1)
data = band.ReadAsArray()
nodata = band.GetNoDataValue()
In [12]:
x = raster.RasterXSize
y = raster.RasterYSize
proj = raster.GetProjection()
gt = raster.GetGeoTransform()

W celu zbudowania

In [13]:
masked_data = np.ma.masked_equal(data,nodata)
new_data = np.where(masked_data>2,1,0)

W celu dokonania dygitalizacji, tworzymy nowy raster w pamięci operacyjnej komputera używamy sterownika MEM zamiast standardowy geotiff. Tworzony raster jest typu Byte (8 bitów), do którego są dodawane wyniki procesu zmiany danych. Ustawiamy 0 aby symbolizował nodata.

In [14]:
driver = gdal.GetDriverByName('MEM')
new_rast= driver.Create("",x,y,1, gdal.GDT_Byte) 
new_rast.SetProjection(proj)
new_rast.SetGeoTransform(gt)

tband=new_rast.GetRasterBand(1)
tband.WriteArray(new_data)
tband.SetNoDataValue(0)
Out[14]:
0
Out[14]:
0
Out[14]:
0
Out[14]:
0

Proces dygitalizacji wymaga utworzenia nowego wektora. Aby utworzyć wektor w tym samym układzie odniesienia należy dokonać konwersji układu odniesienia pliku rastrowego zwracanego w formie WKT na format akceptowany przez wektor, a następnie o utworzenie wektora oraz przynajmniej jedno pole, które będzie przechowywało wartość z wektoryzacji.

In [15]:
srs = osr.SpatialReference()
srs.ImportFromWkt(raster.GetProjection())
Out[15]:
0
In [16]:
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource("polygon.shp")
dst_layer = dst_ds.CreateLayer("layer", srs = srs )
newField = ogr.FieldDefn('class', ogr.OFTInteger)
dst_layer.CreateField(newField)
Out[16]:
0

Ostatnim krokiem jest poligonizacja, która w bibliotece GDAL jest specjalizowaną funkcją, wymaga podania dygitalizowanej warstwy rastrowej, docelowej warstwy wektorowej oraz numer kolumny do której zostaną zapisane wyniki dygitalizacji.

In [17]:
gdal.Polygonize(tband, None, dst_layer, 0, [], callback=None )
dst_ds=None
Out[17]:
0

Zadanie:

W zależności od wielkości rastra zmodyfikować skrypt, tak aby dla dużych rastrów, raster był zapisywany na dysku a następnie po zakończonej wektoryzacji usuwany. Na podstawie poprzedniego przykładu wprowadzić do skryptu więcej przedziałów progowych niż 1.

Przykład 3: Przetwarzanie rastra fragmentami

Procedury geo-przetwarzania bardzo często rozmiar danych rastrowych jest zbyt duży aby przechowywać całość danych w pamięci. Drugim mechanizmem wymuszającym porcjowanie danych jest proces zrównoleglania obliczeń. Z tego powodu stosuje się mechanizmy przetwarzania porcjami (chunks).

In [18]:
os.chdir('/home/jarekj/Downloads/landsat/')
nir_path = "NIR.tif"
red_path = "RED.tif"

Mechanizm porcjowania wymaga wiedzy na temat rozmiaru zbioru rastrowego oraz określenia wielkości porcji przetwarzanych danych. Raster ma rozmiary 5000 x 5000 Tu ustalamy wielkość porcji na 1000x1000 co oznacza 25 niezależnie przetwarzanych porcji. Celem przetwarzania jest wyliczenie indeksu NDVI, wg wzoru: (NIR-RED)/(NIR+RED)

In [19]:
ref = gdal.Open(nir_path)
x = ref.RasterXSize
y = ref.RasterYSize
proj = ref.GetProjection()
gt = ref.GetGeoTransform()
step=1000

W pierwszej kolejności zbudowane zostaną dwie listy: lista przesunięć (od 0 do wielkości wymiaru z przesunięciem o wartość kroku) oraz dwie listy rozmiarów porcji. W zależności czy wielkość rastra dzieli się bez reszty przez wielkość porcji czy nie dodajemy na końcu mniejszą porcję o wielkości tej reszty (modulo z dzielenia)

In [20]:
xoffs = list(range(0,x,step))
yoffs = list(range(0,y,step))
# wyliczenie zakresów
ysizes = [step]*(y//step)
if y%step: ysizes.append(y%step)

xsizes = [step]*(x//step)
if x%step: xsizes.append(x%step)

Otwaramy dostęp do plików i pobieramy z nich pasmo z danymi

In [21]:
nir_rast = gdal.Open(nir_path)
nir_band = nir_rast.GetRasterBand(1)
red_rast = gdal.Open(red_path)
red_band = red_rast.GetRasterBand(1)

Oraz tworzymy nowy plik w którym będziemy zapisywać dane. Plik jest typu Float32

In [22]:
ndvi_rast = gdal.GetDriverByName('GTiff').Create("ndvi.tif", x,y, 1, gdal.GDT_Float32)
ndvi_rast.SetProjection(proj)
ndvi_rast.SetGeoTransform(gt)
nband=ndvi_rast.GetRasterBand(1)
Out[22]:
0
Out[22]:
0
  1. Pętla główna przechodzi po wszystkich porcjach wzdłuż wierszy (A) a następnie dla każdego wiersza porcji wzdłuż kolumn (B).
  2. Tworzona jest pusta tablica danych o wymiarach 2*liczba_komorek_w_porcji typu int16 czyli takiego jakiego są dane źródłowe.
  3. Pobierane są dane dla każdej porcji w postaci spłaszczonej
  4. Tworzone są dwa wiersze selekcji: te które zawierają dane not-null (A) oraz te, dla których obie wartości nie są zerami (B), aby uniknąć dzielenia przez zero w podanej formule.
  5. Po połączeniu danych z tablicy danych pobierane są tylko te, które spełniają oba warunki
  6. Tworzymy jednowymiarową tablicę danych wypełnioną wartościami -9999, typu zmiennoprzecinkowego
  7. Wyliczamy wartość ndvi dla każdej komórki (konwersja do float32 jest niezbędna w pythonie 2) i zapisujemy wyniki tylko do tych wartości które spełniają warunek kompletności
  8. Modyfikujemy wymiary tablicy do tej zgodnej z wielkością porcji i wyniki zapisujemy do rastra docelowego
In [23]:
for xoff,xsize in zip(xoffs,xsizes): #1A
    for yoff,ysize in zip(yoffs,ysizes): #1B
        data=np.empty([2,xsize*ysize],dtype='int16') #2
        data[0] = nir_band.ReadAsArray(xoff,yoff,xsize,ysize).flatten() #3
        data[1] = red_band.ReadAsArray(xoff,yoff,xsize,ysize).flatten()
        complete1 = (data != -9999).any(0) #4A
        complete2 = (data != 0).all(0) #4B
        complete = complete1 & complete2
        complete_data = data[:,complete] #5
        ndvi = np.full(len(complete),-9999,dtype='float32') #6
        values = (complete_data[0]-complete_data[1])/(complete_data[0]+complete_data[1]).astype('float32') #7
        ndvi[complete] = values 
        ndvi.shape = (ysize,xsize) #8
        nband.WriteArray(ndvi,xoff,yoff)
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0
Out[23]:
0

Ostatnim krokiem jest zdefiniowanie wartości pustej i zamknięcie pliku.

In [24]:
nband.FlushCache()
nband.SetNoDataValue(-9999)
del nband
ndvi_rast = None
ref = None
Out[24]:
0

Zadanie:

Dowiedzeić się więcej o funkcji ReadAsArray() a następnie zbudować nowy skrypt, który w jednym kroku wyliczy podgląd ndvi w (dużo) mniejszej rozdzielczości i wyświetli w postaci obrazka z dobraną legendą. Sam skrypt tak zmodyfikować, aby użytkownik mógł wskazać mniejszy zakres rastra a następnie wielkość kroku została dobrana do tego rastra

przykład 4: Przetwarzanie pliku wektorowego

Celem ćwiczenia jest zmodyfikowanie atrybutów pliku wektorowego, w taki sposób, że do każdego wiersza zostaną dodane dwa nowe atrybuty: najbliższy obiekt oraz odległość do niego. W tym celu otworzymy warstwę "liczniki.shp" oraz sprawdzimy jakie pola w danej warstwie są dostępne.

In [25]:
os.chdir('/home/jarekj/Dropbox/dydaktyka/programowanie2/data_vector')
vect = ogr.Open("gauges.shp",0) # 1 dla update
layer = vect.GetLayer()
layer.GetGeomType() == ogr.wkbPoint # test logiczny

defn = layer.GetLayerDefn()
nfields = defn.GetFieldCount()
nfields
Out[25]:
True
Out[25]:
2
In [26]:
fnames = [defn.GetFieldDefn(i).GetName() for i in range(nfields)]
types = [defn.GetFieldDefn(i).GetTypeName() for i in range(nfields)]
fields = dict(zip(fnames, types))
fields # w formie dict
Out[26]:
{'id': 'Integer64', 'nazwa': 'String'}

Następnie zarówno pola jak i geometria zostaną zapisane w odpowiednich tablicach.

In [27]:
layer.ResetReading()
attrs=[]
points=[]
for f in layer: # f - feature
    attrs.append([f.GetField(i) for i in range(nfields)])
    points.append(f.GetGeometryRef().GetPoint())


arr_points = np.array(points) # ograniczenie do pierwszych 10 punktów
arr_points = arr_points[:,0:2] # usunięcie pustego Z
arr_points
Out[27]:
array([[354592.23826591, 503273.00136985],
       [355907.25249054, 507371.46236992],
       [352532.049314  , 511272.67123631],
       [358186.61047988, 511469.92337   ],
       [360356.38395051, 513464.36161068],
       [364213.75900941, 505749.61149289],
       [362219.32076873, 499700.54605962],
       [359545.458512  , 499327.95869598],
       [360597.46989169, 503974.34228965],
       [360882.38964036, 507897.46805977],
       [356696.26102531, 504631.84940196],
       [351567.70554928, 507524.88069613],
       [354877.15801458, 510593.24722025],
       [366756.11984368, 501519.64907035],
       [347863.74881659, 513398.61089945],
       [351984.12672041, 504281.17894206]])
  1. Do obliczenia odległości użyjemy funkcji z modułu scipy.spatial.distance w celu utworzenia macierzy odległości między punktami.
  2. Macierz trzeba zamienić na formę symetryczną (kwadratową), gdyż scipy domyślnie tworzy macierz w postaci jedynie połówkowej. Po konwersji usuwamy przekątną gdzie znajdują się wartości 0 (a więc najmniejsze)
  3. Dla każdego wiersza wyliczamy pozycję wartości minimalnej odległości (używamy funkcji obsługującej wartości nan)
  4. Wydobywamy niezbędne wartości z tabeli atrybutów (przekształcając je funkcją zip() tak aby odwrócić wiersze z kolumnami
  5. Wybieramy te nazwy, których indeksy zostały wskazane
  6. Wyliczamy wartości minimalnych odległości
In [28]:
from scipy.spatial import distance

dst = distance.pdist(arr_points) #1
dst = distance.squareform(dst) #2
np.fill_diagonal(dst,np.nan)  

indexes = np.nanargmin(dst,axis=1) #3 opcja z nan
namelist = list(zip(*attrs))[1] #4
names = [namelist[i] for i in list(indexes)] #5
diststo = np.nanmin(dst,axis=1) #6

Wyniki możemy wyświetlić w postaci podwójnej listy

In [29]:
list(zip(names,diststo))
Out[29]:
[('Łazarz', 2504.671584546102),
 ('Łazarz', 2850.9671485636572),
 ('Strzeszyn', 2441.5470118595254),
 ('Umultowo', 2947.1513041790754),
 ('Piątkowo', 2947.1513041790754),
 ('Stare Miasto', 3963.7494503216917),
 ('Starołęka', 2699.6964109939313),
 ('Krzesiny', 2699.6964109939313),
 ('Stare Miasto', 3933.4584110435662),
 ('Rataje', 3933.4584110435662),
 ('Górczyn', 2504.671584546102),
 ('Junikowo', 3270.3222564645775),
 ('Krzyżowniki', 2441.5470118595254),
 ('Krzesiny', 4887.911784199176),
 ('Krzyżowniki', 5129.585654364245),
 ('Górczyn', 2796.1880929034824)]

Następnie budujemy nowy wektor według znanego schematu, kopując do niego istniejące pola oraz dodając dwa nowe: "colsest" i "dist_to"

In [30]:
os.environ['SHAPE_ENCODING'] = "utf-8"
# create layer based on the exisiting layer

driver =  driver = ogr.GetDriverByName("ESRI Shapefile")
new_vect = driver.CreateDataSource("gauges2.shp")
new_layer = new_vect.CreateLayer("layer", layer.GetSpatialRef(),layer.GetGeomType())

for i in range(nfields): 
    new_layer.CreateField(defn.GetFieldDefn(i)) # kopiowanie pól z istniejącego shp

new_layer.CreateField(ogr.FieldDefn("closest", ogr.OFTString)) #Nowe pola
new_layer.CreateField(ogr.FieldDefn("dist_to", ogr.OFTReal))
Out[30]:
0
Out[30]:
0
Out[30]:
0
Out[30]:
0

Następnie wypełaniamy wektor przez features: dodajemy geometrię kopiując ją z istniejącego rastra oraz kopiujemy zawartość istniejących pól, na końcu dodając nowe. Wartości zmiennych "closest" i "dist_to" pobieramy z tablic names i distto poprzez index, którego wartość zwiększamy w każdym kroku (index+=1)

Po zakińczeniu pętli głównej zamykamy wektor.

In [31]:
index=0
new_layer.ResetReading()
for f in layer:
#f = layer.next()
    new_feature = ogr.Feature(new_layer.GetLayerDefn())
    
    geom = f.GetGeometryRef() # kopiowanie geometrii
    new_feature.SetGeometry(geom)
#i=0    
    for i in range(nfields): # kopiowanie pól
        new_feature.SetField(fnames[i],f.GetField(i))
    
    new_feature.SetField("closest",names[index])
    new_feature.SetField("dist_to",diststo[index])
    new_layer.CreateFeature(new_feature)
    new_feature = None
    index+=1

new_vect=None
new_layer=None
vect=None

Zadanie

Zmodyfikować skrypt tak, aby zidentyfikować stopień izolacji punktów (średnia odległość do 4 najbliższych). Wartość izolacji oraz listę czterech najbliższych (jako jeden string) zapisać jako atrybuty wektora.

Przykład 5: Profil na podstawie linii wektorowej

Celem ćwiczenia jest wykonanie profil ze zbioru danych rastrowych wzdłuż linii zapisanej w formie pliku wektorowego.

W tym celu otwieramy dostęp do zbioru rastrowego i wektorowego

In [32]:
os.chdir('/home/jarekj/Dropbox/dydaktyka/programowanie2')
vect_path = "data_vector/przekroj.shp"
rast_path = "data_raster/pop_low.tif"

Wczytujemy dane wektorowe i dla pierwszej linii (feature) pobieramy jej geometrię. Jest to 6 punktów - wierzchołków linii

In [33]:
vect = ogr.Open(vect_path,0) 
layer = vect.GetLayer()
layer.GetGeomType() == ogr.wkbLineString # test logiczny

layer.ResetReading()
f = layer.GetFeature(0)
geom = f.GetGeometryRef()
geom.GetPoints()
Out[33]:
True
Out[33]:
[(971253.675438176, 1845409.8010865962),
 (977285.452210527, 1851578.6636946825),
 (985373.5165189068, 1855828.3246024752),
 (987087.0894655974, 1860831.9576068118),
 (995792.0400347859, 1865835.5906111484),
 (995792.0400347859, 1865835.5906111484)]

Taka liczba punktów jest zbyt mała aby wykonać przekrój. Zwiększamy liczbę punktów przy pomocy metody .Segmentize(), która dodaje punkty do linii co określony odstęp. Tu użyjemy 200 metrów. Po dodaniu punktów wyniki zapisujemy w tablicy points

In [34]:
geom.Segmentize(200)
points = np.array(geom.GetPoints())

Otwieramy raster i zapisujemy go jako tablicę dwuwymiarową.

In [35]:
raster = gdal.Open(rast_path)
band = raster.GetRasterBand(1)
gtr = raster.GetGeoTransform() # paramtery transformacji rastra
arr = band.ReadAsArray()

Następnie w celu uniknięcia pracy w dwóch wymiarach, współrzędne punktów przeliczamy na numer wiersza i komórki, oraz zamieniamy oczywiście ich typ na wartość całkowitą. W tym samym celu tablicę dwuwymiarową zawierającą dane rastra spłaszczamy a następnie wyliczamy indeksy komórek wg formuły: nr_wiersza * ilosc_kolumn + nr_kolumny. Ta formuła pozwala na obsługę rastra w postaci tablicy jednowymiarowej.

W ostatnim kroku pobieramy wartości rastra poprzez indeks komórki.

In [36]:
px = (points[:,0] - gtr[0]) / gtr[1]#x pixel
py = (points[:,1] - gtr[3]) / gtr[5] #y pixel
px = px.astype(int)
py = py.astype(int)
index = py * raster.RasterXSize + px
values = arr.flatten()[index]

Aby prawidłowo wykonać wykres należy obliczyć odległości pomiędzy kolejnymi punktami. W tym celu wykoamy dwie tabice:

  • tablicę pozbawioną ostatniego elementu
  • tablicę pozbawioną pierwszego elementu, a więc przesuniętą względem poprzedniej tablicy o 1

Dzięki temu zabiegowi możemy obliczyć odległości pomiędzy punktami w formie prostej zwektoryzowanej formuły odległości euklidesowej: $ \sqrt{(x_1-x_2)^2+(y_1-y_2)^2} $. Po wyliczeniu odległości musimy dodać na początku 0 (dla pierwszego punktu) oraz odległości przeliczyć na odległości skumulowane

In [37]:
points1 = points[:-1,:]
points2 = points[1:,:]

dists = np.sqrt((points1[:,0]-points2[:,0])**2 + (points1[:,1]-points2[:,1])**2)
dists = np.insert(dists,0,0)
cum_dists = np.cumsum(dists)

Ostatnim krokiem jest wykreślenie wykresu.

In [38]:
plt.plot(cum_dists,values)
Out[38]:
[<matplotlib.lines.Line2D at 0x7fe135fa8fd0>]

Zadanie:

Zmodyfikować skrypt w taki sposób, aby na wykresie dodatkowo oznaczyć miejsca załamania linii przekrojowej