## Acumulado Anual médio em cinco anos (2014-2018) setwd("C:/Users/eliasmedeiros/Dropbox/UFGD/Projeto de Pesquisa/PIBIC_2019_2020/IC_Hugo/mapa") require(rgdal) require(maps) require(sp) require(gstat) require(sp) require(grid) ctba <- readOGR("27MEE250GC_SIR.shp") plot(ctba) pb_meso <- spTransform(ctba, CRS("+proj=utm +zone=24 +south +ellps=intl +towgs84=-205.57,168.77,-4.12,0,0,0,0 +units=km +no_defs")) rv = list("sp.polygons", pb_meso, fill = "transparent", lwd=2) Dados_Alagoas <- read.csv2("Dados_Al_Total_Mes_Ano.csv") str(Dados_Alagoas) head(Dados_Alagoas) Dados_Alagoas <- subset(Dados_Alagoas, Ano>2013) qt = function(x) quantile(x, c(0.025, 0.25, 0.50, 0.75, 0.975)) aggregate(Prec ~ Mes, FUN = qt, data = Dados_Alagoas) dd = aggregate(Prec ~ Ano + Mes + Lat + Long + EstacaoCodigo, data = Dados_Alagoas, FUN = sum) dd2 = aggregate(Prec ~ Mes + Lat + Long + EstacaoCodigo, data = dd, FUN = mean) require(ggplot2) p1 <- ggplot(data = dd, aes(x = factor(Mes), y = Prec))+ #facet_wrap(~EstacaoCodigo, scales = "free")+ geom_boxplot(fill = '#386cb0', alpha=0.5) + labs(y = "Precipitação mensal (mm)", x = "Meses")+ theme_bw() theme_set(theme_cowplot(font_size=4)) # reduce default font size plist <- list(p1) grid.arrange(grobs = plist, ncol = 1) ## display plot ggsave(file = "boxplotICHugo.jpg", arrangeGrob(grobs = plist, ncol = 1), width = 16, height = 8, units = "cm", dpi = 300) dd1 = subset(dd2, Mes==5) coordinates(dd1) <- ~ Long + Lat proj4string(dd1) <- CRS("+proj=longlat +ellps=WGS84 +towgs84=-205.57,168.77,-4.12,0,0,0,0") ozone.UTM <- spTransform(dd1, CRS("+proj=utm +zone=24 +south +ellps=intl +towgs84=-205.57,168.77,-4.12,0,0,0,0 +units=km +no_defs")) bm = colorRampPalette(c('#d73027', '#f46d43', '#fdae61', '#fee090', '#e0f3f8', '#abd9e9', '#74add1', '#4575b4')) rv = list("sp.polygons", pb_meso, fill = "transparent", lwd=2, col="black", first=FALSE) l2 = list("SpatialPolygonsRescale", layout.north.arrow(), offset = c(600,8838), scale = 20) l3 = list("SpatialPolygonsRescale", layout.scale.bar(), offset = c(840,8855), scale = 50, fill=c("transparent","black")) l4 = list("sp.text", c(840,8850), "0") l5 = list("sp.text", c(840+50,8850), "50 km") p1 <- spplot(ozone.UTM['Prec'], cex=2, col.regions=bm(15), xlim=c(580,930), ylim = c(8835,9027), sp.layout = list(l2, l3, l4, l5, rv), cuts=quantile(ozone.UTM['Prec']$Prec), key.space="right", colorkey=F, region = TRUE) png(filename = "Media_des.png", res=300, width = 2680, height = 1280) plot(p1) dev.off() # Análise do mes de Maio ozoneDF <- data.frame('prec'=ozone.UTM$Prec, 'longx' = ozone.UTM@coords[,1], 'longx2' = ozone.UTM@coords[,1]^2, 'laty' = ozone.UTM@coords[,2]) mod <- lm(prec ~ ., data=ozoneDF) summary(mod) mod <- step(mod) summary(mod) ozone.UTM$estimado <- fitted(mod) ozoneDF$estimado <- fitted(mod) R2 <- cor(ozone.UTM$Prec, ozone.UTM$estimado)^2 cat(R2) #ozoneDF$res <- ozoneDF$prec - ozoneDF$estimado ozoneDF$res <- ozoneDF$prec - ozoneDF$estimado ozone.UTM$res <- ozoneDF$prec - ozoneDF$estimado require(ggplot2) ggplot(ozoneDF, aes(x=res)) + geom_histogram(aes(y=..density..), binwidth=150, colour="black", fill="white") + labs(x="Residuals", y="Density") + #xlim(c(-325,750)) + theme(axis.text.x=element_text(size=16)) + theme(axis.text.y=element_text(size=16)) + theme(axis.title.x=element_text(size=20)) + theme(axis.title.y=element_text(size=20)) + geom_density(alpha=.2, fill="#FF6666") lzn.vgm <- variogram(res~1, ozone.UTM, cutoff=110, width = 20) plot(lzn.vgm) Sph.fit = fit.variogram(lzn.vgm, model = vgm(psill = 2000, model = "Sph", range = 60, nugget = 500)) Sph.fit plot(lzn.vgm, Sph.fit, xlab = "Distância (km)", ylab="Variograma", pch=15, col='black') Exp.fit = fit.variogram(lzn.vgm, model = vgm(psill = 200, model = "Exp", range = 60, nugget = 50)) Exp.fit plot(lzn.vgm, Exp.fit, xlab = "Distância (km)", ylab="Variograma", pch=15, col='black', lwd=2) Gau.fit = fit.variogram(lzn.vgm, model = vgm(psill = 200, model = "Gau", range = 60, nugget = 50)) Gau.fit plot(lzn.vgm, Gau.fit, xlab = "Distância (km)", ylab="Variograma", pch=15, col='black', add=T) x <- krige.cv(res~1, ozone.UTM, Sph.fit, nmax = 10, nfold=2) ozone.UTM$vd <- fitted(mod) + x$var1.pred cor(ozone.UTM$vd, ozone.UTM$Prec) # cross validation require(hydroGOF) x <- krige.cv(Prec~I(Long^2), ozone.UTM, Gau.fit, nmax = 10) summary(lm(x$var1.pred ~ ozone.UTM$Prec)) gof(x$var1.pred, ozone.UTM$Prec, digits=4) Fitted <- data.frame(dist = seq(0.01, max(lzn.vgm$dist), length = 10001)) Fitted$gamma <- variogramLine(Exp.fit, dist_vector = Fitted$dist)$gamma require(reshape) Modeled <- melt(Fitted, id.vars = "dist", measure.vars = c("gamma")) Modeled$modelos <- factor(Modeled$variable, levels = c("gamma"), labels = c("Exponencial")) p1 = ggplot(lzn.vgm, aes(x = dist, y = gamma)) + geom_point() + geom_line(data = Modeled, aes(x = dist, y = value, colour = "darkblue"), size=1) + labs(x = "Distância (km)", y = "Variograma", colour = 'Modelo', subtitle = "\n Pepita = 0 \n Patamar = 3389 \n Alcance = 35 km") + theme(legend.position = c(0.75,0.15)) + scale_color_manual(values = 'darkblue', labels = "Exponencial") theme_set(theme_cowplot(font_size=10)) # reduce default font size plist <- list(p1) grid.arrange(grobs = plist, ncol = 1) ## display plot ggsave(file = "vgms.jpg", arrangeGrob(grobs = plist, ncol = 1), width = 14, height = 8, units = "cm") ptsreg <- spsample(pb_meso, 50000, type = "regular") ddp <- data.frame('Long'=ptsreg@coords[,1], 'Lat' = ptsreg@coords[,2]) coordinates(ddp) <- ~ Long + Lat ddp@proj4string <- ozone.UTM@proj4string Ord.kriged = krige(Prec ~ I(Long^2), ozone.UTM, ddp, model = Exp.fit) y = Ord.kriged['var1.pred']$var1.pred spplot(Ord.kriged['var1.pred'], col.regions=bm(4), xlim=c(580,930), ylim = c(8835,9027), sp.layout = list(l2, l3, l4, l5, rv), colorkey = list( right = list( fun = draw.colorkey, args = list( key = list( at = quantile(y), # colour breaks col = bm(11), # colours labels = list( at = quantile(y), labels = c("0 mm", "88 mm", "176 mm", "236 mm", "363 mm") ) ) ) ) )) y = Ord.kriged['var1.pred']$var1.pred p1 <- spplot(Ord.kriged['var1.pred'], cex=0.5, col.regions=bm(15), xlim=c(580,930), ylim = c(8835,9027), sp.layout = list(l2, l3, l4, l5, rv), cuts=seq(0,375,25), key.space="right", colorkey=T, region = TRUE) png(filename = "krige_maio.png", res=300, width = 2680, height = 1280) plot(p1) dev.off() Ord.kriged$ErroPred <- sqrt(Ord.kriged$var1.var) spplot(Ord.kriged['ErroPred'], col.regions=bm(30), xlim=c(580,930), ylim = c(8835,9027), sp.layout = list(l2, l3, l4, l5, rv), cuts=seq(0,60,5), key.space="right", colorkey=T, region = TRUE) # ##################### ## calculando as probabilidades: abaixo de 50 mm (<= 50) Ord.kriged$A <- pnorm(25.0, Ord.kriged$var1.pred, sqrt(Ord.kriged$var1.var), lower.tail=TRUE) Ord.kriged$B <- pnorm(50.0, Ord.kriged$var1.pred, sqrt(Ord.kriged$var1.var), lower.tail=TRUE) Ord.kriged$C <- pnorm(75.0, Ord.kriged$var1.pred, sqrt(Ord.kriged$var1.var), lower.tail=TRUE) Ord.kriged$D <- pnorm(100.0, Ord.kriged$var1.pred, sqrt(Ord.kriged$var1.var), lower.tail=TRUE) bmprobs <- colorRampPalette( c('#fff7ec', '#fee8c8','#fdd49e','#fdbb84', '#fc8d59', '#ef6548', '#d7301f', '#990000')) l2 = list("SpatialPolygonsRescale", layout.north.arrow(), offset = c(600,8838), scale = 30) l3 = list("SpatialPolygonsRescale", layout.scale.bar(), offset = c(840,8858), scale = 50, fill=c("transparent","black")) l4 = list("sp.text", c(840,8850), "0") l5 = list("sp.text", c(840+50,8850), "50 km") p1 <- spplot(Ord.kriged, c("C", "D", "A", "B"), cex=0.5, names.attr = c("P(Precipitação < 75 mm)", "P(Precipitação < 100 mm)", "P(Precipitação < 25 mm)", "P(Precipitação < 50 mm)"), col.regions=bmprobs(15), xlim=c(580,930), ylim = c(8835,9027), sp.layout = list(l2, l3, l4, l5, rv), cuts=seq(0, 1, 0.20), key.space="right", colorkey=T, region = TRUE) png(filename = "krige_probs.png", res=300, width = 2680, height = 1880) plot(p1) dev.off()