Jakub Nowosad

• Zakład Geoinformacji, WNGiG, UAM
• nowosad@amu.edu.pl
• http://nowosad.github.io

Dane

• http://nowosad.github.io/projekty/pazuR/dane.rar

"Systemy Informacji Geograficznej to zintegrowany sieciowo zestaw specjalistów, danych, metod badawczych, sprzętu komputerowego i oprogramowania, które to elementy działają w kontekście instytucjonalnym." (D.J. Maguire i in. 1991, 1998, 2011)

Do głównych zadań systemów informacji geograficznej należą:

• Pozyskiwanie danych przestrzennych.

• Gromadzenie i udostępnanie danych przestrzennych

• Przetwarzanie danych przestrzennych

• Analiza danych przestrzennych

• Wizualizacja danych przestrzennych

1. Punkt (ang. point) - pojedyncza lokalizacja punktowa, zawierająca zazwyczaj przynajmniej dwie współrzędne.
2. Linia (ang. line) - obiekt składający się z uporządkowanego zbioru punktów, które połączone są prostymi liniami.
3. Poligon (ang. polygon) - obiekt składający się z uporządkowanego zbioru punktów, które połączone są prostymi liniami, a pierwszy i ostatni punkt posiada takie same współrzędne. W związku z tym, obiekty poligonowe posiadają powierzchnie. Obiekty poligonowe często też zawierają luki (ang. holes), symbolizujące inny rodzaj powierzchni, który nie powinien być wliczany do badanego obiektu.
4. Siatka (ang. grid) - zbiór punktów lub prostokątnych komórek (zwanych "oczkami siatki") zorganizowanych w regularną siatkę. Każda komórka może zawierać jedną lub więcej wartości, które ją reprezentują.

Najczęsciej stosowane są dwa rodzaje układów współrzędnych:

• Geograficzne - wykorzystujące długość i szerokość geograficzną (wartości kątowe - w stopniach) służace do określania położenia punktu.

• Prostokątne płaskie - przedstawiane w miarach liniowych (np. w metrach).

Najkompletniejszą bazą informacji dotyczących różnych układów współrzędnych jest strona http://www.spatialreference.org/.

Jedną z możliwości odczytu danych przestrzennych do R jest wykorzystanie funkcji readOGR oraz readGDAL z pakietu rgdal. Opierają się one o otwartą bibliotekę GDAL/OGR i pozwalają na wczytanie danych przestrzennych do pamięci RAM.

Pakiet rgdal wyróżnia możliwość odczytania znacznej liczby wektorowych i rastrowych formatów danych przestrzennych. Dostępne formaty można sprawdzić za pomocą funkcji ogrDrivers (dane wektorowe) oraz gdalDrivers (dane rastrowe):

library('rgdal')
ogrDrivers() # sprawdzenie dostępnych formatów danych wektorowych
gdalDrivers() # sprawdzenie dostępnych formatów danych rastrowych

Funkcja readOGR, pozwalająca na wczytanie danych wektorowych, wymaga zadeklarowania dwóch parametrów:

• dsn - folderu zawierającego plik wektorowy, w przypadku wykorzystywania pliku zawartego w folderze roboczym konieczne jest zadeklarowanie tego za pomocą '.'
• layer - nazwę pliku wektorowego, bez podawania jego rozszerzenia

Przy korzystaniu z danych, które zawierają, np. polskie znaki, należy dodatkowo zdefiniować sposób ich kodowania za pomocą parametru endcoding.

wektor_punkt <- readOGR(dsn = 'dane/bdo250', layer = 'lotn_point', encoding = 'Windows-1250')
wektor_punkt

## OGR data source with driver: ESRI Shapefile 
## Source: "dane/bdo250", layer: "lotn_point"
## with 120 features and 16 fields
## Feature type: wkbPoint with 2 dimensions

##         coordinates AREA PERIMETER LOTN_ LOTN_ID LOT TYP
## 1  (464500, 435700)    0         0     1      17   2   6
## 2  (297700, 516900)    0         0     2      24   2   6
## 3  (529700, 373800)    0         0     3      25   2 -99
## 4  (263800, 546600)    0         0     4      42   2 -99
## 5  (367100, 509400)    0         0     5      68   2   2
## 6  (718000, 253100)    0         0     6      78   2 -99
## 7  (374900, 740900)    0         0     7      81   2   2
## 8  (209800, 623000)    0         0     8      86   2   2
## 9  (344400, 346800)    0         0     9      99   2   2
## 10 (515700, 211600)    0         0    10     107   2   2
## 11 (816200, 373900)    0         0    11  214361   2 -99
## 12 (313200, 511900)    0         0    12  105932   2   5
## 13 (731300, 195000)    0         0    13  306196   2   3
## 14 (584600, 610300)    0         0    14   89229   2   5
## 15 (518600, 436100)    0         0    15       1   2 -99
## 16 (443400, 495200)    0         0    16      32   2   2
## 17 (796900, 323500)    0         0    17     103   2   2
## 18 (713000, 311000)    0         0    18  272785   2   2
## 19 (607000, 295100)    0         0    19  256114   2   5
## 20 (413900, 308000)    0         0    20      60   2   2
##                            NAZ PAS NAW IATA ICAO
## 1                   Goszczanów   1   3  -99   EP
## 2           Kaliska-Międzychód   1   3  -99   EP
## 3           Kamieńsk-Orla Góra   1   1  -99   EP
## 4                Lipki Wielkie   2   3  -99   EP
## 5             Poznań-Kobylnica   2   3  -99 EPPK
## 6                      Rzeszów   2   3  -99 EPRJ
## 7                 Słupsk-Krępa   3   3  -99 EPSR
## 8               Szczecin-Dąbie   3   3  -99 EPSD
## 9          Wrocław-Mirosławice   2   3  -99   EP
## 10 Żar-Miedzybrodzie Żywieckie   2   3  -99 EPZR
## 11                       Chełm -99 -99  -99  -99
## 12                    Pakosław -99 -99  -99  -99
## 13                       Sanok -99 -99  -99  -99
## 14                      Sławka -99 -99  -99  -99
## 15          Aleksandrów Łódzki   1   3  -99  -99
## 16     Konin-Kazimierz Biskupi   2   1  -99 EPKB
## 17                Zamość-Mokre   3   3  -99 EPZA
## 18         Stalowa Wola-Turbia   2   3  QXQ EPST
## 19                     Pińczów   2   3  -99 EPPC
## 20      Opole-Polska Nowa Wieś   3   3  QPM EPOP
##                                          ZARZ WYS POLYGONID SCALE ANGLE
## 1                                 Jerzy Dulas 136         0     1     0
## 2                            Adam Smorawiński  85         0     1     0
## 3                                     Aviaeco 310         0     1     0
## 4  Nadleśnictwo Karwin Z Siedzibą W Drezdenku  45         0     1     0
## 5                          Aeroklub Poznański  86         0     1     0
## 6           Port Lotniczy/Aeroklub Rzeszowski 211         0     1     0
## 7                            Aeroklub Słupski  75         0     1     0
## 8                        Aeroklub Szczeciński   1         0     1     0
## 9                        Aeroklub Dolnośląski 153         0     1     0
## 10                   Górska Szkoła Szybowcowa 385         0     1     0
## 11                                        -98 -99         0     1     0
## 12                                        -98 -99         0     1     0
## 13                                        -98 -99         0     1     0
## 14                                        -98 -99         0     1     0
## 15                    Stowarzyszenie Lotnicze 190         0     1     0
## 16                          Aeroklub Koniński 110         0     1     0
## 17                  Aeroklub Ziemi Zamojskiej 225         0     1     0
## 18                     Aeroklub Stalowowolski 150         0     1     0
## 19                        Aeroklub Pińczowski 186         0     1     0
## 20                           Aeroklub Opolski 188         0     1     0

summary(wektor_punkt)

## Object of class SpatialPointsDataFrame
## Coordinates:
##              min    max
## coords.x1 192548 816191
## coords.x2 177732 756615
## Is projected: TRUE 
## proj4string :
## [+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000
## +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Number of points: 120
## Data attributes:
##       AREA     PERIMETER     LOTN_          LOTN_ID            LOT      
##  Min.   :0   Min.   :0   Min.   :  1.0   Min.   :     1   Min.   :1.00  
##  1st Qu.:0   1st Qu.:0   1st Qu.: 30.8   1st Qu.:    37   1st Qu.:2.00  
##  Median :0   Median :0   Median : 60.5   Median :    72   Median :2.00  
##  Mean   :0   Mean   :0   Mean   : 60.5   Mean   : 35890   Mean   :1.93  
##  3rd Qu.:0   3rd Qu.:0   3rd Qu.: 90.2   3rd Qu.: 54046   3rd Qu.:2.00  
##  Max.   :0   Max.   :0   Max.   :120.0   Max.   :306196   Max.   :2.00  
##                                                                         
##       TYP                                  NAZ           PAS       
##  Min.   :-99.0   Aleksandrów Łódzki          :  1   Min.   :-99.0  
##  1st Qu.:-99.0   Bagicz                      :  1   1st Qu.:  1.0  
##  Median :  2.0   Biała Podlaska              :  1   Median :  1.0  
##  Mean   :-25.3   Białystok-Krywlany          :  1   Mean   :-20.7  
##  3rd Qu.:  3.0   Bielsko-Biała-Aleksandrowice:  1   3rd Qu.:  2.0  
##  Max.   :  7.0   Borne Sulinowo              :  1   Max.   :  7.0  
##                  (Other)                     :114                  
##       NAW             IATA         ICAO                          ZARZ   
##  Min.   :-99.0   -99    :89   -99    :33   -98                     :49  
##  1st Qu.:-99.0   BZG    : 1   EP     : 8   Port Lotniczy           :10  
##  Median :  2.0   CZW    : 1   EPKP   : 2   Aeroklub Krakowski      : 2  
##  Mean   :-36.5   GDN    : 1   EPMM   : 2   Adam Smorawiński        : 1  
##  3rd Qu.:  3.0   IEG    : 1   EPPI   : 2   Aeroklub Białostocki    : 1  
##  Max.   :  3.0   KRK    : 1   EPPK   : 2   Aeroklub Bielsko-Bialski: 1  
##                  (Other):26   (Other):71   (Other)                 :56  
##       WYS          POLYGONID     SCALE           ANGLE  
##  Min.   :-99.0   Min.   :0   Min.   :0.000   Min.   :0  
##  1st Qu.:-99.0   1st Qu.:0   1st Qu.:1.000   1st Qu.:0  
##  Median : 85.0   Median :0   Median :1.000   Median :0  
##  Mean   : 67.8   Mean   :0   Mean   :0.992   Mean   :0  
##  3rd Qu.:156.8   3rd Qu.:0   3rd Qu.:1.000   3rd Qu.:0  
##  Max.   :628.0   Max.   :0   Max.   :1.000   Max.   :0  
## 

plot(wektor_punkt)

W podobny sposób można wczytać inne dane wektorowe - linie i poligony.

wektor_linia <- readOGR(dsn = 'dane/bdo250', layer = 'koleje_arc', encoding = 'Windows-1250')

## OGR data source with driver: ESRI Shapefile 
## Source: "dane/bdo250", layer: "koleje_arc"
## with 7057 features and 14 fields
## Feature type: wkbLineString with 2 dimensions

wektor_poligon <- readOGR(dsn = 'dane/bdo250', layer = 'admin_region_teryt_woj', encoding = 'Windows-1250')

## OGR data source with driver: ESRI Shapefile 
## Source: "dane/bdo250", layer: "admin_region_teryt_woj"
## with 16 features and 4 fields
## Feature type: wkbPolygon with 2 dimensions

Zapisywanie danych wektorowych w pakiecie rgdal odbywa się za pomocą funkcji writeOGR. W zależności od wybranego sterownika (a więc i formatu zapisu) ten proces wygląda trochę inaczej.

W przypadku zapisu pliku wektorowego do formatu shapefile należy zdefiniować nazwę obiektu w R, docelowy folder, docelową nazwę pliku oraz sterownik do obsługi formatu shapefile.

writeOGR(wektor_linia, 'dane/nowe', 'kolej', driver='ESRI Shapefile')

Obiekty rastrowe można wczytywać do R jako obiekty klasy SpatialGrid z wykorzystaniem funkcji readGDAL z pakietu rgdal.

library('rgdal')
r <- readGDAL('dane/srtm/srtm_40_02.tif')
summary(r)

Alternatywnie, pakiet raster pozwala za pomocą funkcji raster zarówno na tworzenie, jak i wczytywanie obiektów w formacie rastrowym. Wczytywanie może odbywać się z obiektów klasy macierz, 'image', Raster*, Spatial* i innych.

library('raster')
r <- raster('dane/srtm/srtm_40_02.tif')
r

## class       : RasterLayer 
## dimensions  : 6001, 6001, 36012001  (nrow, ncol, ncell)
## resolution  : 0.0008333, 0.0008333  (x, y)
## extent      : 15, 20, 50, 55  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/jn/Dropbox/występnienia/AdapteR/dane/srtm/srtm_40_02.tif 
## names       : srtm_40_02 
## values      : -32768, 32767  (min, max)

plot(r)

Do zapisu obiektu rastrowego służy funkcja writeRaster.

writeRaster(r, "dane/dem", "GTiff")

Istnieje wiele możliwych formatów zapisu danych, które można sprawdzić korzystając z funkcji writeFormats().

Najprostszym sposobem nadawania i zmiany układu współrzędnych w R jest wykorzystanie listy EPSG. Jest ona zawarta w bibliotece PROJ.4, dostępnej w R za pomocą pakietu rgdal.

Aby sprawdzić jaki układ współrzędnych jest nadany obiektowi przestrzennemu należy skorzystać z funkcji proj4string. W przypadku braku układu otrzymamy informację - NA.

proj4string(wektor_punkt)

## [1] "+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

proj4string(r)

## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Gdy nasz obiekt przestrzenny nie ma nadanego układu współrzędnych, lub też nadany układ jest błędny - należy to poprawić za pomocą funkcji CRS:

proj4string(wektor_punkt) <- CRS('+init=epsg:2180')
proj4string(r) <- CRS('+init=epsg:4326')

Przypisanie obiektowi układu współrzędnych nie wpływa na zawarte w nim wartości. Aby przetworzyć układ współrzednych w obiektach wektorowych należy skorzystać z funkcji spTransform:

proj4string(wektor_poligon)

## [1] "+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

wektor_poligon_WGS84 <- spTransform(wektor_poligon, CRS('+init=epsg:4326'))
proj4string(wektor_poligon_WGS84)

## [1] "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Do zmiany układu współrzędnych w obiektach rastrowych służy funkcja projectRaster. Należy podać jej nazwę naszego obiektu rastrowego oraz definicję układu współrzędnych w formacie PROJ.4.

r_1992 <- projectRaster(r, crs='+init=epsg:2180')
proj4string(r_1992)

## [1] "+init=epsg:2180 +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"

R pozwala na wizualizację danych przestrzennych za pomocą wielu pakietów, oferujących szerokie możliwości. Kilka najpopularniejszych z nich to:

• Pakiet sp - zawiera on metody tworzenia wizualizacji oparte zarówno o podstawowe możliwości R (plot, points, lines, image, itd.), jak i metody oparte o pakiet lattice (funkcja spplot).
• Pakiet ggplot2 - posiada złożone możliwości tworzenia wykresów; po odpowiednim przetworzeniu potrafi również przedstawiać dane przestrzenne.
• Pakiet rasterVis zawierający, między innymi, funkcję levelplot służącą wizualizacji danych rastrowych.
• Pakiet ggmap - oparty o ggplot2 pakiet pozwalający na tworzenie wizualizacji, w których warstwą podkładową są mapy z Google Maps czy OpenStreetMap.

Większość obiektów klasy Spatial* pozwala na ich wyświetlanie za pomocą funkcji plot. W przypadku, gdy dane nie mają nadanego układu współrzędnych - wielkość jednostek w osi x jest równa wielkości w osi y.

plot(wektor_punkt)
title('lotniska')

plot(wektor_linia)
title('sieć kolejowa')

plot(wektor_poligon)
title('województwa')

Możliwe jest również łączenie obiektów przestrzennych z wykorzystaniem parametru add.

plot(wektor_poligon)
plot(wektor_punkt, add=TRUE)

Zamiast osi możliwe jest umieszczenie podziałki liniowej. Dodatkowo, możliwe jest dodanie strzałki północy czy też ramki mapy.

plot(wektor_poligon, axes = FALSE)
SpatialPolygonsRescale(layout.scale.bar(), offset = c(180000, 150000), scale = 250000, 
    fill = c("transparent", "black"), plot.grid = F)  # dodanie podziałki liniowej
text(x = 180000, y = 180000, "0")
text(x = 430000, y = 180000, "250 km")
SpatialPolygonsRescale(layout.north.arrow(1), offset = c(200000, 700000), scale = 100000, 
    plot.grid = F)  # dodanie strzałki północy
box()  # dodanie ramki danych

Za pomocą kilku linii kodu można również uzyskać legendę mapy.

plot(wektor_poligon, axes=TRUE)
plot(wektor_punkt, add=TRUE, pch=18, col='blue')
legend("bottomleft", legend=c("Lotniska"), title="Legenda", bty="n", pch=18, col=c("blue"))

Zbiór danych dotyczących lotnisk zawiera, między innymi inforamację o liczbie pasów startowych. Można ją wyświetlić definiując odpowiednie grupy w funkcji legend.

wektor_punkt$PAS1 <- as.factor(wektor_punkt$PAS)
wektor_punkt$PAS[wektor_punkt$PAS == -99] <- NA
plot(wektor_poligon, axes = TRUE)
plot(wektor_punkt, add = TRUE, pch = 18, col = wektor_punkt$PAS)
cols <- 1:nlevels(wektor_punkt$PAS1)
legend(x = 150000, y = 250000, legend = c("brak danych", "1", "2", "3", "4", 
    "5", "6", "7"), title = "Liczba pasów startowych: ", bty = "n", pch = 18, 
    col = cols, cex = 0.8, ncol = 2)
title("Lotniska w Polsce")

spplot jest podstawową funkcją prezentacji danych przestrzennych w oparciu o pakiet lattice. Domyślnie tworzy ona szereg wizualizacji zgodnych z liczbą zmiennych w pliku Spatial*. Aby tego uniknąć, należy najpierw zadeklarować, która zmienna jest obiektem zainteresowania.

spplot(wektor_poligon['POP'])

Funkcja spplot wymaga podania jednego parametru, sp.layout, określającego dodatkowe elementy. Może to być skala liniowa, strzałka północy czy też argumenty modyfikujące prezentowane obiekty punktowe, liniowe czy poligonowe. Dostępne argumenty modyfikujące wygląd obiektów można znaleźć poprzez plik pomocy funkcji par - ?par.

strzałka <- list("SpatialPolygonsRescale", layout.north.arrow(), offset = c(200000, 
    700000), scale = 70000)
podziałka <- list("SpatialPolygonsRescale", layout.scale.bar(), offset = c(190000, 
    150000), scale = 250000, fill = c("transparent", "black"))
podziałka_tekst1 <- list("sp.text", c(190000, 180000), "0")
podziałka_tekst2 <- list("sp.text", c(440000, 180000), "250 km")
layout <- list(strzałka, podziałka, podziałka_tekst1, podziałka_tekst2)
spplot(wektor_poligon["POP"], sp.layout = layout)

Dane o podobnych wartościach można grupować. W poniższym przykładzie za pomocą funkcji brewer.pal z pakietu RColorBrewer tworzona jest nowa paleta kolorystyczna, a następnie określany jest podział wartości z wykorzystaniem funkcji seq. W funkcji spplot nadany jest obiekt do wizualizacji, wybrany podział danych (parametr at), oraz paleta kolorystyczna (parametr col.regions).

library('RColorBrewer')
kolor <- brewer.pal(6, 'Oranges')
podział <- seq(0, 6000000, 1000000)
spplot(wektor_poligon, 'POP', colorkey = list(labels = list(podział)), at = podział, col.regions=kolor)

Zamiast ręcznego określania podziału kolorystycznego, funkcja spplot wraz z funkcją classIntervals z pakietu classInt pozwala na zmianę sposobu klasyfikacji danych. Możliwy jest wybór jednej z kilku rodzajów podziału - 'fixed', 'sd', 'equal', 'pretty', 'quantile', 'kmeans', 'hclust', 'bclust', 'fisher' oraz 'jenks'. Wynik funkcji classIntervals może być też stosowany do tworzenia wizualizacji przy pomocy funkcji plot.

library('classInt')
kolor <- brewer.pal(8, 'Blues')
podział.qt <- classIntervals(wektor_poligon$POP, n = 8, style = 'quantile')
podział.eq <- classIntervals(wektor_poligon$POP, n = 8, style = 'equal')

spplot(wektor_poligon, 'POP', at=podział.qt$brks+1, col.regions=kolor)

spplot(wektor_poligon, 'POP', at=podział.eq$brks+1, col.regions=kolor)

Dane rastrowe można przedstawiać w R zarówno za pomocą funkcji plot, image czy spplot. Alternatywnie istnieje również pakiet rasterVis, który zawiera zestaw możliwości pozwalających na wizualizację danych rastrowych. Jego podstawowa funkcja to levelplot.

library('rasterVis')
levelplot(r)

Funkcja levelplot pozwala, między innymi, na dodanie izolini za pomocą parametru contour.

levelplot(r, contour = TRUE)

W łatwy sposób można też modyfikować skalę kolorystyczną. Wystarczy stworzyć wektor numeryczny zawierający kolejne przedziały, a następnie nadać go funkcji levelplot przy pomocy parametru at.

podział = c(0, 100, 200, 300, 400, 500, 700, 1000, 2000, 3000)
levelplot(r, contour = TRUE, margin = FALSE, at = podział)

levelplot posiada również możliwość łatwej zmiany kolorystyki przy pomocy parametru par.settings. Istnieje do wyboru kilka dostępnych schematów, jak również jest możliwość stworzenia własnych.

levelplot(r, contour = TRUE, margin = FALSE, at = podział, par.settings = GrTheme)

Dane wektorowe mogą zawierać ramy danych - wówczas możliwe jest stosowanie funkcji identycznych lub zbliżonych do podstawowych funkcji R. Przykładowo, wydzielenie z danych punktowych tylko tych lotnisk, które mają więcej niż 3 pasy oraz wydzielenie z danych poligonowych województw Wielkopolskiego i Lubuskiego.

lotniska3 <- subset(wektor_punkt, PAS>3)
wielkopolska <- subset(wektor_poligon, NAZ=='WOJ. WIELKOPOLSKIE')
lubuskie <- subset(wektor_poligon, NAZ=='WOJ. LUBUSKIE')

Pakiet rgeos jest interfejsem R do bibliotek GEOS (Geometry Engine - Open Source). Zawiera zbiór funkcji do operacji topologicznych na geometriach. Obejmują one między innymi obliczenia powierzchni (funkcja gArea) czy długości (funkcja gLength). Możliwości tego pakietu są jednak znacznie szersze.

library('rgeos')

## rgeos version: 0.3-8, (SVN revision 460)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Polygon checking: TRUE

gArea(wielkopolska) #w m2

## [1] 29805838235

gLength(wektor_linia) #w m

## [1] 25676871

Funkcja over służy dodaniu obiektowi x (w przykładzie wektor_punkt) atrybutów z obiektu y (wektor_poligon).

## OGR data source with driver: ESRI Shapefile 
## Source: "dane/bdo250", layer: "lotn_point"
## with 120 features and 16 fields
## Feature type: wkbPoint with 2 dimensions

wektor_punkt2 <- wektor_punkt
wektor_punkt2@data <- over(wektor_punkt, wektor_poligon)

Zmiejszenie rozdzielczości obiektu rastrowego można wykonać poprzez funkcję aggregate. Należy w niech zdefiniować parametr dact, mówiący o liczbie oczek siatki w każdym kierunku, które będą brane do zmniejszania rozdzielczości obiektu oraz parametr fun, mówiący o wykorzystanej funkcji (np. średnia).

r

## class       : RasterLayer 
## dimensions  : 6001, 6001, 36012001  (nrow, ncol, ncell)
## resolution  : 0.0008333, 0.0008333  (x, y)
## extent      : 15, 20, 50, 55  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/jn/Dropbox/występnienia/AdapteR/dane/srtm/srtm_40_02.tif 
## names       : srtm_40_02 
## values      : -32768, 32767  (min, max)

r_roz <- aggregate(r, fact = 16, fun = mean)
r_roz

## class       : RasterLayer 
## dimensions  : 376, 376, 141376  (nrow, ncol, ncell)
## resolution  : 0.01333, 0.01333  (x, y)
## extent      : 15, 20.01, 49.99, 55  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /tmp/R_raster_jn/raster_tmp_2014-10-18_160103_21724_05141.grd 
## names       : srtm_40_02 
## values      : -6, 1473  (min, max)

Odczytanie wartości z obiektu rastrowego może nastąpić za pomocą funkcji getValues.

getValues(r_roz)

##  [1]  -4.000  -4.000  -4.000  -4.000  -4.000  -3.535  -1.512   1.406
##  [9]   5.875   7.992   8.523   7.742   7.164   4.891  15.551  22.262
## [17]  31.148  36.949  39.297  45.930  51.145  51.941  54.723  35.094
## [25]  38.691  48.648  80.621  96.367  96.453  96.562  92.527  86.355
## [33]  78.281  71.242  64.344  91.562 106.961 107.027 109.965 109.801
## [41] 105.992  96.461  88.125  74.520  59.895  52.695  51.527  65.934
## [49]  74.246  89.332 105.129  97.707 109.680 105.832  66.660  53.344
## [57]  17.984  22.922  52.836  60.164  70.875  74.008  72.668  83.102
## [65]  75.387  74.055  75.156  86.582  98.230  73.297  54.512  38.422
## [73]  38.402  41.035  33.703  26.520  23.504  22.043  18.629  14.707
## [81]  13.098  13.287  18.213  11.250      NA

Przetwarzanie danych rastrowych jest bardzo zbliżone do wykonywania operacji na matrycach. Funkcje, operacje warunkowe czy arytmetyczne działają w ten sam sposób.

raster_2 <- r*2
raster_3 <- r + raster_2

W momencie, gdy konieczne jest przycięcie warstwy rastrowej do warstwy poligonowej należy skorzystać z funkcji crop lub mask. Wynikiem działania funkcji crop jest raster o takim samym zasięgu jak obiekt poligonowy. Funkcja mask produkuje obiekt dokładnie przycięty do granic warstwy wektorowej, jednakże jest to zadanie bardziej czasochłonne.

wielkopolska_wgs84 <- spTransform(wielkopolska, CRS('+proj=longlat +datum=WGS84'))
wielkopolska_wgs84_srtm <- crop(r, wielkopolska_wgs84)
levelplot(wielkopolska_wgs84_srtm, margin=FALSE, par.settings='BuRdTheme')

wielkopolska_wgs84_srtm2 <- mask(r, wielkopolska_wgs84)
levelplot(wielkopolska_wgs84_srtm2, margin=FALSE, par.settings='BuRdTheme')

wielkopolska_wgs84_srtm <- crop(r, wielkopolska_wgs84)
wielkopolska_wgs84_srtm2 <- mask(wielkopolska_wgs84_srtm, wielkopolska_wgs84)
levelplot(wielkopolska_wgs84_srtm2, margin=FALSE, par.settings='BuRdTheme')

Zbiór narzędzi statystycznych uwzględniających w analizie danych ich przestrzenną lokalizację i czasowy zakres, a opartych o teorię funkcji losowych

Funkcje geostatystyki:

• Identyfikacja i modelowanie struktury przestrzennej/czasowej/czasoprzestrzennej zjawiska
• Estymacja - szacowanie wartości parametru w nieopróbowanym miejscu i/albo momencie czasu
• Symulacja- generowanie równie prawdopodobnych obrazów (realizacji), które honorują wyniki pomiarów i strukturę przestrzenną/czasową zjawiska,
• Optymalizacja próbkowania/sieci pomiarowej minimalizująca koszty przy maksymalizacji ilości i dokładności informacji

W R istnieje szereg pakietów do geostatystyki. Najpopularniejsze to gstat oraz geoR.

library('raster')
library('rasterVis')
library('gstat')
library('geostatsp')

## Loading required package: Matrix

data('swissRain')

cols <- brewer.pal(6, 'Blues')
opad <- spplot(swissRain['rain'], cex=1.5, pch=20, colorkey=TRUE, col.regions=cols, cuts=length(cols)-1)
opad

swissDEM <- crop(swissAltitude, swissBorder)
swissDEM <- mask(swissDEM, swissBorder)
p <- levelplot(swissDEM, main='elevation', par.settings=rasterTheme(region=terrain.colors(5)))
p <- p + opad
p <- p + layer(sp.polygons(swissBorder, lwd=2))
p

hist(swissRain$rain)

rain.svar <- variogram(rain ~ 1, swissRain)
plot(rain.svar)

przekatna <- sqrt(diff(swissBorder@bbox["x",])^2 + 
                          diff(swissBorder@bbox["y",])^2)
rain.ivgm <- vgm(nugget=0, model="Sph", range=przekatna/4,
                 psill=var(swissRain$rain))
rain.vgm <- fit.variogram(rain.svar, model=rain.ivgm)
plot(rain.svar, rain.vgm)

g <- as(swissDEM, 'SpatialGridDataFrame')
OK <- gstat(id = "ordinary_kriging", formula = rain ~ 1, 
                data = swissRain, model = rain.vgm)
prediction <- predict(OK, newdata = g, debug.level = 0)
g$OK_pred <- prediction$ordinary_kriging.pred
g$OK_se <- sqrt(prediction$ordinary_kriging.var)

library('classInt')
cols <- brewer.pal(6, 'Blues')
brks <- classIntervals(swissRain$rain, n = 6, style = 'quantile')
cols.se <- brewer.pal(5, 'Reds')
brks.se <- seq(0, 12.5, 2.5)
library('gridExtra')

## Loading required package: grid

p_OKp <- spplot(g, "OK_pred",  at = brks$brks, col.regions = cols)
p_OKs <- spplot(g, "OK_se",  at = brks.se, col.regions = cols.se)

grid.arrange(p_OKp, p_OKs, ncol=2)

• Statystyka punktów
• Ekonometria przestrzenna
• Analizy czasoprzestrzenne
• Inne

• Bivand, R.S., Pebesma, E.J., Gómez-Rubio, V. (2013). Applied Spatial Data Analysis with R. Springer
• http://www.spatstat.org/
• http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf
• http://oscarperpinan.github.io/rastervis/
• http://www.maths.lancs.ac.uk/~rowlings/Teaching/Sheffield2013/index.html
• http://spatial-analyst.net/book/
• http://cran.r-project.org/web/views/Spatial.html
• http://srtm.csi.cgiar.org
• http://www.codgik.gov.pl/index.php/darmowe-dane/