piotr_anomalie.R

library(sp)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(raster)

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

p95 <- read.asciigrid("~/Pobrane/p95.asc", as.image=FALSE, plot.image=TRUE)
obs <- read.asciigrid("~/Pobrane/obs.asc", as.image=FALSE, plot.image=TRUE)
# 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

wynik <- raster(p95) - raster(obs)
wynik[wynik<0] <- 0 # pozbywamy sie anomalii <0 i zastepujemy "0"
calosc <- mask(wynik, pl,inverse=T) # maskujemy wyniki tylko do granic administracyjnych PL
# ten proces moze chwile zajec bez zrownoleglania...
plot(calosc, main="anomalie", 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"

dane <- na.omit(data.frame(pow=cell_size[],anom=calosc[]))

# obliczenie anomalii 0-1
policz_anomalie <- function(klasa1=0, klasa2=10, interwal=0.5, dane=dane){
  
  brejki <- c(-Inf,seq(from=klasa1,to=klasa2, by=interwal),Inf)
  tmp <- NULL
    for (i in 1:(length(brejki)-1)){
    suma_powierzchni <- dplyr::filter(dane, anom>brejki[i], anom<=brejki[i+1]) %>% pull(pow) %>% sum()
    udzial_powierzchni <- suma_powierzchni/raster_area
    wynik <- data.frame(suma_powierzchni=suma_powierzchni, 
                        udzial_powierzchni=udzial_powierzchni, klasa1=brejki[i], klasa2=brejki[i+1])
    tmp <- rbind.data.frame(tmp, wynik)
  }
  return(tmp)
}

policz_anomalie(klasa1=0,klasa2=1, interwal = 0.5, dane = dane)

##   suma_powierzchni udzial_powierzchni klasa1 klasa2
## 1        18049.692         0.05818924   -Inf    0.0
## 2         9439.851         0.03043253    0.0    0.5
## 3         9490.850         0.03059694    0.5    1.0
## 4       273209.128         0.88078129    1.0    Inf

# albo zmieniajac nieco ustawienia funkcji:
policz_anomalie(klasa1=0,klasa2=10, interwal=1,dane = dane)

##    suma_powierzchni udzial_powierzchni klasa1 klasa2
## 1          18049.69         0.05818924   -Inf      0
## 2          18930.70         0.06102947      0      1
## 3          28794.93         0.09283013      1      2
## 4          35725.69         0.11517376      2      3
## 5          55304.99         0.17829419      3      4
## 6          59818.09         0.19284368      4      5
## 7          75145.68         0.24225732      5      6
## 8          18419.74         0.05938221      6      7
## 9              0.00         0.00000000      7      8
## 10             0.00         0.00000000      8      9
## 11             0.00         0.00000000      9     10
## 12             0.00         0.00000000     10    Inf

calosc[] <- round(calosc[]) # mozemy takze zaokraglic wyniki dla ladniejszego efektu, ale to bardziej do plotowania
plot(calosc, col=c( "white", rev(heat.colors(n=7))))
plot(pl, add=T)