Jakub Nowosad
- Zakład Geoinformacji, WNGiG, UAM
- nowosad@amu.edu.pl
- http://nowosad.github.io
Dane
15.10.2014
Jakub Nowosad
Dane
"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
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:
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:
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)