Integracja z QGIS
Wstęp
Skrypty tworzy się w celu wielokrotnego wykonywania tego samego kodu. Ale nawet najlepiej napisane skrypty są trudne w codziennym użyciu dla osób, które nie są specjalistami danego języka programowania. Dlatego umiejętność napisania skryptu w sposób umożliwiający jego bezproblemową obsługę wszystkim użytkownikom systemów integracji geograficznej stanowi krytyczny krok do odniesienia sukcesu w przestrzeni rynkowej. W tej części kursu zapoznamy się jak przygotować skrypt, aby zintegrować go z popularnym programem QuantumGIS ( QGIS).
Przygotowanie środowiska QuantumGIS i przygotowanie skryptu
- [W systemie Windows: Dodać pusty katalog msys w katalogu: C:Files (x86)2.18]
- Wybrać Geoprocessing -> opcje -> R scripts:
- Aktywować; [W systemie Windows: Ustawić ścieżkę do R (:Files-3-3-2) lub inną. Nie do bin!]
- W Narzędziach geoprocessingu powinien się pojawić R scripts
Skrypty są tworzone i zapisywane w specjalnym oknie dialogowym. Dodatkowo po zapisaniu należy przygotować opis skryptu, w tym opis parametrów:
Po zapisaniu pod własną nazwą, skrypt powinien być gotowy do uruchomienia
Pisanie skryptów
Skrypt w który ma być uruchamiany w środowisku GGIS różni się nieco od czystego skryptu R. Przede wszystkim, skrypt QGIS nie ma części odpowiedzialnej za czytanie i zapis danych, tym zajmują się funkcje brick() i readOGR() oraz writeRaster() i writeOGR() z pakietów raster i rgdal. W związku z tym nie ma potrzeby definiować ich w skrypcie, jak również wczytywać tych dwóch bibliotek. Pozostałe użyte biblioteki muszą być wczytane jawnie poleceniem library(). Użytkownik nie musi mieć zainstalowanej biblioteki. W takiej sytuacji, skrypt doinstaluje bibliotekę przy pierwszym uruchomieniu.
Skrypt składa się z bloku parametrów i ciała skryptu przypominającej ciało funkcji. Ciało skryptu może zawierać własne funkcje. Większość rozwiązań dostępnych w języku R jest również dostępna w skryptach QGIS.
Parametry
Skrypt R posiada zestaw parametrów wejściowych oraz wyjście (może być wiele). Zarówno wejście jak i wyjście ma kilka typów możliwych parametrów. Parametry wypisywane są kolejno linia po linii, niektóre parametry pozwalają na zdefiniowanie wartości domyślnej.
Paramtery skryptu
Nazwa Grupy= group – zawiera nazwę grupy w której pojawi się skrypt. pomiędzy nazwą a “=” nie może być spacji
parametry wejściowe
RasterLayer= raster – uwaga na małą literę
VectorLayer= vector
Table= table
Field_vector1= Field VectorLayer – uwaga na Dużą literę.
Field_table1= Field Table –Field musi pochodzić z warstwy wektorowej lub tabeli już zadeklarowanej
Liczba= Number – tylko jeden typ liczbowy. Brak typu integer, date itp.
Logic= Boolean typ logiczny zwraca TRUE lub FALSE
Napis= String – dowolny łańcuch liter (o ograniczonej długości).
ListaWyboru= Selection opt1;opt2;opt3;… – lista opcji do wyboru
Zasieg= extent – cztery wartości, xmin, xmax, ymin, ymax
Referencja= crs – kod w postaci “EPSG:2180”
parametry wyjściowe
OutputRaster= output raster – zawsze musi być słowo output na początku.
OutputVector= output vector
OutputTable= output table – tworzy plik csv, a nie tabelę ładowaną do środowiska
Każdy parametr składa się kilku opcji opcji:
Nazwa parametru=TypParametru wartosc_domyslna
parametr rozpoczyna się od ## Pierwsza część parametru to nazwa zmiennej; zmienna pod taką nazwą będzie używana w skrypcie. Druga część to typ określa jakiego typu jest zmienna. W przypadku zmiennej Field, należy podać nazwę wektora/tabeli z której dane pole pochodzi. każdy parametr może też mieć opcję optional jeżeli nie jest wymagany.
Skrypty mają jeszcze dwa specjalne typy wyjścia – konsolę i urządzanie graficzne. Konsola nie definiowana w nagłówku ale w ciele skryptu, natomiast wyjście graficzne jest definiowane logicznym parametrem showplots bez żadnych dodatkowych atrybutów.
Przykładowy nagłówek skryptu wygląda tak:
##Moje Skrypty= group
##Input= raster
##Output= output rasterCiało skryptu
Ciało skryptu to ciąg poleceń, które zostaną wykonane w trakcie jego działania. Ciało skryptu nie potrzebuje żadnych funkcji wczytywania danych. Przykładowe ciało funkcji może wyglądać w następujący sposób:
area <- res(Input)[1]*res(Input)[2]
Output<-Input/areaTworzenie skryptów w środowisku R
Niestety środowisko QGIS nie posiada dobrej jakości narzędzi wspomagających proces tworzenia skryptu jak i jego testowania. Z tego powodu złożone skrypty należy zbudować i przetestować w środowisku R, i jeżeli już działają, należy je dopiero przenieść do środowiska QGIS. Jeżeli jednak uruchamiamy skrypty poza środowiskiem QGIS, należy zadbać o biblioteki i procedury wczytywania i zapisu danych. Skrypt musi również sprawdzić parametry wejściowe, czy umożliwią zadziałanie skryptu lub zwrócą sensowny wynik. Taki test należy przeprowadzić przed rozpoczęciem intensywnych obliczeń. Przed opuszczeniem skryptu jeżeli tworzymy nową warstwę, należy też nadać jej CRS, kopiując go z warstwy wejściowej lub podać bezpośrednio jako parametr skryptu.
W celu zilustrowania problemu przygotujemy prosty skrypt, który wylicza środek ciężkości zbioru punktów i zapisuje wynik w postaci pliku shp. Bardzo ważnym elementem skryptu geo-przestrzennego,
setwd("~/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA/poznan") #tymczasowy katalog roboczy
require(rgdal) #musimy wczytać bibliotekę## Loading required package: rgdal## Loading required package: sp## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.0, released 2016/04/25
## Path to GDAL shared files: /usr/share/gdal/2.1
## Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.2-3dsn <- "poznan_shape.sqlite"
Layer <- readOGR(dsn,layer="gauges") #nazwy powinny być ogólne## OGR data source with driver: SQLite
## Source: "poznan_shape.sqlite", layer: "gauges"
## with 16 features
## It has 2 fields#==================== Właściwa część skryptu =============================#
srodek <- apply(Layer@coords,2,mean) #wylicz średnią
Output <- SpatialPointsDataFrame(rbind(srodek), #rbind conwertuje do matrix
data=data.frame(id="center"),
proj4string = Layer@proj4string) #CRS
#==================== Koniec skryptu =====================================#
plot(Layer) + points(Output,pch=1,col=2) # sprawdzenie czy działa
## numeric(0)Końcowy skrypt może wyglądać w sposób następujący:
##Testy= group
##Layer= vector
##Output= output vector
srodek <- apply(Layer@coords,2,mean) #wylicz średnią
Output <- SpatialPointsDataFrame(rbind(srodek), #rbind conwertuje do matrix
data=data.frame(id="center"),
proj4string = Layer@proj4string) #CRS W kolejnej części zostaną przedstawione przykłady ilustrujące przygotowanie skryptów ilustrujące różne zasady przetwarzania danych geo-przestrzennych.
Przetwarzanie raster - konsola
Prosty skrypt, który pobiera wartości wskazanych kwantyli rastra i wypisuje przedziały na konsoli. Kwantyle są przekazywane jako łańcuch rozdzielony przecinkami, skrypt musi dokonać konwersji na liczby i sprawdzić czy nie ma tam czegoś co nie jest liczbą i czy liczby \(\in [0,1]\).
setwd("~/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA") #tymczasowy katalog roboczy
require(raster) #musimy wczytać bibliotekę## Loading required package: rasterLayer <- brick("pop_low.tif") #nazwy powinny być ogólne
quants <- "0,0.1,0.3,0.6,0.9,1" # dostarczamy zmienną tekstową
#==================== Właściwa część skryptu =============================#
split <- as.numeric(unlist(strsplit(quants,","))) #musimy ją podzielić,unlist zamnienia listę na wektor
stopifnot(sum(is.na(split))==0) #kontrola poprawności argumentów
stopifnot(max(split)<=1 & min(split)>=0)
x <- getValues(Layer)
round(quantile(x,split,na.rm=TRUE),3) #> output na konsolę## 0% 10% 30% 60% 90% 100%
## 0.000 0.013 1.479 30.648 125.623 1507.513Gotowy skrypt
##Layer= raster
##quants= String
split <- as.numeric(unlist(strsplit(quants,",")))
stopifnot(sum(is.na(split))==0)
stopifnot(max(split)<=1 & min(split)>=0)
x <- getValues(Layer)
result <- round(quantile(x,split,na.rm=TRUE),3)
>resultPrzetwarzanie raster - obraz
Prosty skrypt dzielący dane na kilka równych grup i drukujący boxplot. W celu przyspieszenia działania skryptu zredukujemy ilość danych poprzez proces losowania. Jako parametry skryptu zostaną przekazane udział wylosowanych elementów i liczba klas. Skrypt musi sprawdzić czy udział losowanych elementów \(\in (0,1)\).
setwd("~/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA") #tymczasowy katalog roboczy
require(raster) #musimy wczytać bibliotekę
Layer <- brick("pop_low.tif") #nazwy powinny być ogólne
Sampling <- 0.1 #albo samopisujące się bo będą wyświetlane w oknie
Nclass <- 4
#==================== Właściwa część skryptu =============================#
stopifnot(Sampling>0 & Sampling<1) #zabezpieczenie przed błędeami
library(ggplot2) #tym razem biblioteka musi być załadowana jawnie
x <- getValues(Layer)
sampling <- round(length(x)*Sampling) #wyliczamy ilość próbek
rnd <- sample(na.omit(x),sampling)
quants <- seq(0,1,by=1/Nclass) #kwantyle z ilości
qt <- quantile(rnd,quants,na.rm=TRUE) # kwantyle, w celu podziału na grupy równej wielkości
cc <- cut(rnd,qt,include.lowest=TRUE) #podział na 4 klasy
data <- data.frame(pop=rnd,class=cc) # musimy utworzyć data frame
p <- ggplot(data,aes(class,pop))
p + geom_boxplot()
#==================== Koniec skryptu =====================================#Gotowy skrypt
##Layer= raster
##Sampling= Number 0.1
##Nclass= Number 4
##showplots
stopifnot(Sampling>0 & Sampling<1)
library(ggplot2)
x <- getValues(Layer)
sampling <- round(length(x)*Sampling)
rnd <- sample(na.omit(x),sampling)
quants <- seq(0,1,by=1/Nclass)
qt <- quantile(rnd,quants,na.rm=TRUE)
cc <- cut(rnd,qt,include.lowest=TRUE)
data <- data.frame(pop=rnd,class=cc)
p <- ggplot(data,aes(class,pop))
p + geom_boxplot()Przetwarzanie raster - wektor
Prosty skrypt, który pobiera raster i następnie zamienia na poligony te obszary, których wartość jest większa niż próg. Skrypt musi sprawdzić, czy próg nie jest większy niż maximum (tzn. czy cokolwiek powstanie). W ramach skryptu trzeba wykonać re-klasyfikację do formatu binarnego [0,1], a następnie wykonać polecenie zamiany na poligony. Należy pamiętać o skopiowaniu georeferencji
setwd("~/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA") #tymczasowy katalog roboczy
require(raster)
Layer <- brick("pop_low.tif")
Threshold <- 100 #albo samopisujące się bo będą wyświetlane w oknie
#==================== Właściwa część skryptu =============================#
stopifnot(Threshold<maxValue(Layer))
library(rgeos) #wymagana do operacji dissolve## rgeos version: 0.3-21, (SVN revision 540)
## GEOS runtime version: 3.5.0-CAPI-1.9.0 r4084
## Linking to sp version: 1.2-3
## Polygon checking: TRUEsel <- function(x) {x==1} #funkcja selekcjonująca poligony
Layer[is.na(Layer)] <- 0 # usunięcie wartości pustych
reclass_matrix <- rbind(c(minValue(Layer),Threshold,0), #macierz reklasyfikacji
c(Threshold,maxValue(Layer),1))
Layer1 <- reclassify(Layer,reclass_matrix)
breaks <- c(minValue(Layer),Threshold,maxValue(Layer)) #szybsza alternatywa i prostrza
Layer2 <- cut(Layer,breaks,right=TRUE)-1 #chcemy mieć tylko 0 i 1
Output <- rasterToPolygons(Layer1,fun=sel,dissolve = TRUE)
proj4string(Output) <- crs(Layer) #dwie różne funkcje bo rgdal==proj4string a raster==crs
#==================== Koniec skryptu =====================================#
plot(Layer1)
plot(Output,col=2)
Właściwy skrypt
##Layer= raster
##Threshold= Number 100
##Output= output vector
stopifnot(Threshold<maxValue(Layer))
library(rgeos)
sel <- function(x) {x==1}
Layer[is.na(Layer)] <- 0
breaks <- c(minValue(Layer),Threshold,maxValue(Layer))
Layer1 <- cut(Layer,breaks,right=TRUE)-1
Output <- rasterToPolygons(Layer1,fun=sel,dissolve = TRUE)
proj4string(Output) <- crs(Layer)Przetwarzanie wektor - pola tabeli
Skrypt znajduje dla każdego punktu najbliższego sąsiada i dodaje pola tabeli: id oraz odległość. Skrypt testuje unikalność pola, które identyfikuje nam punkty, oraz czy warstwa to punktu. Skrypt można rozbudować o własne nazwy dodawanych pól.
setwd("~/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA/poznan") #tymczasowy katalog roboczy
require(rgdal) #musimy wczytać bibliotekę
dsn <- "poznan_shape.sqlite"
Layer <- readOGR(dsn,layer="gauges") #nazwy powinny być ogólne## OGR data source with driver: SQLite
## Source: "poznan_shape.sqlite", layer: "gauges"
## with 16 features
## It has 2 fieldsField <- "id" #nazwa pola id, musimy mieć pole będące indywidualnym identyfikatorem
#==================== Właściwa część skryptu =============================#
stopifnot(class(Layer) == "SpatialPointsDataFrame") #test czy punkty
stopifnot(sum(duplicated(Layer@data[,Field]))==0) # test czy pole jest unikalne
library(rgeos) #wymagana do funkcji gDistance
distances <- gDistance(Layer,byid=TRUE)
diag(distances) <- NA #usuwamy 0 na przekątnej macierzy, aby nie myliły funkcji min. Mamy znaleść min>0
target <- apply(distances,1,which.min) #funkcja which.min znajduje index najmiejszej wartości
distance <- apply(distances,1,min,na.rm=TRUE)
Layer$distance <- distance
Layer$target <- Layer@data[,Field][target] # ta stuczka ma na celu wpisanie nazw pól
Output <- Layer
#==================== Koniec skryptu =====================================#
str(Output)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 16 obs. of 4 variables:
## .. ..$ id : int [1:16] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..$ nazwa : Factor w/ 16 levels "Chyby","Górczyn",..: 2 3 6 9 15 16 5 12 10 11 ...
## .. ..$ distance: num [1:16] 2505 2851 2442 2947 2947 ...
## .. ..$ target : int [1:16] 11 11 13 5 4 10 8 7 10 9 ...
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:16, 1:2] 354592 355907 352532 358187 360356 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
## ..@ bbox : num [1:2, 1:2] 347864 499328 366756 513464
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"Właściwy skrypt
##Layer= vector
##Field=Field Layer
##Output=output vector
stopifnot(class(Layer) == "SpatialPointsDataFrame")
stopifnot(sum(duplicated(Layer@data[,Field]))==0)
library(rgeos)
distances <- gDistance(Layer,byid=TRUE)
diag(distances) <- NA
target <- apply(distances,1,which.min)
distance <- apply(distances,1,min,na.rm=TRUE)
Layer$distance <- distance
Layer$target <- Layer@data[,Field][target]
Output <- LayerPrzetwarzanie danych z docelową projekcją
Celem jest prosty skrypt, który, zakłada, że użytkownik integruje dane z różnego źródła i musi zintegrować je do jednego układu współrzędnych. Moduł pozwala wykonać reprojekcję na bazie crs, wektora i rastra. Każdy jest opcjonalny. Priorytety: 1 CRS; 2 Raster; 3 wektor;
setwd("~/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA") #tymczasowy katalog roboczy
require(raster)
Layer <- brick("pop_low.tif")
TargetRaster <- brick("pop.tif")
TargetVector <- NULL
TargetProj <- crs("+init=epsg:4326") # WGS84
TargetProj <- NULL
#==================== Właściwa część skryptu =============================#
#Sprawdzenie czy jest minimum 1 opcjonalny argument:
stopifnot(is.null(TargetProj) | is.null(TargetRaster) | is.null(TargetVector))
# Stosujemy funcję, aby uniknąć ewaluacji kolejnych elementów gdy znaleźliśmy pierwszy dobry
# Reprojekcję wykonujemy jeżeli choćby jeden jest
find_projection <- function(Layer,TargetProj,TargetRaster,TargetVector) {
if(!is.null(TargetProj)) { # prriorytet
if(!compareCRS(TargetProj,projection(Layer))) return(TargetProj)
return(NULL)
} else if (!is.null(TargetRaster)) { # 2 priorytet
if(!compareCRS(projection(TargetRaster),projection(Layer))) return(projection(TargetRaster))
return(NULL)
}else if (!is.null(TargetVector)) { # 3 prjorytet
if(!compareCRS(projection(TargetVector),projection(Layer))) return(projection(TargetVector))
return(NULL)
} else { # wszystkie równe nie ma potrzeby robić reprojekcji
return(NULL)
}
}
Target_projection <- find_projection(Layer,TargetProj,TargetRaster,TargetVector)
if(is.null(Target_projection)) { # nie można stosować ifelse, jeżeli prowadzi do kopiowania złożonych obiektów
Output <- Layer
} else {
Output <- projectRaster(Layer,crs=Target_projection)
}
#==================== Koniec skryptu =====================================# Właściwy skrypt
##Layer= raster
##TargetProj= optional crs
##TargetRaster= optional raster
##TargetVector= optional vector
##Output= output raster
stopifnot(is.null(TargetProj) | is.null(TargetRaster) | is.null(TargetVector))
find_projection <- function(Layer,TargetProj,TargetRaster,TargetVector) {
if(!is.null(TargetProj)) {
if(!compareCRS(TargetProj,projection(Layer))) return(TargetProj)
return(NULL)
} else if (!is.null(TargetRaster)) {
if(!compareCRS(projection(TargetRaster),projection(Layer))) return(projection(TargetRaster))
return(NULL)
}else if (!is.null(TargetVector)) {
if(!compareCRS(projection(TargetVector),projection(Layer))) return(projection(TargetVector))
return(NULL)
} else {
return(NULL)
}
}
Target_projection <- find_projection(Layer,TargetProj,TargetRaster,TargetVector)
if(is.null(Target_projection)) {
Output <- Layer
} else {
Output <- projectRaster(Layer,crs=Target_projection)
}