Wstęp (Projekt w fazie rozwoju)

R jest językiem zaprojektowanym do przetwarzania danych znajdujących się w pamięci operacyjnej. Niestety dane geo-przestrzenne często są zbyt duże, aby mogły być przetwarzane bezpośrednio w pamięci. Systemy informacji geograficznej wykorzystują algorytmy, które pozwalają ładować do pamięci jedynie niezbędną część danych wykonać obliczenia a następnie zapisać wynik na dysk. Pakiet R pozwala na podobne podejście. W tym zadaniu zapoznamy się z metodami przetwarzania danych rastrowych w taki sposób, aby do pamięci w czytywać jedynie niezbędną część danych. W drugiej części opracujemy metodę wczytywania optymalnej ilości danych.

Własności zbioru geoprzestrzennego i bezpośredni dostęp do danych

R oprócz licznych funkcji będących procedurami GIS dostarcza też wiele funkcji, które możemy uznać za funkcje niższego poziomu. Pierwszą grupę stanowią funkcje pobierające informację o topologii rastra i służą one do pracy na geometrycznej części danych, druga grupa funkcji pozwala na pobieranie danych. Część tych funkcji przewinęła się w innych częściach kursu, tu zostaną omówione pod względem funkcjonalnym.

Funkcje opisujące i modyfikujące topologię.

Rozmiar rastra definiowany jest funkcjami nrow() ncol() oraz ncell(), która jest iloczynem dwóch poprzednich. Funkcje te pozwalają na zmodyfikowanie parametrów rastra ale nie modyfikują zbioru. Ich użycie jako modyfikatorów ma sens tylko w sytuacji, gdy jest to część złożonej procedury. Funkcje modyfikujące raster wywołują je w sposób niejawny. Podobną rolę pełnią funkcje res() oraz xres() i yres() zwracające i modyfikujące rozdzielczość.

library(raster)
## Loading required package: sp
rast <- raster("DATA/pop_low.tif")
ncol(rast) #zwraca wartości
## [1] 201
ncol(rast) <- 100
plot(rast) # zgłosi błąd, parametry są niespójne
## Error in .plotraster2(x, col = col, maxpixels = maxpixels, add = add, : no values associated with this RasterLayer
rast <- raster("DATA/pop_low.tif") # reset
res(rast)
## [1] 240 240
res(rast) <- c(120,120)
plot(rast) # zgłosi błąd, parametry są niespójne
## Error in .plotraster2(x, col = col, maxpixels = maxpixels, add = add, : no values associated with this RasterLayer
rast <- raster("DATA/pop_low.tif") # reset

Modyfikacji rozdzielczości możemy dokonać przy pomocy funkcji aggregate(), lub przy pomocy funkcji resample(), definiując wcześniej nową topologię rastra:

rast <- raster("DATA/pop_low.tif") # reset
new_rast <- raster(nrows=nrow(rast)/2.2, ncols=ncol(rast)/2.2, ext=extent(rast)) # uwaga na niekonsekwencję nazw
new_rast <- resample(rast,new_rast,method="bilinear") #przepróbkowanie
res(new_rast) #dość nietypowe i niesymetryczne, ale jak najbardziej prawidłowe
## [1] 530.1099 525.5172
plot(new_rast) # raster o innej topologii.

Zasięg geograficzny opisywany jest funkcją extent(), funkcja ta może przyjmować dwie postacie: zasięg geograficzny definiowany przez cztery wartości, lub w przypadku powiązania z obiektem raster pozwala na zdefiniowanie zasięgu geograficznego na podstawie wierszy i kolumn. Ta druga funkcja jest potrzebna, gdy przetwarzamy raster jak macierz (tj bez kontekstu geograficznego) i potrzebujemy ten kontekst dla wybranego fragmentu zbioru danych. Na przykład gdy dzielimy raster na kilka obszarów. Extent jest przedstawiany w jednostkach geograficznych (najczęściej metrach) lub w stopniach dla układów LonLat.

rast <- raster("DATA/pop_low.tif")
rast
## class       : RasterLayer 
## dimensions  : 127, 201, 25527  (nrow, ncol, ncell)
## resolution  : 240, 240  (x, y)
## extent      : 953085, 1001325, 1838625, 1869105  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : /home/jarekj/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA/pop_low.tif 
## names       : pop_low 
## values      : 0, 1507.513  (min, max)
ext1 <- extent(rast,20,50,10,100)
ext1 # zasięg mniejszego rastra
## class       : Extent 
## xmin        : 955245 
## xmax        : 977085 
## ymin        : 1857105 
## ymax        : 1864545
rast1 <- crop(rast,ext1)
rast1 #wypisuje wszystkie parametry, które zostały zmienione
## class       : RasterLayer 
## dimensions  : 31, 91, 2821  (nrow, ncol, ncell)
## resolution  : 240, 240  (x, y)
## extent      : 955245, 977085, 1857105, 1864545  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       : pop_low 
## values      : 0, 266.0363  (min, max)
plot(rast1) #wycięty fragment

projection(rast)
## [1] "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
crs(rast)
## CRS arguments:
##  +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0
## +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
isLonLat(rast) # test czy LatLon, a więc inne jednostki
## [1] FALSE

Pobieranie wartości i obliczenia

Przykład: zastosowanie operacji na topologii do porównania fragmentów rastra ze wzorcem.

Celem zadania jest analiza porównawcza zbioru rastrowego podzielonego na fragmenty o stałej wielkości. To zilustrowania zagadnienia użyjemy testu statystycznego t-Studenta a więc typowej operacji statystycznej. Ograniczeniem tego postępowania jest wymóg dostosowania ilości wierszy i kolumn rastra do wielkości okna próbkowania. W pierwszym kroku wyznaczamy parametry: wielkość okna (w pikselach), współrzędne środka okna próbkowania oraz przetwarzany raster.

library(raster)
size <- 30 #rozmiar kafelka
r <- 70 #wspórzędne środka
c <- 130
rast <- raster("DATA/pop_low.tif")

W następnym kroku modyfikujemy raster tak aby był zgodny z siatką $ 30 30$ oraz budujemy data.frame zawierający indeksy siatki.

nrows <- round(nrow(rast)/size)*size #zaokrąglenie do size
ncols <- round(ncol(rast)/size)*size
new_rast <- raster(nrows=nrows, ncols=ncols, ext=extent(rast))  #topologia po zaokrągieniu
new_rast <- resample(rast,new_rast,method="bilinear") #resampling
new_rast #nowy raster, jego wymiary są wielokrotnością 30
## class       : RasterLayer 
## dimensions  : 120, 210, 25200  (nrow, ncol, ncell)
## resolution  : 229.7143, 254  (x, y)
## extent      : 953085, 1001325, 1838625, 1869105  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : pop_low 
## values      : -0.3076951, 966.7551  (min, max)
#budujemy data.frame zawierający indexy siatki
rows <- seq(1,nrow(new_rast),by=size) 
cols <- seq(1,ncol(new_rast),by=size)
blocks <- expand.grid(rows=rows,cols=cols) #data frame z indeksami

Wyłącznie w celu wizualizacji pokazujemy lokalizację próbki, wprowadzając w to miejsce absurdalne wartości

#tylko do wizualizacji lokalizacji zaznaczamy pobierany obraz
new_rast2 <- new_rast
new_rast2[(r-size/2+1):(r+size/2),(c-size/2+1):(c+size/2)] <- -1000
image(new_rast2)

#koniec wizuajizacji

Proces analizy zostanie przedstawiony w sposób bardzo zwięzły, wykorzystujący własności wektoryzacji. W pierwszym kroku pobieramy wartości z próbki: $ 30 30 = 900 $ wartości a następnie na podstawie data.frame wykonujemy macierz o wymiarach $ ilość .komórek .siatki (35) ilość .pikseli .w .komórce (900) $ stosując funkcję apply po wierszach data.frame blocks

sample <- getValuesBlock(new_rast,r-size/2,size,c-size/2,size) #próbka, którą testujemy
values <- apply(blocks,1, # poszczególne komórki
                function(blocks,rast,size) {getValuesBlock(rast,blocks[1],size,blocks[2],size)}, new_rast,30)

W ostatnim kroku stosujemy funkcję anonimową po kolumnach macierz values Istotą tej funkcji, oprócz wyliczenia statystyki t-Studenta, jest sprawdzenie (iflese()), aby nie wyliczać statystyk, jeżeli ilość wartości pustych w komórce jest powyżej pewnej wartości (tu więcej niż 90%). W ostatnim kroku wizualizujemy macierz.

statistics <- apply(values,2,function(values,sample) 
  {ifelse(sum(is.na(values))>0.9*length(values),NA, t.test(values,sample)$statistic)},  sample) 
#musimy zabezpieczć się przed sytuacją, że brakuje danych do testów

# Ostatni krok wizualizacja macierzy:
mat <- matrix(statistics,length(rows),length(cols)) #wynik zamieniamy na macierz i drukujemy....
mat <- t(mat[length(rows):1,]) # zgodnie z kierunkami geograficznymi, obraz musi być odbity względem y a następnie przeprowadzona transpozycja

library(ggplot2)
library(reshape2)
mmat <- melt(mat)
ggplot(mmat,aes(x=Var1,y=Var2,fill=value)) + geom_tile()

image(mat,col=terrain.colors(20)) # lub funkcją image(), która jest prostrza niż ggplot

Następnym krokiem może być przekształcenie operacji w funkcję.

Przetwarzanie danych w pamięci i na dysku

Czytanie danych.

Problem zostanie zilustrowany poprzez napisanie prostej funkcji, która przeliczy liczbę ludności w komórce na liczbę ludności na metr kwadratowy. Do wykonania ćwiczenia użyjemy uproszczonej wersji mapy ludnościowej hrabstwa Hamilton

library(raster)
pop <- raster("DATA/pop_low.tif")

Prosta funkcja działająca w pamięci

. W pierwszej kolejności przyjrzymy się charakterystyce obiektu. Interesuje nas rozdzielczość: tu jest to 240 na 240 metrów.

pop
## class       : RasterLayer 
## dimensions  : 127, 201, 25527  (nrow, ncol, ncell)
## resolution  : 240, 240  (x, y)
## extent      : 953085, 1001325, 1838625, 1869105  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : /home/jarekj/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA/pop_low.tif 
## names       : pop_low 
## values      : 0, 1507.513  (min, max)

Do pozyskania rozdzielczości służy funkcja res, a powierzchnię komórki możemy obliczyć, mnożąc przez siebie obie liczby. W przypadku projekcji płaskich nie należy stosować funkcji area, która służy do obliczenia warstwy rastrowej zawierającej powierzchnie wszystkich komórek, co jest użyteczne w przypadku danych w projekcji lat-lon.

res(pop)
## [1] 240 240
area <- res(pop)[1]*res(pop)[2]

Tak więc aby wyliczyć gęstość zaludnienia, należy podzielić liczbę ludności przez powierzchnię komórki:

pop.to.dens.ram <- function(pop) {
  area <- res(pop)[1]*res(pop)[2]
  stopifnot(area>0) #zatrzymanie funkcji jeżeli area była by 0, w zasadzie niemożliwe
  dens <- pop/area
  dens
}
#wykonanie funkcji i wydruk wyniku
dens <- pop.to.dens.ram(pop)
par(mfrow=c(1,2))
plot(pop)
plot(dens)

par(mfrow=c(1,1))

Funkcje przechowująca dane na dysku

Powyższa funkcja realizuje wszystkie czynności w poleceniu dens <- pop/area tworzy dane w pamięci operacyjnej, a więc w przypadku dużego zbioru i niewielkich zasobów może okazać się że obliczenia się nie powiodły. Z drugiej strony operacje dyskowe są zawsze wolniejsze niż operacje w pamięci i nie ma sensu spowalniać obliczeń jeżeli mogą być one z powodzeniem przeprowadzone w pamięci. Kolejne kroki pokażą jak przygotować funkcję, która będzie mogła pracować zarówno z danymi w pamięci jak i na dysku. W tym celu należy zapoznać się z funkcjami update(), writeStart(), writeStop(), writeValues(), getValues() oraz blockSize() w pliku pomocy.

pop.to.dens.disk <- function(pop,filename) { #filename musi być podane jako ścieżka dostępu
  area <- res(pop)[1]*res(pop)[2]
  dens <- raster(pop)
  bs <- blockSize(dens) #wielkość bloku danych przechowywanego w pamięci
  dens <- writeStart(dens,filename,overwrite=T) #musimy włączyć prawo do nadpisywania wyniku
  for(i in 1:bs$n) { #dla każdego bloku
    tmp <- getValues(pop,row=bs$row[i],nrows=bs$nrows[i])
    tmp <- tmp/area
    dens <- writeValues(dens,tmp,bs$row[i])
  }
  dens <- writeStop(dens)
  return(dens)
}

#wykonanie funkcji i wydruk wyniku
dens <- pop.to.dens.disk(pop,"DATA/tmp.tif")
dens #zwrócić uwagę na data source
## class       : RasterLayer 
## dimensions  : 127, 201, 25527  (nrow, ncol, ncell)
## resolution  : 240, 240  (x, y)
## extent      : 953085, 1001325, 1838625, 1869105  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : /home/jarekj/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA/tmp.tif 
## names       : tmp 
## values      : 0, 0.0261721  (min, max)

Funkcja wybierająca sposób zapisu w zależności od dostępnej pamięci

Poprzednia funkcja będzie przechowywać wyniki na dysku niezależnie od tego czy jest to niezbędne czy nie. Kolejne rozszerzenie funkcji pozwoli wybrać sposób przechowywania danych w zależności od ilości dostępnej pamięci. Należy zapoznać się z funkcjami: canProcessInMemory oraz rasterTempFile. użytkownik może wymusić zapis na dysk niezależnie od wielkości danych poprzez podanie nazwy pliku (domyślnie jest pusta)

pop.to.dens.auto <- function(pop,filename='') { #filename musi być podane jako ścieżka dostępu
  area <- res(pop)[1]*res(pop)[2]
  dens <- raster(pop)
  mem <- canProcessInMemory(dens,n=2) #n - number of data copies in memory
  filename <- trim(filename) #usuwa białe znaki
  
  if(filename!='') mem <- FALSE #jeżeli jest filename; wymuszamy zapis na dysku, niezależnie czy konieczny czy nie
  if(mem==FALSE & filename=='') filename <- rasterTmpFile() #tworzy plik tymczasowy w katalogu roboczym
  
  if(mem==TRUE) {
    dens <- pop/area
  } else {
    bs <- blockSize(dens) 
    dens <- writeStart(dens,filename,overwrite=T)
    for(i in 1:bs$n) { 
      tmp <- getValues(pop,row=bs$row[i],nrows=bs$nrows[i])
      tmp <- tmp/area
      dens <- writeValues(dens,tmp,bs$row[i])
    }
    dens <- writeStop(dens)
  }
  return(dens)
}

#wykonanie funkcji i wydruk wyniku
dens <- pop.to.dens.auto(pop,"DATA/tmp.tif")
dens #zwrócić uwagę na data source
## class       : RasterLayer 
## dimensions  : 127, 201, 25527  (nrow, ncol, ncell)
## resolution  : 240, 240  (x, y)
## extent      : 953085, 1001325, 1838625, 1869105  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : /home/jarekj/Dropbox/dydaktyka/prog_geo/PUBLISH/DATA/tmp.tif 
## names       : tmp 
## values      : 0, 0.0261721  (min, max)
dens <- pop.to.dens.auto(pop)
dens #zwrócić uwagę na data source
## class       : RasterLayer 
## dimensions  : 127, 201, 25527  (nrow, ncol, ncell)
## resolution  : 240, 240  (x, y)
## extent      : 953085, 1001325, 1838625, 1869105  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## data source : in memory
## names       : pop_low 
## values      : 0, 0.0261721  (min, max)