XCO2 - Aquisição de dados de XCO2
Panosso AR; Costa LM; Lima LR
09/04/2021
CO2 Virtual Science Data Environment
O objetivo desse material é apresentar os procedimentos básicos para aquisição de dados do satélite OCO-2 e processamento inicial em R.
Aquisição de dados
1) Acesse o endereço https://co2.jpl.nasa.gov/
2) Acesse o Browse OCO-2 Data em Level 2 Data Set OCO-2.
3) Role a página para baixo e acesse CUSTOMIZE PRODUCT ON BUILD PAGE.
4) No menu à esquerda estarão as 9 categorias para personalizar o banco de dados.
5) Em DATA TYPE selecione OCO-2 Satellite.
6) Em PRODUCTS selecione OCO-2 Full.
7) Em DATA VERSION selecione 10.
8) Em SPATIAL + TEMPORAL selecione Customize Your Spatial + Temporal Coverage.
OPTION 01 Utilize para selecionar a área para aquisição dos dados.
OPTION 02 Utilize para selecionar o período para aquisição dos dados.
9) Em DOWN SAMPLE PRODUCT selecione Yes, I want a Level 3 data product. Altere os valores das células e o passo temporal desejado.
10) Em DATA VARIABLES selecione as variáveis desejadas.
11) Abaixo são apresentadas as opções para os filtros e o tipo de arquivo, selecione CSV FILE.
12) Forneça um endereço de e-mail para onde os links serão direcionados.
13) Acesse o seu e-mail, será enviado uma mensagem com o endereço dos arquivos onde você poderá acompanhar o progresso do processamento de seus dados. Ao final dessa etapa um novo e-mail é enviado informando que os dados estão disponíveis.
14) Acesse o link enviado em seu e-mail e você será direcionado a página.
15) Role a página para baixo e selecione a opção WGET File List.
16) Salve o arquivo .txt na pasta data-raw.
Download e processamento dos dados no R.
Baixando os dados solicitados em CO2 Virtual Science Data Environment
Antes de realizar a leitura, recomendamos a criação da estrutra de diretórios como apresentado abaixo:
pasta data-raw/CVS que receberá o download dos arquivos CSV;
pasta data-raw/txt que receberá os arquivos TXT após a faxina inicial;
pasta data que receberá o arquivos tratados em RDS;
# Download dos dados
# buscar o nome do arquivo contendo o termo CUSTON
%>%` <- magrittr::`%>%`
nome <- dir("data-raw/", pattern = "CUSTOM")
# ler os endereços armazenados no arquivo
enderecos<-read.table(paste0("data-raw/", nome))
n_links <- nrow(enderecos)
n_split <- length(stringr::str_split(enderecos[1,1],"/",simplify = TRUE))
nomes_arquivos_csv <- stringr::str_split_fixed(enderecos[,1],"/",n=Inf)[,n_split]
nomes_arquivos_txt<-sub(".csv",".txt",nomes_arquivos_csv)
### download dos arquivos - aproximadamente 2 GB ao todo
# for(i in 1:n_links){
# download.file(enderecos[i,1],
# paste0("data-raw/CSV",nomes_arquivos_csv[i]),
# mode="wb")
# }Aguarde até todos os arquivos serem baixados. Verifique os arquivos na pasta data-raw/CSV. No repositório do GitHub esses arquivos não serão fornecidos por serem muito grandes.
Pronto, agora temos todos os arquivos na pasta data-raw aguardando a faxina.
Faxina dos arquivos
Os arquivos .csv são, de maneira geral, maiores que 120 \(MB\), portanto, vamos abrir cada arquivo individualmente, filtrar valores faltantes (NA) para a coluna de \(CO_2\) atmosférico (xco2), e salvar em novos arquivos txt na pasta data-raw/txt. Os valos de concentração não observados têm como marcação -999999.
for(i in 1:n_links){
da<-read.csv(paste0("data-raw/CSV/", nomes_arquivos_csv[i]), sep=",")
da <- da %>%
filter(xco2..Moles.Mole...1.. != -999999)
write.table(da,
paste0("data/TXT/",nomes_arquivos_txt[i]),
quote = FALSE,
sep="\t",row.names = FALSE)
}Vamos agora juntar todos os arquivos em um único denominado XCO2_20142020.rds para posterior processamento.
# Lendo os arquivo totais
for(i in 1:n_links){
if(i == 1){
tab<-read.table(paste0("data-raw/TXT/",nomes_arquivos_txt[i]),h=TRUE,sep="\t")
}else{
da<-read.table(paste0("data-raw/TXT/",nomes_arquivos_txt[i]),h=TRUE,sep="\t")
tab<-rbind(tab,da)
}
}
dplyr::glimpse(tab)
readr::write_rds(tab,"data/XCO2_20142020.rds")Processamento
Verificar o banco de dados com a função glimpse() do pacote dplyr.
Carregando os pacotes requeridos
library(magrittr)
library(ggplot2)
library(dplyr)
library(tidyr)
source("../r/graficos.R")
source("../r/my_fun.R")xco2 <- readr::read_rds("../data/XCO2_20142020.rds")
glimpse(xco2)## Rows: 364,529
## Columns: 11
## $ longitude <dbl> -74.58225, -74.58225, -74.58...
## $ longitude_bnds <chr> "-74.70703125:-74.4574652778...
## $ latitude <dbl> -30.22572489, -29.97654828, ...
## $ latitude_bnds <chr> "-30.3503131952:-30.10113658...
## $ time..YYYYMMDDHHMMSS. <dbl> 2.014091e+13, 2.014091e+13, ...
## $ time_bnds..YYYYMMDDHHMMSS. <chr> "20140909000000:201409100000...
## $ altitude..km. <dbl> 3307.8, 3307.8, 3307.8, 3307...
## $ alt_bnds..km. <chr> "0.0:6615.59960938", "0.0:66...
## $ fluorescence_offset_relative_771nm_idp <dbl> 0.012406800, 0.010696600, -0...
## $ fluorescence_offset_relative_757nm_idp <dbl> -3.58630e+00, 8.81219e-02, -...
## $ xco2..Moles.Mole...1.. <dbl> 0.000394333, 0.000395080, 0....Selecionar somente as variáveis necessárias para a sua análise, para isso, vamos criar um outro denominado oco_2, a partir da seleção de longitude, latitude, tempo e a concentração de \(CO_2\) atmosférico:
oco_2<-xco2 %>%
select(longitude,
latitude,
time..YYYYMMDDHHMMSS.,
xco2..Moles.Mole...1..)Agora precisamos transformar os dados de \(CO_2\) para \(ppm\), e criar uma coluna para da data, o dia mês e ano de cada observação, a partir da coluna time..YYYYMMDDHHMMSS.
oco_2 <- oco_2 %>%
dplyr::mutate(
xco2_obs = xco2..Moles.Mole...1.. *1e06,
ano = time..YYYYMMDDHHMMSS.%/%1e10,
mês = time..YYYYMMDDHHMMSS.%/%1e8 %%100,
dia = time..YYYYMMDDHHMMSS.%/%1e6 %%100,
data = as.Date(stringr::str_c(ano,mês,dia,sep="-"))
) %>%
glimpse()## Rows: 364,529
## Columns: 9
## $ longitude <dbl> -74.58225, -74.58225, -74.58225, -74.58225, ...
## $ latitude <dbl> -30.22572489, -29.97654828, -29.72737167, -2...
## $ time..YYYYMMDDHHMMSS. <dbl> 2.014091e+13, 2.014091e+13, 2.014091e+13, 2....
## $ xco2..Moles.Mole...1.. <dbl> 0.000394333, 0.000395080, 0.000394728, 0.000...
## $ xco2_obs <dbl> 394.333, 395.080, 394.728, 394.420, 391.963,...
## $ ano <dbl> 2014, 2014, 2014, 2014, 2014, 2014, 2014, 20...
## $ mês <dbl> 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,...
## $ dia <dbl> 9, 9, 9, 9, 9, 9, 9, 18, 18, 11, 11, 11, 11,...
## $ data <date> 2014-09-09, 2014-09-09, 2014-09-09, 2014-09...Vamos filtrar os dados que estão dentro de uma área desejada, para isso vamos utlizar os polígonos e funções dos pacotes geobr e sp.
Inicialmente, devemos plotar o mapa do Brasil, e ler os polígonos dos municípios de um estado (SP por exmplo), cada estado, cada regiões e do país como um todo.
muni <- geobr::read_municipality(code_muni = "SP",showProgress = FALSE)
estados <- geobr::read_state(code_state = "all",showProgress = FALSE)
regiao <- geobr::read_region(showProgress = FALSE)
br <- geobr::read_country(showProgress = FALSE)
sp <- geobr::read_state(code_state = "SP",showProgress = FALSE)
ms <- geobr::read_state(code_state = "MS",showProgress = FALSE)Plotando as delimitações do país.
br %>%
ggplot() +
geom_sf(fill="#2D3E50", color="#FEBF57",
size=.15, show.legend = FALSE) +
tema_mapa()
Plotando as regiões.
regiao %>%
ggplot() +
geom_sf(fill="#2D3E50", color="#FEBF57",
size=.15, show.legend = FALSE)+
tema_mapa()
Plotando as delimitações do país com os estados.
estados %>%
ggplot() +
geom_sf(fill="#2D3E50", color="#FEBF57",
size=.15, show.legend = FALSE)+
tema_mapa()
Plotando as delimitações do país com os estados específicos.
# plotar os estados
br %>%
ggplot() +
geom_sf(fill="#2D3E50", color="#FEBF57",
size=.15, show.legend = FALSE)+
geom_sf(data=sp,fill="#2D3E50", color="#FEBF57")+
geom_sf(data=ms,fill="#2D3E50", color="#FEBF57")+
tema_mapa()
Todos os munícipios de um estados.
# Fazer o plot
muni %>%
ggplot() +
geom_sf(fill="#2D3E50", color="#FEBF57",
size=.15, show.legend = FALSE)+
tema_mapa()
A partir do mapa base, vamos sobrepor os registros do objeto oco_2.
br %>%
ggplot() +
geom_sf(fill="#2D3E50", color="#FEBF57",
size=.15, show.legend = FALSE)+
tema_mapa() +
geom_point(data=oco_2 %>% filter(ano == 2014) ,
aes(x=longitude,y=latitude),
shape=3,
col="red",
alpha=0.2)
É necessário selecionar somente aqueles pontos dentro do polígono da área. Inicialmente devemos extrair o polígono de cada shape do pacote geobr.
pol_sp <- as.matrix(sp$geom[[1]]) # pol estado de são paulo
pol_ms <- as.matrix(ms$geom[[1]]) # pol estado de ms
pol_br <- as.matrix(br$geom[[1]]) # pol do brasil
pol_br <- pol_br[pol_br[,1]<=-34,] # CORREÇÃO do polígono
pol_br <- pol_br[!((pol_br[,1]>=-38.8 & pol_br[,1]<=-38.6) &
(pol_br[,2]>= -19 & pol_br[,2]<= -16)),]
pol_norte <- as.matrix(regiao$geom[[1]]) # pol região norte
pol_sudeste <- as.matrix(regiao$geom[[3]]) # pol região sudeste
pol_sudeste <- pol_sudeste[pol_sudeste[,1]<=-34,] # CORREÇÃO do polígono
pol_sul <- as.matrix(regiao$geom[[4]]) # pol região sul
pol_centro_oeste <- as.matrix(regiao$geom[[5]]) # pol centro oeste
pol_nordeste <- as.matrix(regiao$geom[[2]]) # pol nordeste
pol_nordeste <- pol_nordeste[pol_nordeste[,1]<=-34,]# CORREÇÃO do polígono
pol_nordeste <- pol_nordeste[!((pol_nordeste[,1]>=-38.7 & pol_nordeste[,1]<=-38.6) &
pol_nordeste[,2]<= -15),]# CORREÇÃO do polígonoSelecionamos os pontos que estão dentro do polígono do estado de São Paulo e Mato Grosso do Sul, no nosso exemplo. Utilizaremos a função point.in.polygon() do pacote sp.
oco_2 <- oco_2 %>%
mutate(
point_in_pol = as.logical(sp::point.in.polygon(point.x = longitude,
point.y = latitude,
pol.x = pol_sp[,1],
pol.y = pol_sp[,2])
|
sp::point.in.polygon(point.x = longitude,
point.y = latitude,
pol.x = pol_ms[,1],
pol.y = pol_ms[,2]))
)Vamos refazer o mapa anterior com os dados dentro do polígono.
br %>%
ggplot() +
geom_sf(fill="#2D3E50", color="#FEBF57",
size=.15, show.legend = FALSE)+
tema_mapa() +
geom_point(data=oco_2 %>% filter(point_in_pol),
aes(x=longitude,y=latitude),
shape=3,
col="red",
alpha=0.2)
Selecionamos os pontos que estão dentro do Brasil.
oco_2 <- oco_2 %>%
mutate(
point_in_pol = as.logical(sp::point.in.polygon(point.x = longitude,
point.y = latitude,
pol.x = pol_br[,1],
pol.y = pol_br[,2])))Vamos refazer o mapa anterior com os dados dentro do polígono.
br %>%
ggplot() +
geom_sf(fill="#2D3E50", color="#FEBF57",
size=.15, show.legend = FALSE)+
tema_mapa() +
geom_point(data=oco_2 %>% filter(ano==2014, point_in_pol),
aes(x=longitude,y=latitude),
shape=3,
col="red",
alpha=0.2)
Observa-se um problema na região nordeste, pontos não alocados, então vamos adicionar a esse mapa os pontos dessa região.
oco_2 <- oco_2 %>%
mutate(
point_in_pol = as.logical(sp::point.in.polygon(point.x = longitude,
point.y = latitude,
pol.x = pol_br[,1],
pol.y = pol_br[,2])
|
sp::point.in.polygon(point.x = longitude,
point.y = latitude,
pol.x = pol_nordeste[,1],
pol.y = pol_nordeste[,2]))
)Vamos refazer o mapa anterior com os dados dentro do polígono.
br %>%
ggplot() +
geom_sf(fill="#2D3E50", color="#FEBF57",
size=.15, show.legend = FALSE)+
tema_mapa() +
geom_point(data=oco_2 %>% filter(point_in_pol),
aes(x=longitude,y=latitude),
shape=3,
col="red",
alpha=0.2)
Vamos olhar os outliers.
oco_2 %>%
filter(point_in_pol) %>%
ggplot(aes(y=xco2..Moles.Mole...1..))+
geom_boxplot(fill="orange",outlier.colour = "red")+
coord_cartesian(xlim=c(-1,1))+
theme_classic()
O padrão de coef é 1,5 e define o comprimento dos bigodes com base no coeficiente e na distância interquartil, com todos os pontos fora dos bigodes considerados outliers.
filtrando os outliers.
oco_2 %>%
filter(point_in_pol, xco2..Moles.Mole...1.. > limites_outlier(oco_2$xco2..Moles.Mole...1..)[1] &
xco2..Moles.Mole...1.. < limites_outlier(oco_2$xco2..Moles.Mole...1..)[2]) %>%
ggplot(aes(y=xco2..Moles.Mole...1..))+
geom_boxplot(fill="orange",outlier.colour = "red")+
coord_cartesian(xlim=c(-1,1))+
theme_classic()
oco_2 %>%
mutate(ano_mes = paste0(ano,mês)) %>%
group_by(ano,mês) %>%
ggplot(aes(y=xco2_obs,x=as.factor(mês),
fill=as.factor(mês)))+
geom_violin(trim = FALSE) +
theme(legend.position = "none")+
facet_wrap(~ano,scales = "free")+
labs(x="Mês",y="XCO2 (ppm)")+
theme_bw()
Observe a emissão de \(CO_2\) ao longo do tempo para esses dados.
oco_2 %>%
filter(point_in_pol) %>%
ggplot(aes(x=data,y=xco2_obs)) +
geom_point(color="blue") +
geom_line(color="red")
Agora devemos retirar a tendência ao longo do tempo, para isso, dentro do período específico, faremos a retirada por meio de um ajuste linear:
oco_2 %>%
arrange(data) %>%
mutate(x= 1:nrow(oco_2)) %>%
ggplot(aes(x=x,y=xco2_obs)) +
geom_point(shape=21,color="black",fill="gray") +
geom_smooth(method = "lm") +
ggpubr::stat_regline_equation(aes(
label = paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~~")))
Extrair os coeficientes \(a\) e \(b\) da análise de regressão linear (\(y = a + bX\)).
d_aux<-oco_2 %>%
arrange(data) %>%
mutate(x= 1:nrow(oco_2)) %>%
filter(point_in_pol) %>%
select(x,xco2_obs)
mod <- lm(d_aux$xco2_obs~d_aux$x)
summary.lm(mod)##
## Call:
## lm(formula = d_aux$xco2_obs ~ d_aux$x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -51.611 -1.547 0.450 2.030 47.133
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.931e+02 1.717e-02 22898.6 <2e-16 ***
## d_aux$x 4.353e-05 8.165e-08 533.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.307 on 147890 degrees of freedom
## Multiple R-squared: 0.6578, Adjusted R-squared: 0.6578
## F-statistic: 2.843e+05 on 1 and 147890 DF, p-value: < 2.2e-16a<-mod$coefficients[1]
b<-mod$coefficients[2]Criando a variável xco2_est a partir da retirada da tendência.
oco_2 <- oco_2 %>%
arrange(data) %>%
mutate(
xco2_est = a + b * as.numeric(data),
delta = xco2_est - xco2_obs,
XCO2 = (a-delta) - (mean(xco2_obs) - a)
) %>%
glimpse()## Rows: 364,529
## Columns: 13
## $ longitude <dbl> -72.58572, -72.33615, -72.33615, -72.08659, ...
## $ latitude <dbl> 6.154060, 5.157354, 5.406530, 3.163941, 3.41...
## $ time..YYYYMMDDHHMMSS. <dbl> 2.014091e+13, 2.014091e+13, 2.014091e+13, 2....
## $ xco2..Moles.Mole...1.. <dbl> 0.000391368, 0.000389822, 0.000388482, 0.000...
## $ xco2_obs <dbl> 391.368, 389.822, 388.482, 390.302, 386.635,...
## $ ano <dbl> 2014, 2014, 2014, 2014, 2014, 2014, 2014, 20...
## $ mês <dbl> 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9,...
## $ dia <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,...
## $ data <date> 2014-09-06, 2014-09-06, 2014-09-06, 2014-09...
## $ point_in_pol <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FA...
## $ xco2_est <dbl> 393.8232, 393.8232, 393.8232, 393.8232, 393....
## $ delta <dbl> 2.45517575, 4.00117575, 5.34117575, 3.521175...
## $ XCO2 <dbl> 382.8391, 381.2931, 379.9531, 381.7731, 378....oco_2 %>%
# filter(point_in_pol) %>%
ggplot(aes(x=data,y=XCO2)) +
geom_point(shape=21,color="black",fill="gray")
Buscando os dados de duas regiões.
# Criando o filtro
oco_2 <- oco_2 %>%
mutate(
point_in_pol_sp_ms = as.logical(sp::point.in.polygon(point.x = longitude,
point.y = latitude,
pol.x = pol_sp[,1],
pol.y = pol_sp[,2])
|
sp::point.in.polygon(point.x = longitude,
point.y = latitude,
pol.x = pol_ms[,1],
pol.y = pol_ms[,2]))
)
# Plotando os dados
ggplot(sp) +
geom_sf(fill="#2D3E50", color="#FEBF57",
size=.15, show.legend = FALSE)+
geom_sf(data=ms,fill="#2D3E50", color="#FEBF57",
size=.15, show.legend = FALSE)+
tema_mapa() +
geom_point(data=oco_2 %>% filter(point_in_pol_sp_ms),
aes(x=longitude,y=latitude,color=xco2_est))
Pronto, agora seu projeto pode continuar…