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)