# załadowanie bibliotek programistycznych
# w pierwszej kolejności najprawdopodobniej trzeba będzie je zainstalować
library("sp") # biblioteka do wizualizacji danych przestrzennych
library("raster") # biblioteka do obsługi plików rastrowych
library("rgdal") # obsługa rozszerzeń gdal
## rgdal: version: 1.1-3, (SVN revision 594)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
## Path to GDAL shared files: /usr/share/gdal/1.11
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.2-1
library("rgeos") # obsługa bibliotek geoprzestrzennych
## rgeos version: 0.3-12, (SVN revision 498)
## GEOS runtime version: 3.5.0-CAPI-1.9.0 r4084
## Linking to sp version: 1.2-0
## Polygon checking: TRUE
library("RColorBrewer") # automatyczne generowanie skal kolorystycznych
kolory = colorRampPalette(brewer.pal(7,"OrRd")[-1])(50)
sp.theme(set = TRUE, regions = list(col = kolory))
# Zastosujemy w obliczeniach prosty model dyspersji zanieczyszczen
# Calkowita koncentracja zanieczyszczeń w dystansie 'd' i dla kąta θ z punktowego źródła zanieczyszczeń
# dana jest poniższym wzorem:
#f=α*exp(−(d/β)2e−^(κ*cos(θ−ϕ)^2))
# gdzie:
# a - poziom zanieczyszczeń w punkcie emisji
# β - współczynnik odległości
# ϕ - dominujący kierunek smugi zanieczyszczeń
# κ - ekscentryczność (mimośródowość)
# zapiszmy teraz powyższy wzór w postaci funkcji programistycznej:
pfunc <- function(d2, ang, a, b, k, phi)
{ a * exp(-(d2*exp(-k*cos(ang-phi))^2/b^2)) }
# a następnie zastosujmy ją do obliczeń w regularnej siatce punktów
plume <- function(src,dst,a,b,k,phi)
{ src = coordinates(src)
dst = coordinates(dst)
d2 = (dst[,1]-src[,1])^2 + (dst[,2]-src[,2])^2
ang = atan2(dst[,2]-src[,2],dst[,1]-src[,1])
pfunc(d2,ang,a,b,k,phi) }
# modyfikacja ukladu wspolrzednych do postaci kompatybilnej z naszymi wzorami fizycznymi (zapisanymi w układzie SI)
toOSGB <- function(s){spTransform(s,CRS("+init=epsg:27700"))}
cumbria <- toOSGB(readOGR(dsn="cumbria.shp",layer="cumbria",verbose=FALSE)) # warstwa GIS: granice administracyjne
elektrownia <- toOSGB(readOGR(dsn="nuclear.shp",layer="nuclear",verbose=FALSE)) # warstwa GIS: elektrownia
miasta <- readOGR(dsn="ukgaz.shp",layer="ukgaz",verbose=FALSE) # warstwa GIS: miasta
# wyrysujmy zatem wczytane pliki:
plot(cumbria)
plot(elektrownia,add=TRUE,pch=19,cex=3.5,col="red") #Definiuje kolor oraz ksztalt punkt emitującego zanieczyszczenia.
plot(miasta,add=TRUE)
# tworzymy pusty raster o zasięgu odpowiadającym zasięgowi naszego poligonu
plevel <- raster(extent(cumbria),nrows=500,ncols=500)
# wykonajmy obliczenia w utworzonej siatce punktów:
plevel[] <- plume(src=elektrownia, dst=SpatialPoints(plevel), a=9000, b=400, k=5, phi=70*pi/180)
projection(plevel) <- proj4string(cumbria) # nadanie układu współrzędnych
plot(plevel,col=kolory) # wyrysowanie wyniku
plot(cumbria,add=TRUE, lty=4, lwd=2) # dodanie granic administracyjnych
plot(elektrownia,add=TRUE,pch=19,cex=3.5,col="blue")
plot(miasta,add=TRUE)
title("Rozprzestrzenianie chmury radioaktywnej")
#Wykaz miast zagrozonych zanieczyszczeniami powyzej zdefiniowanie wielkosci (tu: 6000).
# 1: ekstrakcja warstwy rastrowej do punktow:
miasta$exposure <- extract(plevel, miasta)
# 2: wybranie tylko wartosci wiekszych od wartosci progowej:
over6000 <- miasta$exposure>6000
# 3: utworzenie ramki danych
zagrozone_miasta <- data.frame(miasta@data[over6000,c("Name","pop","exposure")])
zagrozone_miasta
## Name pop exposure
## 3 Cockermouth 7551 7154.426
## 10 Dearham 1391 6125.387
## 13 Great Clifton 1101 6631.521
## 14 Dean 1077 7721.381
## 15 Brigham 1065 7129.519
## 18 Broughton Moor 726 6235.179
## 20 Bridekirk 636 6793.254
## 33 Bassenthwaite 412 6176.067
## 35 Papcastle 406 7078.424
## 36 Little Clifton 391 7087.436
## 37 Plumbland 367 6035.050
## 40 Gilcrux 303 6161.290
## 41 Embleton 297 7024.913
## 42 Blindcrake 287 6563.121
## 43 Lorton 250 7301.323
## 46 Loweswater 209 7701.603
## 47 Winscales 186 6876.589
## 48 Camerton 172 6322.226
## 50 Buttermere 127 6316.868
## 51 Setmurthy 108 6427.336
## 52 Wythop 51 6890.644
## 87 Egremont 6496 7923.451
## 89 Cleator Moor 6347 8186.493
## 90 Arlecdon 3678 8201.449
## 91 Frizington 2415 8306.656
## 92 Distington 2283 6526.419
## 93 St John Beckermet 1925 8619.762
## 96 Moresby 1280 6027.835
## 100 Lamplugh 763 8190.845
## 104 St Bridget Beckermet 385 8619.762
## 109 Ennerdale Bridge 240 8591.672
## 110 Haile 237 8901.306
## 112 Ponsonby 92 8825.822
# zadania do wykonania:
# Zadanie 1:
# Jaki kształt ma smuga zanieczyszczeń jeśli κ = 0 , a jaki przy wartościach dodatnich>
# Który z parametrów (w uproszczeniu) spośród parametrów κ oraz ϕ możemy traktować jako siłę i kierunek wiatru?
# Zadanie 2:
# Korzystając z pomocy do funkcji 'colorRampPalette' zmień stosowany schemat barwny
# Zadanie 3:
# Jakie wartości parametru tego parametru odpowiadają kątowi:
# 0, 45, 90, 180, 270 stopni? (należy traktować analogicznie jak kierunek wiatru)
# Zadanie 4:
# Wymień 3 miasta o populacji powyżej 300 mieszkańców, o największym zagrożeniu chmurą radioaktywną jeśli kierunek napływu mas
# powietrza zmieni się na zachodni