Determinação de variáveis derivadas
Roilan Hernández Valdés
Lendo os dados observados
# all.vars <- readRDS("~/Dropbox/Dissertacao/EstSP.MG_2001-2014_FINAL.rds")
info <- readRDS("~/Dropbox/Dissertacao/1_EstSP.MG_2001-2014_info.rds")Praparação dos dados
basin_name <- "mogu"
# basin_name <- "pauliceia"Reprojeção do raster em coordenadas tipo Lat/Lon
# Carregando Modelo Digital de Elevação do Terreno (MDET)
wsdem <- raster(paste("../MODIS/wsdem_",basin_name,".asc",sep ="" ))
# wsdem <- raster("../MODIS/wsdem.asc")
# Corrigindo extent do wsdem (as coordenadas devem ser centradas)
fake.ext <- extent(wsdem)
fake.ext <- c(-(fake.ext[2]-fake.ext[1])/2,(fake.ext[2]-fake.ext[1])/2,
-(fake.ext[4]-fake.ext[3])/2,(fake.ext[4]-fake.ext[3])/2)
extent(wsdem) <- fake.ext
# Carregando as projeções do arquivo
# Lat/Lon
basin_lonlat <- raster(paste("../MODIS/mdet_",basin_name,".tif",sep ="" ))
plot(basin_lonlat,
main = "Domínio da Bacia e linhas centrais")
abline(v = mean(extent(basin_lonlat)[1:2]))
abline(h = mean(extent(basin_lonlat)[3:4]))
lonlat <- projection(basin_lonlat)
# Lambert Azimutal
basin_laea <- raster(paste("../MODIS/mdet_",basin_name,"_laea_int.tif",sep ="" ))
plot(basin_laea)
# abline(v = 0.0, h = 0.0)
laea <- projection(basin_laea)
# Atribuindo a projeção laea ao arquivo raster da bacia
projection(wsdem) <- laea
# Reprojectando a Lat/Lon
wsdem <- projectRaster(from = wsdem,crs = laea )
wsdem <- projectRaster(from = wsdem,crs = lonlat )
plot(wsdem/1000, main = "Elevação Digital do Terreno (km)")
# Criando Shapefile da bacia MoGu
wsdem_shp <- reclassify(wsdem , c(0,2000,0)) %>%
rasterToPolygons(dissolve = TRUE)## Loading required package: rgeos## rgeos version: 0.3-15, (SVN revision 515)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Linking to sp version: 1.2-1
## Polygon checking: TRUE##
## Attaching package: 'rgeos'## The following object is masked from 'package:Hmisc':
##
## translate plot(wsdem_shp, add = TRUE)
Criação do raster com localização das estações
# Raster com a mesma resolução que "wsdem"
r <- raster(resolution = res(wsdem))
# determinando cobertura espacial das estações INMET
extent.estacoes <- extent(min(info$lon), max(info$lon) , min(info$lat), max(info$lat) )
# recortando raster a área espacial das estações
r <- crop(r, extent.estacoes)
# rasterizando coordenadas das estações
p <- info[,c("lon","lat")]
coordinates(p) <- ~lon+lat
estacoes <- rasterize(p,r)
# domínio da bacia MoGu e localização das estações
plot(wsdem)
text(estacoes,col = "black")
plot(estacoes, col = "black", legend = FALSE)
plot(wsdem, add = TRUE)
text(estacoes, col = "black")
# Codigo da Bacia (ANA)
if(basin_name %in% c("mogu", "pauliceia")) basin_id <- 6
# Diretótio da Bacia
basin_path <- gsub(pattern = "X",
replacement = basin_id,
x = "~/DBHM/data_process/geo_data/data_base/ANA/BaciaX")
# Listando os arquivos shape
shapefiles_list <- list.files(basin_path,
pattern = "shp$",
recursive = T,
full.names = T) # Arquivo das Grandes Bacias Hidrogŕaficas
basin_file <- grep(x = shapefiles_list,
pattern = paste("Bacia", basin_id),
value = T)
# Shapfile de estações fluviométricas
flu_file <- grep(x = shapefiles_list,
pattern = "Est_Fluviometricas_ANA_2010",
value = T)
# Shapefile de Rede de Drenagem
dren_file <- grep(x = shapefiles_list,
pattern = "Hidrografia 1000000",
value = T)
# Shapefiles de Ottobacias (níveis 1,2 e 4)
otto1_file <- "~/DBHM/data_process/geo_data/data_base/ANA/Ottobacias_Nivel1/Ottobacias_Nivel1.shp"
otto2_file <- "~/DBHM/data_process/geo_data/data_base/ANA/Ottobacias_Nivel2/Ottobacias_Nivel2.shp"
otto4_file <- "~/DBHM/data_process/geo_data/data_base/ANA/Ottobacias_Nivel4/Ottobacias_Nivel4.shp"
# nome do arquivo com delimitação dos estados brasileiros
br_states_file <- "~/DBHM/data_process/geo_data/data_base/br_states.rds"
# importa delimitação dos estados brasileiros
br_states <- readRDS(file = br_states_file)
# estacoes flu (vazao)
flu <- shapefile(flu_file, encoding="latin1")
# projetando flu para lonlat para salvar como kml
flu <- spTransform(x = flu, CRSobj = CRS(proj4string(br_states)))
# bacia brasileira
basin <- shapefile(basin_file, encoding="latin1")
# projetando basin para lonlat para salvar como kml
basin <- spTransform(x = basin, CRSobj = CRS(proj4string(br_states)))
# rede de drenagem
dren <- shapefile(dren_file, encoding="latin1")
# selecionando somente as variáveis de interesse no slot de dados
# do objeto dren (SpatialLinesDataFrame)
dren@data <- subset(dren@data, sel = c("NORIO", "NORIOCOMP", "COBACIA", "CORIO"))
# projetando dren para lonlat para salvar como kml
dren <- spTransform(x = dren, CRSobj = CRS(proj4string(br_states)))
# gerando arquivo KML para visualização no google-earth
KML(x = crop(dren,wsdem_shp), file = paste0("./","dren.kml"), overwrite = TRUE)
# ottobacias
otto1 <- shapefile(otto1_file, encoding="latin1")
otto2 <- shapefile(otto2_file, encoding="latin1")
otto2 <- spTransform(x = otto2, CRSobj = CRS(proj4string(br_states)))
otto4 <- shapefile(otto4_file, encoding="latin1")
otto4 <- spTransform(x = otto4, CRSobj = CRS(proj4string(br_states)))
# gerando arquivo KML para visualização no google-earth
KML(x = otto4, file = paste0("./","otto4.kml"), overwrite = TRUE)
# Bacia Hidrográfica
wsdem_shp_br <- spTransform(x = wsdem_shp, CRSobj = CRS(proj4string(br_states)))
# gerando arquivo KML para visualização no google-earth
KML(x = wsdem_shp, file = paste0("./",basin_name,".kml"), overwrite = TRUE) # estados
plot(crop(br_states ,extent(estacoes)), border = "white", col = "gray", lwd = 2, axes = T)
# ottobacias nível de divisão 4
plot(crop(otto4 ,extent(estacoes)), add = T, border = "blue")
# ottobacias nível de divisão 2
plot(crop(otto2 ,extent(estacoes)), add = T, border = "red", lwd = 3)
# ottobacias nível de divisão 1
plot(crop(otto1 ,extent(estacoes)), add = T, border = "black", lwd= 4)
# estações em numero na tabelo info
text(estacoes, col = "black")
# contorno da bacia
plot(wsdem_shp, add = TRUE, border = "green")
polys <- as.data.frame(wsdem_shp@polygons[[1]]@Polygons[[1]]@coords)
# Mapa interativo com localização das estações
gmap(lat = mean(extent(basin_lonlat)[3:4]),
lng = mean(extent(basin_lonlat)[1:2]), zoom = 7,
width = 900, height = 800,
# map_type = "satellite",
map_type = "hybrid",
tools = c("pan","wheel_zoom","box_zoom","reset","resize")) %>%
ly_points(lon, lat, data = info, hover = c(site,alt,name), col = "red") %>%
ly_polygons(xs = x, ys = y, data = polys, color = "yellow", fill_alpha = 0 ) %>%
ly_abline(v = (info[info$site == "X001","lon" ])) %>%
ly_abline(h = (info[info$site == "X001","lat" ])) %>%
tool_lasso_select(select_every_mousemove = TRUE) %>%
tool_save()## xlim not specified explicitly... calculating...## ylim not specified explicitly... calculating... # Extraindo coordenadas Lat/Lon de cada celula do raster
how.many.station.in.each.cell <- data.frame( xyFromCell(cell = 1:ncell(wsdem),object = wsdem) )
cell.distance <- sapply(info$site, function(i){
b <- distGeo(p1 = how.many.station.in.each.cell[,1:2],
p2 = info[info$site == i, c("lon","lat")])
b
})
cell.distance <- as.data.frame(cell.distance)plot.raio.est <- function(raio = 100000,
cell.distance = cell.distance
){
cat("Raio seleçoado:", raio,"m\n")
cat("Construindo fatores \n")
# Juntando o codigo das estações com distância menor que o raio à célula
factor <- sapply(1:length(cell.distance[,1]),function(i){
paste0(names(cell.distance)[which(cell.distance[i,] <= raio)], collapse = "")
})
# Determinando as combinações únicas de estações para cada célula.
unique.factor <- data.frame(id = 1:length(unique(factor)),
code = unique(factor))
cat("Determinados", length(unique.factor$id),"fatores")
# Contando quantas estações vão ser utilizadas na interpolação para cada célula
for(i in 1:length(unique.factor$code)){
unique.factor$nstn[i] <- sum(str_count(unique.factor$code[i],LETTERS))
}
plot(data.frame(stations = min(unique.factor$nstn):max(unique.factor$nstn),
fatores = as.vector(table(unique.factor$nstn))),
type = "l")
factor.data <- data.frame(factor = factor,
id = rep(NA,times=length(factor)),
nstn = rep(NA,times=length(factor) ) )
for(i in as.vector(unique.factor$code[unique.factor$code != ""] )){
factor.data$id[which(factor.data$factor == i)] <- unique.factor$id[which(unique.factor$code == i)]
factor.data$nstn[which(factor.data$factor == i)] <- unique.factor$nstn[which(unique.factor$code == i)]
}
est.by.cell <- code.by.cell <- raster(wsdem)
values(est.by.cell) <- factor.data$nstn
code.by.cell <- raster(wsdem) %>%
setValues(factor.data$id) #%>%
# rasterToPolygons(dissolve = TRUE)
return(list(est.by.cell,code.by.cell,unique.factor))
} # Lendo dados observados
all.vars <- readRDS("~/Dropbox/Dissertacao/EstSP.MG_2001-2014_FINAL.rds")
# Organizando data frame por variáveis
vars <- names(select(all.vars[[1]], -date))
all.vars <- lapply(vars, function(i){
cat("\n",i, "\n")
d <- sapply(info$site, function(j){
cat(j,".")
all.vars[[j]][,i]
})
d<- cbind(all.vars[[1]][,"date"], as.data.frame(d))
names(d) <- c("date", info$site)
d
}); names(all.vars) <- vars##
## prec
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## tair
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## tair.orig
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## rh
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## rh.orig
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## ws
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## ws.orig
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## rg
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## rg.orig
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## es
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## es.orig
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## potrad
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## daytime
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## zlwd
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .
## zlwd.orig
## A701 .A705 .A706 .A707 .A708 .A711 .A712 .A713 .A714 .A715 .A716 .A718 .A725 .A726 .A727 .A728 .A729 .A733 .A734 .A735 .A736 .A737 .A738 .A739 .A740 .A741 .A745 .A746 .A747 .A748 .A753 .A755 .A502 .A505 .A506 .A507 .A508 .A509 .A510 .A511 .A512 .A513 .A514 .A515 .A516 .A517 .A518 .A519 .A520 .A521 .A522 .A523 .A524 .A525 .A526 .A527 .A528 .A529 .A530 .A531 .A532 .A533 .A534 .A535 .A536 .A537 .A538 .A539 .A540 .A541 .A542 .A543 .A544 .A545 .A546 .A547 .A548 .A549 .A550 .A551 .A552 .A553 .A554 .A555 .A556 .A557 .F501 .X001 .Raio 150 km
# Determinando áreas baixo a mesma combinação de estações
est.n.code.by.cell <- plot.raio.est(raio = 150000, cell.distance = cell.distance)## Raio seleçoado: 150000 m
## Construindo fatores
## Determinados 391 fatores
# Convertendo a poligonos
code.poly <- est.n.code.by.cell[[2]] %>% rasterToPolygons(dissolve = TRUE)
# Plot das áreas geradas
plot(est.n.code.by.cell[[1]], col = tim.colors(8), interpolate = TRUE)
plot(code.poly, add = TRUE, border = "red")
plot(wsdem_shp, add = TRUE, border = "black", lwd = 3)
text(estacoes, col = "red" )
# Levando a máscara da bacia
code.150 <- est.n.code.by.cell[[2]]
code.150[is.na(wsdem)] <- NA
# Restringindo na área da Bacia Hidrográfica
plot(code.150, col = tim.colors(length(unique(code.150))),
main = "Áreas de iguais combinações de estações\n para interpolação")
plot(code.poly, add = TRUE, border = "white")
plot(wsdem_shp, add = TRUE, border = "black", lwd = 3)
text(estacoes, col = "red" )
#Frecuência de células por combinação
freq.factor <- freq(code.150)
plot(freq.factor[!is.na(freq.factor[,1]),],
type = "h", lwd = 6)
classes <- est.n.code.by.cell[[3]]setEPS()
postscript(file = "~/Dissertação/imagens/stn.code.area.150km.eps",
onefile = TRUE,pointsize = 7,
width = 6,height = 6)
plot(code.150, col = tim.colors(length(unique(code.150))),
main = "Área de igual combinação de estações\n para interpolação")
plot(code.poly, add = TRUE, border = "white")
plot(wsdem_shp, add = TRUE, border = "black", lwd = 3)
text(estacoes, col = "red" )
for(set.id in unique(code.150)){
cat(set.id,"..")
est.names <- info$site[str_detect(string = classes$code[set.id] ,
pattern = info$site)]
plot(code.150 == set.id,
main = paste("Região",set.id,"\n",paste0(est.names,collapse = "-")),
legend = FALSE,
axes = FALSE)
}## 8 ..## 9 ..## 34 ..## 40 ..## 43 ..## 60 ..## 63 ..## 64 ..## 67 ..## 68 ..## 69 ..## 70 ..## 71 ..## 79 ..## 83 ..## 89 ..## 98 ..## 99 ..## 116 ..## 117 ..## 120 ..## 122 ..## 123 ..## 124 ..## 125 ..## 126 ..## 131 ..## 137 ..## 139 ..## 141 ..## 143 ..## 144 ..## 146 ..## 151 ..## 154 ..## 155 ..## 159 ..## 163 ..## 165 ..## 166 ..## 167 ..## 169 ..## 171 ..## 175 ..## 177 ..## 178 ..## 179 ..## 184 ..## 185 ..## 187 ..## 188 ..## 190 ..## 191 ..## 192 ..## 194 ..## 195 ..## 199 ..## 200 ..## 201 ..## 202 ..## 204 ..## 210 ..## 215 ..## 216 ..## 222 ..## 223 ..## 224 ..## 225 ..## 228 ..## 230 ..## 233 ..## 234 ..## 235 ..## 236 ..## 237 ..## 238 ..## 239 ..## 240 ..## 242 ..## 246 ..## 247 ..## 249 ..## 250 ..## 251 ..## 253 ..## 254 ..## 255 ..## 260 ..## 261 ..## 262 ..## 264 ..## 265 ..## 267 ..## 268 ..## 269 ..## 270 ..## 271 ..## 272 ..## 273 ..## 274 ..## 275 ..## 276 ..## 278 ..## 279 ..## 280 ..## 281 ..## 282 ..## 283 ..## 284 ..## 286 ..## 287 ..## 288 ..## 289 ..## 290 ..## 291 ..## 292 ..## 295 ..## 296 ..## 297 ..## 299 ..## 301 ..## 302 ..## 303 ..## 304 ..## 305 ..## 309 ..## 314 ..## 315 ..## 317 ..## 324 ..## 328 ..## 331 ..## 336 ..## 340 ..dev.off()## png
## 2 ex.regioes <- classes[order(classes$nstn),] %>% filter(id %in% unique(code.150)) %>% head
ex.regioes## id code nstn
## 1 249 A711A726A738A739A509A530X001 7
## 2 216 A711A726A738A739A747A516A530X001 8
## 3 224 A711A726A738A739A741A747A530X001 8
## 4 225 A711A726A738A739A509A516A530X001 8
## 5 240 A711A726A738A739A747A509A530X001 8
## 6 250 A706A711A726A738A739A509A530X001 8 plot(code.150 %in% ex.regioes$id, col = c("white","red"),legend = FALSE,
main = "Áreas de iguais combinações de estações\n para interpolação")
plot(code.poly, add = TRUE, border = "blue")
plot(wsdem_shp, add = TRUE, border = "black", lwd = 3)
text(estacoes, col = "black")
Séries de dados para as regiões com menor quantidade de estações
Série observada SEM preenchimento
for(set.id in ex.regioes$id){
cat(set.id,"..")
est.names <- info$site[str_detect(string = classes$code[set.id] ,
pattern = info$site)]
nstn.at.hour <- data.frame(date = all.vars[[1]][,"date"],
tair = apply(!is.na(all.vars[["tair.orig"]][,est.names]),1,sum),
prec = apply(!is.na(all.vars[["prec"]][,est.names]),1,sum),
rg = apply(!is.na(all.vars[["rg.orig"]][,est.names]),1,sum),
es = apply(!is.na(all.vars[["es.orig"]][,est.names]),1,sum),
ws = apply(!is.na(all.vars[["ws.orig"]][,est.names]),1,sum),
zlwd = apply(!is.na(all.vars[["zlwd.orig"]][,est.names]),1,sum)) %>%
mutate(tair = ifelse(tair == 0, -1,tair),
prec = ifelse(prec == 0, -1,prec),
rg = ifelse(rg == 0, -1,rg),
es = ifelse(es == 0, -1,es),
wts = ifelse(ws == 0, -1,ws),
zlwtd = ifelse(zlwd == 0, -1,zlwd)
)
timePlot(nstn.at.hour, c("tair","prec","rg","es","wts","zlwtd"),
ref.y = 0,
plot.type = "l",lwd = 3,key.columns = 3,
main = paste("Região", set.id,"\n", paste(est.names, collapse = " ")),
date.breaks = 32,date.format = "%m\n%y",
cols = tim.colors(6) )
}## 249 ..
## 216 ..
## 224 ..
## 225 ..
## 240 ..
## 250 ..
Série observada COM preenchimento
for(set.id in ex.regioes$id){
cat(set.id,"..")
est.names <- info$site[str_detect(string = classes$code[set.id] ,
pattern = info$site)]
nstn.at.hour <- data.frame(date = all.vars[[1]][,"date"],
tair = apply(!is.na(all.vars[["tair"]][,est.names]),1,sum),
prec = apply(!is.na(all.vars[["prec"]][,est.names]),1,sum),
rg = apply(!is.na(all.vars[["rg"]][,est.names]),1,sum),
es = apply(!is.na(all.vars[["es"]][,est.names]),1,sum),
ws = apply(!is.na(all.vars[["ws"]][,est.names]),1,sum),
zlwd = apply(!is.na(all.vars[["zlwd"]][,est.names]),1,sum)) %>%
mutate(tair = ifelse(tair == 0, -1,tair),
prec = ifelse(prec == 0, -1,prec),
rg = ifelse(rg == 0, -1,rg),
es = ifelse(es == 0, -1,es),
wts = ifelse(ws == 0, -1,ws),
zlwtd = ifelse(zlwd == 0, -1,zlwd)
)
timePlot(nstn.at.hour, c("tair","prec","rg","es","wts","zlwtd"),
ref.y = 0,
plot.type = "l",lwd = 3,key.columns = 3,
main = paste("Região", set.id,"\n", paste(est.names, collapse = " ")),
date.breaks = 32,date.format = "%m\n%y",
cols = tim.colors(6) )
}## 249 ..
## 216 ..
## 224 ..
## 225 ..
## 240 ..
## 250 ..
Raio 100 km
# Determinando áreas baixo a mesma combinação de estações
est.n.code.by.cell <- plot.raio.est(raio = 100000, cell.distance = cell.distance)## Raio seleçoado: 1e+05 m
## Construindo fatores
## Determinados 224 fatores
# Convertendo a poligonos
code.poly <- est.n.code.by.cell[[2]] %>% rasterToPolygons(dissolve = TRUE)
# Plot das áreas geradas
plot(est.n.code.by.cell[[1]], col = tim.colors(8), interpolate = TRUE)
plot(code.poly, add = TRUE, border = "red")
plot(wsdem_shp, add = TRUE, border = "black", lwd = 3)
text(estacoes, col = "red" )
# Levando a máscara da bacia
code.150 <- est.n.code.by.cell[[2]]
code.150[is.na(wsdem)] <- NA
# Restringindo na área da Bacia Hidrográfica
plot(code.150, col = tim.colors(length(unique(code.150))),
main = "Áreas de iguais combinações de estações\n para interpolação")
plot(code.poly, add = TRUE, border = "white")
plot(wsdem_shp, add = TRUE, border = "black", lwd = 3)
text(estacoes, col = "red" )
#Frecuência de células por combinação
freq.factor <- freq(code.150)
plot(freq.factor[!is.na(freq.factor[,1]),],
type = "h", lwd = 6)
classes <- est.n.code.by.cell[[3]]setEPS()
postscript(file = "~/Dissertação/imagens/stn.code.area.100km.eps",
onefile = TRUE,pointsize = 7,
width = 6,height = 6)
plot(code.150, col = tim.colors(length(unique(code.150))),
main = "Área de igual combinação de estações\n para interpolação")
plot(code.poly, add = TRUE, border = "white")
plot(wsdem_shp, add = TRUE, border = "black", lwd = 3)
text(estacoes, col = "red" )
for(set.id in unique(code.150)){
cat(set.id,"..")
est.names <- info$site[str_detect(string = classes$code[set.id] ,
pattern = info$site)]
plot(code.150 == set.id,
main = paste("Região",set.id,"\n",paste0(est.names,collapse = "-")),
legend = FALSE,
axes = FALSE)
}## 3 ..## 5 ..## 6 ..## 26 ..## 30 ..## 34 ..## 36 ..## 39 ..## 46 ..## 48 ..## 49 ..## 54 ..## 56 ..## 57 ..## 60 ..## 62 ..## 64 ..## 66 ..## 67 ..## 69 ..## 70 ..## 71 ..## 72 ..## 73 ..## 74 ..## 82 ..## 86 ..## 89 ..## 91 ..## 101 ..## 102 ..## 103 ..## 104 ..## 105 ..## 107 ..## 109 ..## 111 ..## 116 ..## 117 ..## 118 ..## 120 ..## 121 ..## 122 ..## 123 ..## 124 ..## 135 ..## 137 ..## 144 ..## 145 ..## 146 ..## 147 ..## 148 ..## 149 ..## 161 ..## 162 ..## 163 ..## 165 ..## 167 ..## 179 ..## 181 ..## 185 ..## 186 ..## 187 ..## 188 ..## 189 ..## 194 ..## 195 ..dev.off()## png
## 2 ex.regioes <- classes[order(classes$nstn),] %>% filter(id %in% unique(code.150)) %>% head
ex.regioes## id code nstn
## 1 3 A736A747A748 3
## 2 107 A738A739A530 3
## 3 181 A739A509A530 3
## 4 26 A736A747A748A753 4
## 5 39 A736A737A747A748 4
## 6 67 A708A711A747X001 4 plot(code.150 %in% ex.regioes$id, col = c("white","red"),legend = FALSE,
main = "Áreas de iguais combinações de estações\n para interpolação")
plot(code.poly, add = TRUE, border = "blue")
plot(wsdem_shp, add = TRUE, border = "black", lwd = 3)
text(estacoes, col = "black")
Séries de dados para as regiões com menor quantidade de estações
Série observada SEM preenchimento
for(set.id in ex.regioes$id){
cat(set.id,"..")
est.names <- info$site[str_detect(string = classes$code[set.id] ,
pattern = info$site)]
nstn.at.hour <- data.frame(date = all.vars[[1]][,"date"],
tair = apply(!is.na(all.vars[["tair.orig"]][,est.names]),1,sum),
prec = apply(!is.na(all.vars[["prec"]][,est.names]),1,sum),
rg = apply(!is.na(all.vars[["rg.orig"]][,est.names]),1,sum),
es = apply(!is.na(all.vars[["es.orig"]][,est.names]),1,sum),
ws = apply(!is.na(all.vars[["ws.orig"]][,est.names]),1,sum),
zlwd = apply(!is.na(all.vars[["zlwd.orig"]][,est.names]),1,sum)) %>%
mutate(tair = ifelse(tair == 0, -1,tair),
prec = ifelse(prec == 0, -1,prec),
rg = ifelse(rg == 0, -1,rg),
es = ifelse(es == 0, -1,es),
wts = ifelse(ws == 0, -1,ws),
zlwtd = ifelse(zlwd == 0, -1,zlwd)
)
timePlot(nstn.at.hour, c("tair","prec","rg","es","wts","zlwtd"),
ref.y = 0,
plot.type = "l",lwd = 3,key.columns = 3,
main = paste("Região", set.id,"\n", paste(est.names, collapse = " ")),
date.breaks = 32,date.format = "%m\n%y",
cols = tim.colors(6) )
}## 3 ..
## 107 ..
## 181 ..
## 26 ..
## 39 ..
## 67 ..
Série observada COM preenchimento
for(set.id in ex.regioes$id){
cat(set.id,"..")
est.names <- info$site[str_detect(string = classes$code[set.id] ,
pattern = info$site)]
nstn.at.hour <- data.frame(date = all.vars[[1]][,"date"],
tair = apply(!is.na(all.vars[["tair"]][,est.names]),1,sum),
prec = apply(!is.na(all.vars[["prec"]][,est.names]),1,sum),
rg = apply(!is.na(all.vars[["rg"]][,est.names]),1,sum),
es = apply(!is.na(all.vars[["es"]][,est.names]),1,sum),
ws = apply(!is.na(all.vars[["ws"]][,est.names]),1,sum),
zlwd = apply(!is.na(all.vars[["zlwd"]][,est.names]),1,sum)) %>%
mutate(tair = ifelse(tair == 0, -1,tair),
prec = ifelse(prec == 0, -1,prec),
rg = ifelse(rg == 0, -1,rg),
es = ifelse(es == 0, -1,es),
wts = ifelse(ws == 0, -1,ws),
zlwtd = ifelse(zlwd == 0, -1,zlwd)
)
timePlot(nstn.at.hour, c("tair","prec","rg","es","wts","zlwtd"),
ref.y = 0,
plot.type = "l",lwd = 3,key.columns = 3,
main = paste("Região", set.id,"\n", paste(est.names, collapse = " ")),
date.breaks = 32,date.format = "%m\n%y",
cols = tim.colors(6) )
}## 3 ..
## 107 ..
## 181 ..
## 26 ..
## 39 ..
## 67 ..