Preparação dos dados de forçantes atmosféricas, bacia Pauliceia

1 Introdução

O Modelo Hidrológico Distribuído da Biosfera (DBHM) é um modelo espacialmente distribuído que integra processos hidrológicos e processos de transferência solo-vegetação-atmosfera (SVAT) na escala bacia hidrográfica (Tang et al 2006).

Para simular o ciclo hidrológico terrestre o DBHM requer um conjunto de dados de entrada para domínio espacial da bacia que inclui: dados de estações meteorológicas de superfície, de sensoriamento remoto e informações pedológicas.

Nesse documento descreve-se o procedimento para interpolação dos dados meteorológicos diários para a grade da região de estudo da bacia hidrográfica. A bacia hidrográfica usada como exemplo é a bacia do Rio Mogi-Guaçú (MoGu) que abrange parte dos estados de MG e SP.

O passo de tempo de cálculo do DBHM é horário, mas dados nessa escala de tempo são disponíveis apenas a partir de 2008 no Brasil, quando as estações meteorológicas automáticas (EMA) do INMET foram implementadas efetivamente. Essas séries temporais são relativamente curtas e com baixa densidade espacial. Conjuntos de dados diários são mais facilmente encontrados, a rede de estações convencionais tem maior densidade e cobertura espacial; e as séries temporais são mais extensas.

Devido a esse fato, comum à outros países, o DBHM utiliza dados diários como entrada. A desagregação temporal dos dados diários para escala horária é realizada por métodos de downscaling temporal disponíveis em sub-rotinas do DBHM durante uma simulação.

2 Estrutura de diretórios e organização dos dados

Para reproduçao desse tutorial assume-se que o Sistema Operacional é Linux e que os dados estejam organizados segundo a estrutura de diretórios mostrada na Figura 1.

Figura 1. Estrutura de diretórios e organização dos dados.

O processamento dos dados atmosféricos de entrada para o DBHM serão ralizados no diretório /home/user/DBHM/data_process/atm_data.

2.1 Diretório de trabalho, scripts e pacotes

O diretório /home/user/DBHM/data_process/atm_data/doc/ será usado como diretório de trabalho para as operações realizadas no RStudio. Este documento (dbhm_atm.Rmd) foi escrito em R Markdown usando o RStudio e encontra-se nesse mesmo diretório.

rm(list = ls()) 
# nome do usuário para construção do caminho de diretórios
user_name <- Sys.info()[["user"]]
(wd <- paste0("/home/"
              ,user_name
              ,"/DBHM/data_process/atm_data/doc"))

[1] "/home/hidro2roilan/DBHM/data_process/atm_data/doc"

              #,"/DBHM_interp/data_process/atm_data/doc"))
# verifica estrutura de diretórios e define dir de trabalho
if(dir.exists(wd)) setwd(wd) else stop("Estrutura de diretórios não existente ou incorreta!")
getwd()

[1] "/home/hidro2roilan/DBHM/data_process/atm_data/doc"

Carregando pacotes necessários e definindo opções de configuração de alguns pacotes:

# definindo globalmente tz = "GMT"
Sys.setenv(TZ = 'GMT')
# data manipulation packages
# install.packages(
#    'printr',
#    type = 'source',
#    repos = c('http://yihui.name/xran', 'http://cran.rstudio.com')
# )
pcks <- c("knitr", "printr", "knitcitations",
          "R.utils", "magrittr", "lubridate","stringr", 
          "plyr", "dplyr","raster", "lattice","printr",
          "rasterVis")
invisible(sapply(pcks, require, character.only = TRUE, quietly = TRUE))

# configuraçoes knitcitations
cleanbib()
cite_options(citation_format = "text",
            cite.style = "authoryear", 
            style = "html", 
            hyperlink = "to.bib")
# configuraçoes knitr
opts_chunk$set(cache = FALSE)

Chunk de configuração:

    # Nome da bacia
    bacia <- "pauliceia"
    # Vetor com o intervalo de tempo a interpolar
    year.init <- c(2009,7,1)
    year.fin  <- c(2014,8,31)
    # Resolução da simulação
    resolution <- 0.5

Carregando scripts:

# script para criaçao de lista de diretorios de uma simulaçao
source("../R/set_sim_paths.R")
source("../R/start_end_sim.R")
source("../R/est.to.dbhm.R")
source("../R/getBdmepInmetData_v3.R")

2.2 Estrutura de diretórios do DBHM

Para interpolação dos dados diários das estações hidrometeorológicas é necessário entender a estrutura de organização dos diretórios do DBHM (Figura 2). Os 2 subdiretórios principais são: MODEL e ATM_INPUT_DATA. A estrutura de diretórios abaixo do diretório DBHM deve estar contida no diretório /home/user/DBHM/simulations/mogu/1km/2004/ da Figura 1.

Figura 2. Estrutura de diretórios do DBHM.

O diretório DBHM/simulations/pauliceia/0.5/tutorial_interp da Figura 1 é equivalente ao diretório DBHM/MODEL da Figura 2.

O diretório DBHM/simulations/pauliceia/0.5/2009-2014/tutorial_interp é definido pela variável sim_d (diretório de simulação):

getwd()

[1] "/home/hidro2roilan/DBHM/data_process/atm_data/doc"

sim_d <- "../../../simulations/bacia/resolutionkm/year/tutorial_interp_intel10"
sim_d <- gsub(pattern = "bacia", x = sim_d, replacement = bacia)
sim_d <- gsub(pattern = "year", x =  sim_d, replacement = paste(year.init[1],year.fin[1],sep = "-") )
sim_d <- gsub(pattern = "resolution", x =  sim_d, replacement = resolution)
sim_d

[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10"

Para proceder a interpolação o diretório sim_d deve existir, se não existir é criado.

if(!dir.exists(sim_d)) mkdirs(sim_d)
stopifnot(dir.exists(sim_d))

2.3 Tipos de dados entrada

O DBHM possui dois tipos de dados entrada:

• Variantes no tempo: armazenados no diretório DBHM/ATM_INPUT_DATA e incluem os arquivos binários de:

  • forçantes meteorológicas interpoladas para o domínio da bacia armazenados em DBHM/ATM_INPUT_DATA/Interpolate2;

  • evapotranspiração de referência (ET0) armazenados em DBHM/ATM_INPUT_DATA/ETPATH2

  • Índice de área foliar (LAI) e fração de radiação fotossintéticamente ativa (PAR) absorvida obtidos do MODIS DBHM/ATM_INPUT_DATA/AVHRR/[LAI, FPAR]

• Estáticos (invariantes no tempo): armazenados no diretório DBHM/MODEL/GEO_INPUT_DATA
  • arquivos raster que caracterizam a geomofologia da bacia
  • tabelas com informações geográficas das estações hidrometeorológicas

3 Fonte de dados

Os dados meteorológicos podem ser obtidos de diversas fontes:

• INMET (Brasil)
• ANA (Brasil)
• DAEE (SP)
• CIIAGRO (SP)
• FEPAGRO (RS)
• SIMEPAR (PR)
• EPAGRI/CIRAM (SC)

Para bacia MoGu foram obtidos dados do INMET da ANA e do DAEE.

4 Arquivos de entrada

4.1 Dados das estações hidrometeorológicas

As variáveis meteorológicas necessárias para o DBHM são:

• Precipitação
• Temperatura do ar mínima
• Temperatura do ar máxima
• Umidade relativa média do ar
• Velocidade media do vento
• Número de horas de brilho solar
• Fração de cobertura de nuvens (opcional)

Os dados de cada estação devem ser armazenados em um arquivo texto ASCII com 21 colunas e sem cabeçalho. Exemplos de arquivos estão no diretório ../data_base/mogu/inmet/atm_input_data/[0-9].*txt$. As colunas correspondem as seguintes variáveis:

• stn_id: código identificador da estação meteorológica
• year: ano
• month: mês
• day: dia
• tm: temperatura média diária do ar (0.1 °C)
• tmax: temperatura máxima do ar (0.1 °C)
• tmin: tempratura mínima do ar (0.1 °C)
• um: umidade relativa média diária do ar (%)
• umin: umidade relativa mínima do ar (%)
• n_sum: fração de cobertura de nuvens (0.01)
• n_lowm: fração de cobertura de nuvens baixas (0.01)
• fsm: velocidade média diária do vento (0.1 m s^-1)
• fmaxx: direção da velocidade máxima do vento (°)
• fmaxs: velocidade máxima do vento (0.1 m s^-1)
• rsum: precipitação total diária (0.1 mm dia^-1)
• d0m: temperatura de relva (0.1 °C)
• sun: número de horas de brilho solar (0.1 h)
• E01: evaporação diária (mm dia^-1)
• snow: neve (cm)
• mslp: pressão média ao nível do mar (mb)
• patm: pressão atmosférica (mb)

As variáveis em negrito são essenciais para o DBHM. Todas variáveis estão salvas no arquivo como números inteiros. Os fatores de conversão para valores reais são indicados entre parênteses.

O trecho de código R abaixo mostra uma amostra dos dados de um arquivo. Para mostrar o número de colunas foi criado um cabeçalho definido pelo prefixo col e o sufixo nº da coluna, mas este cabeçalho não consta no arquivo de dados original. O código 32766 é usado como rótulo para identificação de dado faltante.

data.dir <- "../data_base/bacia/dbhm_atm_input_data"
    data.dir <- gsub(pattern = "bacia", replacement = bacia, x = data.dir)
    data.dir

[1] "../data_base/pauliceia/dbhm_atm_input_data"

    data.files <- list.files(path = data.dir, pattern = ".txt$")
    
stn_data <- read.table(paste(data.dir,data.files[1],sep = "/"))

names(stn_data) <- paste0("col", 1:ncol(stn_data))
head(stn_data)

col1	col2	col3	col4	col5	col6	col7	col8	col9	col10	col11	col12	col13	col14	col15	col16	col17	col18	col19	col20	col21
83669	2000	1	1	215	237	151	94	32766	32766	32766	33	36	32766	279	32766	97	32766	32766	32766	896
83669	2000	1	2	207	249	184	93	32766	32766	32766	37	36	32766	140	32766	100	32766	32766	32766	896
83669	2000	1	3	183	200	171	96	32766	32766	32766	33	360	32766	721	32766	100	32766	32766	32766	895
83669	2000	1	4	197	208	178	96	32766	32766	32766	23	32	32766	396	32766	100	32766	32766	32766	898
83669	2000	1	5	207	236	182	89	32766	32766	32766	27	360	32766	769	32766	97	32766	32766	32766	901
83669	2000	1	6	221	271	181	88	32766	32766	32766	13	36	32766	26	32766	63	32766	32766	32766	902

# estrutura e formato dos dados
glimpse(tbl_df(stn_data))

Observations: 5479
Variables:
$ col1  (int) 83669, 83669, 83669, 83669, 83669, 83669, 83669, 83669, ...
$ col2  (int) 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 20...
$ col3  (int) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
$ col4  (int) 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1...
$ col5  (int) 215, 207, 183, 197, 207, 221, 232, 236, 248, 241, 252, 2...
$ col6  (int) 237, 249, 200, 208, 236, 271, 275, 284, 287, 288, 296, 2...
$ col7  (int) 151, 184, 171, 178, 182, 181, 178, 176, 192, 181, 193, 1...
$ col8  (int) 94, 93, 96, 96, 89, 88, 86, 62, 64, 68, 60, 73, 66, 61, ...
$ col9  (int) 32766, 32766, 32766, 32766, 32766, 32766, 32766, 32766, ...
$ col10 (int) 32766, 32766, 32766, 32766, 32766, 32766, 32766, 32766, ...
$ col11 (int) 32766, 32766, 32766, 32766, 32766, 32766, 32766, 32766, ...
$ col12 (int) 33, 37, 33, 23, 27, 13, 27, 13, 17, 13, 17, 23, 13, 10, ...
$ col13 (int) 36, 36, 360, 32, 360, 36, 27, 9, 360, 360, 36, 360, 36, ...
$ col14 (int) 32766, 32766, 32766, 32766, 32766, 32766, 32766, 32766, ...
$ col15 (int) 279, 140, 721, 396, 769, 26, 161, 18, 0, 0, 0, 0, 242, 0...
$ col16 (int) 32766, 32766, 32766, 32766, 32766, 32766, 32766, 32766, ...
$ col17 (int) 97, 100, 100, 100, 97, 63, 80, 73, 47, 23, 23, 40, 70, 4...
$ col18 (int) 32766, 32766, 32766, 32766, 32766, 32766, 32766, 32766, ...
$ col19 (int) 32766, 32766, 32766, 32766, 32766, 32766, 32766, 32766, ...
$ col20 (int) 32766, 32766, 32766, 32766, 32766, 32766, 32766, 32766, ...
$ col21 (int) 896, 896, 895, 898, 901, 902, 902, 902, 903, 904, 903, 9...

Abaixo mostra-se o nome da variável correspondente a cada coluna do arquivo.

var_names <- c("stn_id", "year","month","day","tm","tmax","tmin","um","umin",
               "n_sum","n_lowm","fsm","fmaxx","fmaxs","rsum","d0m","sun","E01",
               "snow","mslp","patm")
names(stn_data) <- var_names
head(stn_data)

stn_id	year	month	day	tm	tmax	tmin	um	umin	n_sum	n_lowm	fsm	fmaxx	fmaxs	rsum	d0m	sun	E01	snow	mslp	patm
83669	2000	1	1	215	237	151	94	32766	32766	32766	33	36	32766	279	32766	97	32766	32766	32766	896
83669	2000	1	2	207	249	184	93	32766	32766	32766	37	36	32766	140	32766	100	32766	32766	32766	896
83669	2000	1	3	183	200	171	96	32766	32766	32766	33	360	32766	721	32766	100	32766	32766	32766	895
83669	2000	1	4	197	208	178	96	32766	32766	32766	23	32	32766	396	32766	100	32766	32766	32766	898
83669	2000	1	5	207	236	182	89	32766	32766	32766	27	360	32766	769	32766	97	32766	32766	32766	901
83669	2000	1	6	221	271	181	88	32766	32766	32766	13	36	32766	26	32766	63	32766	32766	32766	902

Abaixo mostra-se a conversão das variáveis unidades SI.

# dados reais
stn_data_r <- stn_data
# substituição de 32766 por NA
stn_data_r[stn_data_r == 32766] <- NA
# conversão para valores
stn_data_r <- mutate(stn_data_r
                     ,tmax = tmax/10
                     ,tm = tm/10
                     ,tmin = tmin/10
                     ,um = um
                     ,fsm = fsm/10
                     ,sun = sun/10
                     ,rsum = rsum/10)
head(stn_data_r)

stn_id	year	month	day	tm	tmax	tmin	um	umin	n_sum	n_lowm	fsm	fmaxx	fmaxs	rsum	d0m	sun	E01	snow	mslp	patm
83669	2000	1	1	21.5	23.7	15.1	94	NA	NA	NA	3.3	36	NA	27.9	NA	9.7	NA	NA	NA	896
83669	2000	1	2	20.7	24.9	18.4	93	NA	NA	NA	3.7	36	NA	14.0	NA	10.0	NA	NA	NA	896
83669	2000	1	3	18.3	20.0	17.1	96	NA	NA	NA	3.3	360	NA	72.1	NA	10.0	NA	NA	NA	895
83669	2000	1	4	19.7	20.8	17.8	96	NA	NA	NA	2.3	32	NA	39.6	NA	10.0	NA	NA	NA	898
83669	2000	1	5	20.7	23.6	18.2	89	NA	NA	NA	2.7	360	NA	76.9	NA	9.7	NA	NA	NA	901
83669	2000	1	6	22.1	27.1	18.1	88	NA	NA	NA	1.3	36	NA	2.6	NA	6.3	NA	NA	NA	902

tail(stn_data_r)

	stn_id	year	month	day	tm	tmax	tmin	um	umin	n_sum	n_lowm	fsm	fmaxx	fmaxs	rsum	d0m	sun	E01	snow	mslp	patm
5474	83669	2014	12	26	27.3	32.3	21.2	54	NA	NA	NA	1.3	5	NA	0.0	NA	3.7	NA	NA	NA	NA
5475	83669	2014	12	27	27.7	32.5	21.6	52	NA	NA	NA	1.3	360	NA	0.0	NA	1.7	NA	NA	NA	NA
5476	83669	2014	12	28	27.6	31.5	22.4	58	NA	NA	NA	1.0	360	NA	0.0	NA	2.3	NA	NA	NA	NA
5477	83669	2014	12	29	24.7	29.9	20.2	78	NA	NA	NA	3.0	360	NA	2.1	NA	5.3	NA	NA	NA	NA
5478	83669	2014	12	30	24.9	30.4	18.1	69	NA	NA	NA	1.7	14	NA	0.0	NA	5.0	NA	NA	NA	NA
5479	83669	2014	12	31	24.1	28.6	19.4	76	NA	NA	NA	2.0	360	NA	0.2	NA	6.7	NA	NA	NA	NA

Note que a precisão dos dados é de 3 dígitos (apenas 1 casa decimal) para a maioria das variáveis, exceto a umidade relativa em que a precisão é de 2 dígitos (nenhuma casa decimal).

4.2 Informações das estações hidrometeorológicas

Além dos dados meteorológicos é necessário um arquivo auxiliar com as informações geográficas de cada estação hidrometeorológica:

• código identificador da estação hidrometeorológica
• longitude
• latitude
• altitude

O código da estação faz a relação com o arquivo de dados da estação que também esta variável (stn_id).

O arquivo de informações das estações é chamado allstn.txt.

O arquivo de exemplo para a bacia do Rio MoGu com as informações das estações hidrometeorológicas disponíveis para o ano de 2004 está disponível no diretório ../data_base/mogu/dbhm_atm_input_data_2004, aqui referenciado por meteo_data_d.

Esse arquivo não possui cabeçalho, conforme mostrado abaixo.

meteo_data_d <- "../data_base/bacia/dbhm_atm_input_data"
    meteo_data_d <- gsub(x = meteo_data_d, pattern = "bacia", replacement = bacia)
    meteo_data_d

[1] "../data_base/pauliceia/dbhm_atm_input_data"

(allstn_file <- paste0(meteo_data_d,"/allstn.txt"))

[1] "../data_base/pauliceia/dbhm_atm_input_data/allstn.txt"

head(readLines(allstn_file))

[1] "       83669   -47.5500      -21.48000        617.390"
[2] "       83726   -47.8700      -21.97000        856.000"
[3] "       83630   -47.3700      -20.58000       1026.200"
[4] "       83851   -47.4333      -23.48333        645.000"
[5] "       99010   -47.6332      -21.61945        705.000"
[6] "       99009   -47.6398      -21.64970        591.000"

A leitura do arquivo allstn_file com a função read.fwf() permite verificar a largura reservada para cada coluna do arquivo através dos valores informados para o argumento widths. Dados faltantes de altitude da estação são rotulados pelo valor -999.000. Note que não pode haver valores faltantes para o código e as coordenadas das estações nesse arquivo, pois o modelo precisa dessas informações para calcular as distâncias entre as estações.

# leitura de arquivo com formato de largura fixa
info_stn <- read.fwf(allstn_file 
                     ,widths = c(12, 11, 15, 15)
                     ,na.strings = "       -999.000")
# criando cabeçalho
names(info_stn) <- c("stn_id", "lon", "lat", "alt")
head(info_stn)

stn_id	lon	lat	alt
83669	-47.5500	-21.48000	617.39
83726	-47.8700	-21.97000	856.00
83630	-47.3700	-20.58000	1026.20
83851	-47.4333	-23.48333	645.00
99010	-47.6332	-21.61945	705.00
99009	-47.6398	-21.64970	591.00

# estrutura e formato dos dados
 glimpse(tbl_df(info_stn))

Observations: 7
Variables:
$ stn_id (int) 83669, 83726, 83630, 83851, 99010, 99009, 83683
$ lon    (dbl) -47.5500, -47.8700, -47.3700, -47.4333, -47.6332, -47.6...
$ lat    (dbl) -21.48000, -21.97000, -20.58000, -23.48333, -21.61945, ...
$ alt    (dbl) 617.39, 856.00, 1026.20, 645.00, 705.00, 591.00, 873.35

4.3 Arquivos para a bacia de Pé de Gigante.

Para a criação dos dados da bacia de Pé de Gigante temos os dados das estações do municipio no qual está localizada, Luis Antonio. Estão disponibilizados os da torre meteorologica no interior do dominio da encerra a área de estudo.

As estações convencionais a ser utilizadas são as código:

• 83669 São Simão
• 83726 São Carlos
• 83630 França
• 83676 Catanduva
• 83851 Sorocaba
• 83683 Machado

 est.list <- c(83669,83726, 83630,83851,83683)
 # 83676 da problemas na hora de ler o arquivo baixado
 atual.wd <- getwd() 
 setwd("~/Dissertação/DADOS/")
    estacoes <- lapply(est.list, function(i){
                                    file1 <- paste(i,".dbhm",sep = "")
                                if(!file.exists(file1)){ 
                                    file <- paste(i,".txt",sep = "")
                                if(file.exists(file)){
                                    est.to.dbhm(readInmet(file, diario = T))
                                }else{
                                    est.to.dbhm(readInmet(getBdmepInmetData(i, 
                                                                            sdate = "01/01/2000", 
                                                                            edate = "31/12/2014"), 
                                                          diario = T))
                                }
                                }
        })
    
   files <-  list.files(pattern = ".dbhm$" )
   
   lapply(files, function(x) file.copy(paste (".", x , sep = "/"),  
          paste (atual.wd,meteo_data_d,gsub("dbhm","txt",x), sep = "/"), recursive = FALSE,  copy.mode = TRUE, overwrite = T))

## [[1]]
## [1] TRUE
## 
## [[2]]
## [1] TRUE
## 
## [[3]]
## [1] TRUE
## 
## [[4]]
## [1] TRUE
## 
## [[5]]
## [1] TRUE

 setwd(atual.wd)
 lapply(est.list, FUN = function(i){info_stn[info_stn$stn_id == i,]} )

## [[1]]
##   stn_id    lon    lat    alt
## 1  83669 -47.55 -21.48 617.39
## 
## [[2]]
##   stn_id    lon    lat alt
## 2  83726 -47.87 -21.97 856
## 
## [[3]]
##   stn_id    lon    lat    alt
## 3  83630 -47.37 -20.58 1026.2
## 
## [[4]]
##   stn_id      lon       lat alt
## 4  83851 -47.4333 -23.48333 645
## 
## [[5]]
##   stn_id    lon    lat    alt
## 7  83683 -45.92 -21.67 873.35

Este arquivo contém informações de todas estações disponíveis para bacia MoGu, ou seja tanto as do INMET quanto da ANA. Uma forma simples de verificar se todas estações estão no arquivo allstn.txt é mostrada abaixo.

# Código das estações da tabela de informações
stn_id_info <- sort(info_stn$stn_id)
# Código das estações obtidos dos nomes dos arquivos de estação
stn_files <- list.files(meteo_data_d
                        ,pattern = "[0-9]{5}.*txt$"
                        ,full.names = TRUE)
head(stn_files)

[1] "../data_base/pauliceia/dbhm_atm_input_data/83630.txt"
[2] "../data_base/pauliceia/dbhm_atm_input_data/83669.txt"
[3] "../data_base/pauliceia/dbhm_atm_input_data/83683.txt"
[4] "../data_base/pauliceia/dbhm_atm_input_data/83726.txt"
[5] "../data_base/pauliceia/dbhm_atm_input_data/83851.txt"
[6] "../data_base/pauliceia/dbhm_atm_input_data/99009.txt"

stn_id_files <- sort(as.integer(gsub("\\.txt", "", basename(stn_files))))
stn_id_files

[1] 83630 83669 83683 83726 83851 99009 99010

# verificando se todos elementos dos dois vetores são iguais
#all(stn_id_files %in% stn_id_info)
all.equal(stn_id_files, stn_id_info)

[1] TRUE

Para a bacia Mogu serão utilizadas para interpolação 7 estações hidrometeorológicas. Entretanto esse número de estações é válido apenas para interpolação da precipitação, uma vez que as estações da ANA utilizadas só possuem dados de precipitação diária. A interpolação das demais variáveis foi realizada somente com os dados das estações do INMET (7 estações).

4.4 Arquivos raster

Os arquivos raster em formato ARC/INFO ASCII GRID da geomorfologia necessários para interpolação são:

• demhdr.asc: arquivo texto ASCII com o cabeçalho;
• wsdemirr.asc: modelo digital de elevação do terreno da bacia, incluindo as àreas de irrigação fora da bacia (se não há areas de irrigação, como no caso da bacia MoGu, esse arquivo é equivalente ao arquivo wsdem.asc);
• fracirr.asc: fração de cobertura da área da bacia incluindo as àreas de irrigação fora da bacia (se não há areas de irrigação, como no caso da bacia MoGu, esse arquivo é equivalente ao arquivo frac.asc);

Esses arquivos foram gerados numa etapa prévia de preparação dos dados geomorfológicos da bacia e encontram-se no diretório: /home/user/DBHM/data_process/geo_data (Figura 1).

O conteúdo do arquivo demhdr.asc para a bacia MoGu é mostrado abaixo:

      ncols         282
      nrows         213
      xllcorner     -175432.9434359
      yllcorner     -86300.9824231 
      cellsize      1000"           
      NODATA_value  -9999

• xllcorner: é a coordenada x do canto inferior à esquerda do domínio espacial (ou mínimo valor das coordenadas em x).

• yllcorner: é a coordenada y do canto inferior à esquerda do domínio espacial (ou mínimo valor das coordenadas na direção y).

• cellsize: é o espaçamento de grade (resolução espacial) em metros.

O arquivos ASCII GRID podem ser visualizados em qualquer editor de texto, inclusive no RStudio.

 raster_dir <- "~/DBHM/data_process/geo_data/"
 
 wsdem_file <- paste(raster_dir,"wsdem_",bacia,".asc", sep = "")
    file.exists(wsdem_file)

[1] TRUE

        # file.show(wsdem_file)
 frac_file <- paste(raster_dir,"frac_",bacia,".asc", sep = "")

A visualização desses arquivos no R pode ser feita conforme trecho de código abaixo.

    wsdem <- raster(wsdem_file)
        wsdem

class       : RasterLayer 
dimensions  : 14, 12, 168  (nrow, ncol, ncell)
resolution  : 500, 500  (x, y)
extent      : -2901.52, 3098.48, -3322.855, 3677.145  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /home/hidro2roilan/DBHM/data_process/geo_data/wsdem_pauliceia.asc 
names       : wsdem_pauliceia 
values      : -2147483648, 2147483647  (min, max)

        plot(wsdem)

    frac <- raster(frac_file)
        plot(frac)

4.5 Sistema de Projeção Cartográfica do DBHM

Pelos gráficos anteriores observa-se que as coordenadas x e y da região não são no sistema de coordenadas de latitude e longitude. A projeção cartográfica usada no DBHM é a projeção cartográfica Azimutal de Área Equivalente de Lambert (LAEA). Ela preserva a área constante em toda superfície do mapa e representa realisticamente as direções a partir de seu centro de projeção, mas as deformações na escala aumentam com o distanciamento do ponto central da projeção. A escala diminui com a distância do centro ao longo dos raios e aumenta com a distância do centro em uma direção perpendicular ao raios. A base de dados hidrográficos globais HYDRO1k usa esta projeção. Portanto para usar dados em projeções diferentes (p.ex.: a geográfica também conhecida por longlat ou latlon) no DBHM é necessário projetá-los na projeção LAEA. As Equações de transformação da LAEA para o sistema de referência geográfica estão disponíveis em http://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html.

Os parâmetros da projeção LAEA variam por continente. Os parâmetros para o continente sul americano são 60º W para longitude de origem e 15ºS para latitude de origem e o raio da esfera de influência é de 6370997 m. Para as bacias hidrográficas, podemos definir esses dois parâmetros como as coordenadas do centróide da bacia. Entretanto no caso da bacia MoGu os parâmetros utilizados foram lon_0 = -47.33 ºW e lat_0 = -21.81 ºS.

centroide <- c(lon = -47.63, lat = -21.64)
# string com argumentos do sist. coords de ref DBHM
laea_string <- "+proj=laea +lat_0=LAT +lon_0=LON +x_0=0 +y_0=0 +ellps=sphere"
laea <- gsub("LON", as.character(centroide[["lon"]]), laea_string) 
laea <- gsub("LAT", as.character(centroide[["lat"]]), laea)
laea

[1] "+proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +ellps=sphere"

projection(wsdem)

[1] NA

# definindo sist. de coords 
projection(wsdem) <- projection(frac) <- laea
wsdem

class       : RasterLayer 
dimensions  : 14, 12, 168  (nrow, ncol, ncell)
resolution  : 500, 500  (x, y)
extent      : -2901.52, 3098.48, -3322.855, 3677.145  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +ellps=sphere 
data source : /home/hidro2roilan/DBHM/data_process/geo_data/wsdem_pauliceia.asc 
names       : wsdem_pauliceia 
values      : -2147483648, 2147483647  (min, max)

frac

class       : RasterLayer 
dimensions  : 14, 12, 168  (nrow, ncol, ncell)
resolution  : 500, 500  (x, y)
extent      : -2901.52, 3098.48, -3322.855, 3677.145  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +ellps=sphere 
data source : /home/hidro2roilan/DBHM/data_process/geo_data/frac_pauliceia.asc 
names       : frac_pauliceia 
values      : -2147483648, 2147483647  (min, max)

5 Interpolação

Nessa etapa serão criados os arquivos binários com os campos espaciais interpolados das variáveis meteorológicas de entrada do DBHM na escala diária para o domínio da bacia MoGu. A interpolação será feita para o ano de 2004. O procedimento para interpolação de dados para outros anos é o mesmo.

year_sim <- c(2009,2014)

Os dados meteorológicos interpolados são usados para derivar os binários de ET0 segundo as recomendações da FAO. A ET0 é usada como estimativa da evaporação para os pontos de grade dentro da bacia classificados como corpos d’água.

5.1 Métodos de Interpolação

O DBHM disponibiliza três métodos de interpolação:

• ponderada pela distância angular - ADW (New, Hulme, and Jones (2000))

• média da área dos Polígonos de Thiessen (também conhecido por regiões de Dirichlet ou polígonos de Voronoi) (THIESSEN (1911))

• Thin-plate splines (HUTCHINSON (1995))

Nesse tutorial será usado o método ADW.

5.2 Colocando os dados necessários para Interpolação no DBHM

Os arquivos de entrada devem ser adequadamente alocados nos diretórios do DBHM. A lista de diretórios para um sim_d pode ser obtida pela função set_sim_paths:

stopifnot(existsFunction("set_sim_paths"))
sim_dirs <- set_sim_paths(loc = sim_d
                          ,abs.path = FALSE
                          ,dirs.only = TRUE
                          ,AVHRR = TRUE
                          ,create.graph.dir = TRUE)

sim_dirs

$sim_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10"

$source_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/SOURCE"

$parameter_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/PARAMETER"

$include_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/INCLUDE"

$geo_input_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/GEO_INPUT_DATA"

$atm_input_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA"

$atm_input_gauge_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/gauge_13item"

$atm_input_interp_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2"

$et_input_interp_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/ETPATH2"

$atm_input_lai_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/AVHRR/LAI"

$atm_input_fpar_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/AVHRR/FPAR"

$derive_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/Derive_D"

$output_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/HYDRO_OUTPUT"

$graph_output_d
[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/SOURCE/graphics"

lapply(sim_dirs, FUN = dir.exists)

$sim_d
[1] TRUE

$source_d
[1] TRUE

$parameter_d
[1] TRUE

$include_d
[1] TRUE

$geo_input_d
[1] TRUE

$atm_input_d
[1] TRUE

$atm_input_gauge_d
[1] TRUE

$atm_input_interp_d
[1] TRUE

$et_input_interp_d
[1] TRUE

$atm_input_lai_d
[1] TRUE

$atm_input_fpar_d
[1] TRUE

$derive_d
[1] TRUE

$output_d
[1] TRUE

$graph_output_d
[1] TRUE

5.2.1 Arquivos de dados estáticos

1. Gerando arquivo demhdr.asc e salvando-o no diretório de dados estáticos do DBHM (GEO_INPUT_DATA).

# cabeçalho do arquivo wsdem (modelo de elavação do terreno para a área da bacia)
demhdr_content <- readLines(wsdem_file, n = 6)
# local do arquivo no DBHM
(demhdr_file <- file.path(sim_dirs$geo_input_d, "demhdr.asc"))

[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/GEO_INPUT_DATA/demhdr.asc"

# escrevendo arquivo
if(dir.exists(sim_dirs$geo_input_d)){
   writeLines(text = demhdr_content 
              ,con = demhdr_file)
} else { stop("Diretório de dados estáticos (" 
              ,sim_dirs$geo_input_d
              ,")" 
              ," não encontrado:") 
}# end if
# verificação da existência do demhdr_file 
stopifnot(file.exists(demhdr_file))

2. Copiando arquivos raster wsdemirr.asc e fracirr.asc para o diretório de dados estáticos

file.copy(from = wsdem_file
          ,to = file.path(sim_dirs$geo_input_d, "wsdemirr.asc")
          ,overwrite = TRUE)

[1] TRUE

file.copy(from = frac_file
          ,to = file.path(sim_dirs$geo_input_d, "fracirr.asc")
          ,overwrite = TRUE)

[1] TRUE

3. Copiando arquivos de informações das estações hidrometeorológicas para o diretório de dados estáticos;

file.copy(from = allstn_file
          ,to = file.path(sim_dirs$geo_input_d
                          ,basename(allstn_file))
          ,overwrite = TRUE)

[1] TRUE

5.2.2 Arquivos de dados variantes no tempo

4. Copiando arquivos para o diretório de dados variantes no tempo estações de medidas meteorológicas (DBHM/ATM_INPUT_DATA)

 stopifnot(dir.exists(sim_dirs$atm_input_gauge_d))
 cp_stn_files <- file.copy(from = stn_files
                           ,to =  file.path(sim_dirs$atm_input_gauge_d
                                            ,basename(stn_files))
                           ,overwrite = TRUE)
 # verificação de que os arquivos foram copiados
 all(cp_stn_files)

[1] TRUE

5.3 Definição de parâmetros do DBHM para interpolação

Após a preparação dos arquivos de entrada necessários deve-se configurar os parâmetros necessários para interpolação com o DBHM.

5.3.1 Arquivo Input.par

No arquivo ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/PARAMETER/Input.par devem ser definidos os seguintes parâmetros, conforme a bacia e o período de dados disponíveis:

• startyear: ano inicial da simulação (2009, 2014);
• startmonth: mês inicial da simulação (1);
• startday: dia inicial da simulação (1);
• endyear: ano dinal da simulação (2009, 2014);
• endmonth: mês final da simulação (12);
• endday: dia final da simulação (31);

• DEM_FD: caminho para o arquivo wsdemirr.asc,

  ../GEO_INPUT_DATA/wsdemirr.asc

• FRAC_FD: caminho para o arquivo wsdemirr.asc,

  ../GEO_INPUT_DATA/fracirr.asc

• ETPATH: caminho para o diretório dos arquivos binários da ET0 ../../, ../../simulations/pauliceia/0, ../../ATM_INPUT_DATA/ETPATH2

• INTPLT: método de interpolação do inverso da distância (1: ADW, 2: Polígonos de Thiessen, 3: Thin-plate splines);

• ATM_ITP: caminho para o diretório dos arquivos binários dos dados meteorológicos interpolados

  ../../, ../../simulations/pauliceia/0, ../../ATM_INPUT_DATA/Interpolate2

• GEO_hrd: caminho para o arquivo demhdr.asc;

  ../GEO_INPUT_DATA/demhdr.asc

• lambda0: longitude central da projeção LAEA (em graus decimais)

  -47.63

• phi1: latitude do paralelo padrão do sistema de coordenadas (em graus decimais)

  -21.64

• stn_info_file: caminho para o arquivo allstn.txt;

  ../GEO_INPUT_DATA/allstn.txt

• stn_data_dir: caminho para o arquivos ASCII de dados das estações hidrometeorológicas;

  ../../, ../../simulations/pauliceia/0, ../../ATM_INPUT_DATA/gauge_13item

• NP: número de estações vizinhas (5);

• R0: raio de influência, 40 x resolução espacial (40 km);

• CM: coeficiente de ajuste, valor recomendado 4, ver New, Hulme, and Jones (2000);

• NRPP:, parâmetro usado no método Thin Plate Spline (valor recomendado 10);

Note que os caminhos dos arquivos são relativos ao diretório ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/SOURCE. A maioria dos parâmetros já estão pré-definidos no arquivos. Entretanto alguns precisam ser alterados.

# alterando ano inicial e final do Input.par
cat(" Alterando Input.par para o ano:", year_sim, "\n")

 Alterando Input.par para o ano: 2009 2014 

input_par_file <- file.path(sim_dirs$parameter_d, "Input.par")
ip <- readLines(input_par_file)
# atualizando startyear
input_par <- gsub("startyear = [0-9]{4}"
                 ,paste0("startyear = ", year.init[1])
                 ,ip)
input_par <- gsub("startmont = [0-9]{1}"
                 ,paste0("startmont =  ", year.init[2])
                 ,input_par)
input_par <- gsub("startday  = [0-9]{1}"
                 ,paste0("startday  =  ", year.init[3])
                 ,input_par)
# atualizando endyear
input_par <- gsub("endyear = [0-9]{4}"
                  ,paste0("endyear = ", year.fin[1])
                  ,input_par)
input_par <- gsub("endmont = [0-9]{1}"
                  ,paste0("endmont = ", year.fin[2])
                  ,input_par)
input_par <- gsub("endday  = [0-9]{2}"
                  ,paste0("endday  = ", year.fin[3])
                  ,input_par)

# atualizando coordenadas de referência da projeção
input_par <- gsub("lambda0=-?[0-9]{1,}\\.[0-9]{1,}"
                  ,paste0("lambda0=", centroide[[1]])
                  ,input_par)
input_par <- gsub("phi1=-?[0-9]{1,}\\.[0-9]{1,}"
                  ,paste0("phi1=", centroide[[2]])
                  ,input_par)
writeLines(input_par, input_par_file)

5.3.2 Arquivo commom.inc

No arquivo ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/INCLUDE/common.inc devem ser definidos os seguintes parâmetros:

• max_nc: número máximo de colunas do domínio da bacia
• max_nr: número máximo de linhas do domínio da bacia
• max_stn: número máximo de estação meteorológicas disponíveis no ano a ser interpolado (7);
• max_time: número máximo de dias simulados (365, 365);
• max_year: número máximo de anos simulados (2014);

(max_nc <- ncol(wsdem))
(max_nr <- nrow(wsdem))
(max_stn <- nrow(info_stn))
(max_time <- sum(sapply(year_sim[1]:year_sim[2],FUN = function(i){ifelse(leap_year(i), 366, 365)}) )   )
(max_year <- max(year_sim[2] - year_sim[1] + 1))
cat(" Alterando common.inc:", "\n")
common_file <- file.path(sim_dirs$include_d, "common.inc")
common <- readLines(common_file)
common <- gsub("max_nc=[0-9]{1,}"
               ,paste0("max_nc=", max_nc)
               ,common)
common <- gsub("max_nr=[0-9]{1,}"
               ,paste0("max_nr=", max_nr)
               ,common)
# grep("max_stn=[0-9]{1,3},", comm)
common <- gsub("max_stn=[0-9]{1,}"
               ,paste0("max_stn=", max_stn)
               ,common)
common <- gsub("max_time=[0-9]{3}"
               ,paste0("max_time=", max_time)
               ,common)
writeLines(common, common_file) 

5.4 Compilação e execução da interpolação

Para realizar a interpolação os diretórios ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2 e ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/ETPATH2 devem existir e não conterem arquivos binários.

stopifnot(length(list.files(sim_dirs$atm_input_interp_d, pattern = "bin$")) == 0 & 
          length(list.files(sim_dirs$et_input_interp_d, pattern = "bin$")) == 0)

Se esses arquivos já existirem ocorrerá um erro ao tentar compilar o modelo.

A interpolação é realizada com o código ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/SOURCE/Preprocess/interpolate.f que deve ser copiado para o diretório ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/SOURCE.

Nesse arquivo, os anos de início e fim do laço temporal devem ser alterados de acordo com o período de dados; neste exemplo o laço é de:

do i=2009, 2014,2009, 2014

A seguir salve, compile e rode o código interpolate.f.

cat(" Alterando Interpolation.f", year_sim, "\n")

 Alterando Interpolation.f 2009 2014 

interp_code_file <- file.path(sim_dirs$source_d, "interpolation.f")
interp_code <- readLines(interp_code_file)
# grep("i=[0-9]{4},[0-9]{4}", interp.for)
interp_code <- gsub("i=[0-9]{4},[0-9]{4}"
                    ,paste0("i=", year_sim[1], ",", year_sim[2])
                    ,interp_code)
writeLines(interp_code, interp_code_file)  

Para compilá-lo pode-se usar o comando:

ifort interpolation.f -assume byterecl -o interp2004

Para executar:

./interp2004

# caminho para compilador intel fortran
IFORT <- "/opt/intel10/fc/10.1.018/bin/ifort"
## compilando e executando intepolation.f
setwd(sim_dirs$source_d)
getwd()
cmd <- "ifort interpolation.f -assume byterecl -o interpYYYY"
(cmd <- gsub("ifort", IFORT, gsub("YYYY",year_sim[1], cmd)))
#system("source ~/.bashrc")
system(cmd, intern = FALSE, wait = TRUE)
cat("Executando Interpolation.f", year_sim[1], "...", "\n")
system(gsub("YYYY", year_sim[1],"./interpYYYY > telaYYYY.txt")
       ,intern = FALSE
       ,wait = TRUE)
setwd(wd)

Para acompanhar a execução da interpolação abra um terminal e no diretório de simulação utilize o comando linux tail -f tela2009.txt.

No caso da bacia Mogu há várias estações apenas com dados de precipitação. Um aviso será impresso na tela (que foi redirecionado para aquele arquivo texto) para essas estações. As variáveis meteorológicas dessas estações terão valor 32766.

Durante a execução do programa pode aparecer a mensagem:

NP should be less than NPI Error occur in interpolate

Isso Indica que que alguma variavel (p.ex. n_summ - cobertura de nuvens) não tem valores válidos (somente dados faltantes).

Caso ocorra o erro:

undefined reference to ana_dataf

comente todas linhas com CALL ana_DATAF(...) no código ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/SOURCE/SUBROUTINE/Get_Grid_ATM.f da subrotina Get_Grid_ATM_d.

Após a finalização da interpolação verifique se os arquivos binários em ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2 e ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/ETPATH2 foram gerados.

6 Checklist para interpolação dos dados

1. Gerar dados meteorológicos de estações;

2. Colocar os arquivos geográficos da bacia (demhdr.asc, fracirr.asc, wsdemirr.asc) no diretório ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/GEO_INPUT_DATA;

3. Colocar os arquivos de estação no diretório ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/gauge_13item;

4. Colocar o arquivo allstn.txt no diretório ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/GEO_INPUT_DATA;

5. Definir os parâmetros em ../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/PARAMETER/Input.par

6. Definir os parâmetros em r

7. Alterar o período de anos a ser interpolado no código interpolate.f, compilá-lo e rodá-lo;

7 Visualização dos dados interpolados

Com os campos espaciais da precipitação é possível calcular a precipitação média diaria sobre a área da bacia e campo espacial da precipitação total anual.

# lista de arquivos binários interpolados
bin_interp_files <- list.files(sim_dirs$atm_input_interp_d
                               ,pattern = "^[0-9]{4}.*\\.bin"
                               ,full.names = TRUE)
et_bin_files <- list.files(sim_dirs$et_input_interp_d
                           ,pattern = "^[0-9]{4}.*\\.bin"
                           ,full.names = TRUE)
bin_interp_files <- c(bin_interp_files, et_bin_files)
bin_interp_files

 [1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2009_fsm_IDW.bin"   
 [2] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2009_n_summ_IDW.bin"
 [3] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2009_rsum_IDW.bin"  
 [4] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2009_sun_IDW.bin"   
 [5] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2009_tmax_IDW.bin"  
 [6] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2009_tm_IDW.bin"    
 [7] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2009_tmin_IDW.bin"  
 [8] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2009_um_IDW.bin"    
 [9] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2010_fsm_IDW.bin"   
[10] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2010_n_summ_IDW.bin"
[11] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2010_rsum_IDW.bin"  
[12] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2010_sun_IDW.bin"   
[13] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2010_tmax_IDW.bin"  
[14] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2010_tm_IDW.bin"    
[15] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2010_tmin_IDW.bin"  
[16] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2010_um_IDW.bin"    
[17] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2011_fsm_IDW.bin"   
[18] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2011_n_summ_IDW.bin"
[19] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2011_rsum_IDW.bin"  
[20] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2011_sun_IDW.bin"   
[21] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2011_tmax_IDW.bin"  
[22] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2011_tm_IDW.bin"    
[23] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2011_tmin_IDW.bin"  
[24] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2011_um_IDW.bin"    
[25] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2012_fsm_IDW.bin"   
[26] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2012_n_summ_IDW.bin"
[27] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2012_rsum_IDW.bin"  
[28] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2012_sun_IDW.bin"   
[29] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2012_tmax_IDW.bin"  
[30] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2012_tm_IDW.bin"    
[31] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2012_tmin_IDW.bin"  
[32] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2012_um_IDW.bin"    
[33] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2013_fsm_IDW.bin"   
[34] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2013_n_summ_IDW.bin"
[35] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2013_rsum_IDW.bin"  
[36] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2013_sun_IDW.bin"   
[37] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2013_tmax_IDW.bin"  
[38] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2013_tm_IDW.bin"    
[39] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2013_tmin_IDW.bin"  
[40] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2013_um_IDW.bin"    
[41] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2014_fsm_IDW.bin"   
[42] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2014_n_summ_IDW.bin"
[43] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2014_rsum_IDW.bin"  
[44] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2014_sun_IDW.bin"   
[45] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2014_tmax_IDW.bin"  
[46] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2014_tm_IDW.bin"    
[47] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2014_tmin_IDW.bin"  
[48] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2014_um_IDW.bin"    
[49] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/ETPATH2/2009_ReferET_IDW.bin"    
[50] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/ETPATH2/2010_ReferET_IDW.bin"    
[51] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/ETPATH2/2011_ReferET_IDW.bin"    
[52] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/ETPATH2/2012_ReferET_IDW.bin"    
[53] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/ETPATH2/2013_ReferET_IDW.bin"    
[54] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/ETPATH2/2014_ReferET_IDW.bin"    

# binário de precipitação
(prec_bin_file <- bin_interp_files[grep("rsum", basename(bin_interp_files))])

[1] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2009_rsum_IDW.bin"
[2] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2010_rsum_IDW.bin"
[3] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2011_rsum_IDW.bin"
[4] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2012_rsum_IDW.bin"
[5] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2013_rsum_IDW.bin"
[6] "../../../simulations/pauliceia/0.5km/2009-2014/tutorial_interp_intel10/../ATM_INPUT_DATA/Interpolate2/2014_rsum_IDW.bin"

# valores do arquivo binário
prec_values <- readBin(prec_bin_file[2]
                       ,what = double()
                       ,n = ncell(wsdem)*365
                       ,size = 4)
# histograma da prec
histogram(prec_values)

histogram(prec_values[prec_values > 0.5])

# intervalo de variação da precipitação
range(prec_values)

[1]  0.00000 54.48429

Para sobrepor a localização das estações meteorológicas sobre o campos interpolados, vamos importar as coordenadas das estações hidrometeorológicas.

# pontos das estaçoes
stn_loc <- read.table(file.path(sim_dirs$atm_input_interp_d, paste("AVG_",year_sim[1],"_STATION.asc", sep = ""))
                      ,head = TRUE)
stn_loc %<>% dplyr::select(ID, X, Y)
coordinates(stn_loc) <- c("X", "Y")
stn_loc

class       : SpatialPointsDataFrame 
features    : 7 
extent      : -24748.55, 176696.2, -204973, 117842.5  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
variables   : 1
names       :    ID 
min values  : 83630 
max values  : 99010 

Para realizar os cálculos é melhor combinar os campos espaciais da precipitação diária em um objeto da classe RasterBrick do pacote raster.

# construindo brick com campos diários da precipitação
r <- raster(file.path(sim_dirs$geo_input_d, "wsdemirr.asc"))
prec_b <- brick(r, nl =365, values = FALSE)
# matriz com valores de todos pontos de grade para cada dia ao longo das colunas
prec_mat <- matrix(prec_values, ncol = 365)
# atribuindo valores ao brick
prec_b <- setValues(prec_b, prec_mat)
# aplicando máscara da bacia
prec_b_mask <- mask(prec_b, r)
# definindo labels para 3a dimensão dos dados (tempo)
se_dates <- start_end_sim(sim_d, seq.dates = TRUE, warn = FALSE)
se_dates <- as.Date(se_dates[se_dates == format(se_dates, paste(2010,"-%m-%d",sep = ""))])[1:365]
prec_b_mask <- setZ(prec_b_mask, z = se_dates, name = "days")
names(prec_b_mask) <- format(getZ(prec_b_mask), "dia_%Y-%m-%d")
projection(prec_b_mask) <- laea
prec_b_mask

class       : RasterBrick 
dimensions  : 14, 12, 168, 365  (nrow, ncol, ncell, nlayers)
resolution  : 500, 500  (x, y)
extent      : -2901.52, 3098.48, -3322.855, 3677.145  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +ellps=sphere 
data source : in memory
names       : dia_2010.01.01, dia_2010.01.02, dia_2010.01.03, dia_2010.01.04, dia_2010.01.05, dia_2010.01.06, dia_2010.01.07, dia_2010.01.08, dia_2010.01.09, dia_2010.01.10, dia_2010.01.11, dia_2010.01.12, dia_2010.01.13, dia_2010.01.14, dia_2010.01.15, ... 
min values  :     3.45650959,     0.25574738,     1.75369632,     0.00000000,     0.87645936,    14.02123070,     6.58049583,    17.20668602,     4.11970043,     0.30478013,     1.38694382,     8.19934082,     0.35043973,     8.30917358,     6.99771118, ... 
max values  :     4.71741772,     0.41761327,     2.86363411,     0.00000000,     1.10710585,    16.81094551,     7.11306715,    22.92314529,     5.51317787,     0.35245880,     1.90359390,    11.37146187,     0.54608876,     9.47978401,     7.22873831, ... 
days        : 2010-01-01, 2010-12-31 (min, max)

plot(prec_b_mask)

# animate(prec_b, n = 1,pause=0.1)
# prec total anual
prec_anual <- sum(prec_b_mask)
# media anual da prec na bacia
(prec_anual <- cellStats(prec_anual, mean))
# media da chuva diária na bacia
sp_avg_dly_prec <- cellStats(prec_b_mask, mean)
plot(x = se_dates
    ,y = sp_avg_dly_prec
    ,type = "h"
    ,ylab = "Prec (mm)"
    ,xlab = "date"
    ,las = 1)
# pol <- shapefile("../../geo_data/softs/TauDEM-master/mogu/mascara_bacia_alta_pol.shp")
plot(prec_anual)
contour(prec_anual, add= TRUE)
plot(stn_loc
     ,add = TRUE
     ,col = ifelse(as.integer(substr(stn_loc@data$ID, 1, 1)) == 8
                   ,2
                   ,1))

# totais mensais na bacia
prec_mensal <- zApply(prec_b_mask
                      ,by = month
                      ,fun = sum
                      ,name = ""
                      ,na.rm = TRUE)
prec_mensal <- mask(prec_mensal, r)
names(prec_mensal) <- month.abb
plot(prec_mensal)

# media da prec mensal total na bacia
sp_avg_mly_prec <- cellStats(prec_mensal, mean)
barplot(c(t(sp_avg_mly_prec)), 
names.arg = names(sp_avg_mly_prec), las = 1, xlab = "mês", ylab = "Prec(mm)")
box()

8 Referências

[1] M. F. HUTCHINSON. “Interpolating mean rainfall using thin plate smoothing splines”. In: International journal of geographical information systems 9.4 (Jul. 1995), pp. 385–403. DOI: 10.1080/02693799508902045. URL: http://dx.doi.org/10.1080/02693799508902045.

[2] M. New, M. Hulme and P. Jones. “Representing Twentieth-Century SpaceClimate Variability. Part II: Development of 190196 Monthly Grids of Terrestrial Surface Climate”. In: Journal of Climate 13.13 (Jul. 2000), pp. 2217–2238. DOI: 10.1175/1520-0442(2000)013<2217:rtcstc>2.0.co;2. URL: http://dx.doi.org/10.1175/1520-0442(2000)013<2217:RTCSTC>2.0.CO;2.

[3] A. H. THIESSEN. “PRECIPITATION AVERAGES FOR LARGE AREAS”. In: Mon. Wea. Rev. 39.7 (Jul. 1911), pp. 1082–1089. DOI: 10.1175/1520-0493(1911)39<1082b:pafla>2.0.co;2. URL: http://dx.doi.org/10.1175/1520-0493(1911)39<1082b:PAFLA>2.0.CO;2.