library(sp)
library(raster)
p95 <- read.asciigrid("~/Pobrane/p95.asc", as.image=FALSE, plot.image=TRUE)
obs <- read.asciigrid("~/Pobrane/obs.asc", as.image=FALSE, plot.image=TRUE)
wynik <- raster(p95) - raster(obs)
plot(wynik, main="anomalie #1", col=terrain.colors(255))
wynik[wynik<0] <- 0 # pozbywamy sie anomalii <0
plot(wynik, main="anomalie - tylko dodatnie", col=terrain.colors(255)) #
# wczytanie szejpa z granicami administracyjnymi
# akurat w zasiegu mialem tylko taki dziwny z mapa Europy i "z dziura" dla PL
pl <- maptools::readShapeSpatial("/home/bartosz/R/skryptydoR/gis/dziurapl.shp") ;
## Warning: use rgdal::readOGR or sf::st_read
## Warning: use rgdal::readOGR or sf::st_read
calosc <- mask(wynik, pl,inverse=T) # maskujemy wyniki tylko do granic administracyjnych PL
# ten proces moze chwile zajec bez zrownoleglania...
# i dostajemy takie cos:
plot(calosc, col=terrain.colors(255))
calosc[] <- round(calosc[]) # mozemy takze zaokraglic wyniki dla ladniejszego efektu
plot(calosc, col=terrain.colors(255))
#####################################
#####################################
### jesli brac pod uwage rzeczywista powierzchnie gridow:
# 1. oblicz rozmiar wszystkich gridow z rastra [km2]
cell_size<-area(calosc, na.rm=TRUE, weights=FALSE)
## Warning in couldBeLonLat(out): CRS is NA. Assuming it is longitude/latitude
plot(cell_size, main="Powierzchnia grida w km^2")
raster_area <- sum(cell_size[], na.rm=T)
print(paste("obszar PL z rastra:",round(raster_area, digits=1),"km^2"))
## [1] "obszar PL z rastra: 310189.5 km^2"
# obliczenie anomalii 0-1
anomalie01 <- calosc
anomalie01[anomalie01 >0] <- 1
anomalie_w_km2 <- anomalie01*cell_size
sum(anomalie_w_km2[], na.rm=T) # powierzchnia z anomaliami dodatnimi
## [1] 282700
sum(anomalie_w_km2[], na.rm=T)/raster_area # udzial powierzchni z anomaliami dodatnimi
## [1] 0.9113782
# stara czesc ze statystykami...
# pobierzmy z rastra wartosci do wektora
a <- getValues(calosc)
summary(a)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0 2.0 4.0 3.7 5.0 7.0 550285
length(which(is.na(a))) # tyle gridow jest poza PL (tylko w ramach radociekawostki)
## [1] 550285
grids_pl <- length(a)-length(which(is.na(a))) # tyle gridow jest ogolnie w PL
grids_pos <- length(which(a>0)) # a tyle gridow ma dodatnie anomalie
grids_pos/grids_pl * 100 # ile obszaru (%) ma dodatnie anomalie
## [1] 91.53704
# a jesli chcesz dokladniejsze staty to masz tabelke
data.frame(klasy_anomalii=paste0("<=",names(table(a)/grids_pl)),
udzial_procentowy=paste0( round(table(a)/grids_pl,3) *100,"%")
)
## klasy_anomalii udzial_procentowy
## 1 <=0 8.5%
## 2 <=1 6.7%
## 3 <=2 11.2%
## 4 <=3 12%
## 5 <=4 20.3%
## 6 <=5 20.9%
## 7 <=6 19.5%
## 8 <=7 1%