Pré-Processamento

🛠️ Pré-processameto script

🛠 Preparação dos dados para análise.

library(tidyverse)
library(ggridges)
library(ggpubr)
library(geobr)
library(gstat)
library(vegan)
library(sf)
library(dplyr)
library(stringi)
source("../R/my-function.R")

Definindo estados

my_states <- c("MS","MT","GO","DF")

💨 Entrada com a Base: `emissions-sources.rds`

Essa é a base antiga do Climate TRACE, que compreende o período de 2015 a 2020

emissions_sources <- read_rds("../data/emissions_sources.rds")

💨 Entrada com a Base: `climate-trace-br.xlsx`

Essa é a nova base do Climate TRACE, que compreenderá o período de 2021 a 2023

base_ct_att <- readxl::read_excel("../data-raw/climate-trace-br.xlsx")

💨 Entrada com a Base: `nasa-xco2.rds`

nasa_xco2 <- read_rds("../data/nasa-xco2.rds") |> 
  filter(state %in% my_states)
glimpse(nasa_xco2)

Filtrando os polígonos do municípios

munici_state <- municipality |> 
  filter(abbrev_state %in% my_states)
pol_ms <- states |> filter(abbrev_state == "MS") |> 
  pull(geom) |> pluck(1) |> as.matrix()
pol_mt <- states |> filter(abbrev_state == "MT") |> 
  pull(geom) |> pluck(1) |> as.matrix()
pol_go <- states |> filter(abbrev_state == "GO") |> 
  pull(geom) |> pluck(1) |> as.matrix()
pol_df <- states |> filter(abbrev_state == "DF") |> 
  pull(geom) |> pluck(1) |> as.matrix()

Gerar mapa

my_year = 2021
municipality |> 
  filter(abbrev_state %in% my_states) |> 
  left_join(
    nasa_xco2 |> 
      group_by(year, city_ref) |> 
      filter(year == my_year) |> 
      summarise(
        xco2 = mean(xco2,na.rm=TRUE),
        .groups = "drop"
      ) |> 
      rename(  name_muni = city_ref),
    by = c("name_muni")
  ) |> 
  ggplot()  +
  geom_sf(aes(fill=xco2), color="transparent",
          size=.05, show.legend = TRUE)  +
  geom_point(data = nasa_xco2 |> 
               filter(year==my_year), 
               aes(longitude, latitude, #size = emission,
                   color="red"))+
  theme_bw() +
  theme(
    axis.text.x = element_text(size = rel(.9), color = "black"),
    axis.title.x = element_text(size = rel(1.1), color = "black"),
    axis.text.y = element_text(size = rel(.9), color = "black"),
    axis.title.y = element_text(size = rel(1.1), color = "black"),
    legend.text = element_text(size = rel(1), color = "black"),
    legend.title = element_text(face = 'bold', size = rel(1.2))
  ) +
  labs(fill = 'xco2',
       x = 'Longitude',
       y = 'Latitude') +
  scale_fill_viridis_c()

Criando o grid (malha de pontos) para valores não amostrados - nasa-xco2

# vetores para coordenadas x e y selecionadas da base do IBGE1
# Extrair coordenadas da base nasa_xco2 para definir extensão do grid
x<-nasa_xco2$longitude
y<-nasa_xco2$latitude
dis <- 0.5 # distância (do grid) para adensamento de pontos nos estados
grid_geral <- expand.grid( # Criar malha regular
  X=seq(min(x),max(x),dis), # Combinar cada x e y
  Y=seq(min(y),max(y),dis)) |>
  mutate( #teste ponto-no-polígono, TRUE se (X,Y) caem no polígono
    flag_ms = def_pol(X, Y, pol_ms),
    flag_mt = def_pol(X, Y, pol_mt),
    flag_go = def_pol(X, Y, pol_go),
    flag_df = def_pol(X, Y, pol_df)
  ) |>
  filter(flag_ms | flag_go | flag_mt | flag_df) |> # manter pontos correspondentes ao menos 1 dos limites
  select(-c(flag_ms,flag_mt,flag_go,flag_df)) # remover colunas criadas
plot(grid_geral)
sp::gridded(grid_geral) = ~ X + Y

Interpolação por Krigagem Ordinária

Estimar valores médios de concentração de CO2 entre 2015 e 2023

df_aux <- nasa_xco2 |> 
  filter(year == my_year) |> 
  group_by(longitude, latitude) |> 
  summarise(
    xco2 = mean(xco2,na.rm=TRUE), # média xco2
    .groups = "drop"
  ) |> sample_n(10000) 
sp::coordinates(df_aux) = ~ longitude + latitude # Converter data frame para objeto espacial - atribuir as colunas longitude/latitude como coordenadas. Várias funções de geoestatística do pacote gstat (como variogram() e krige()) não aceitam um data.frame, mas um objeto com coordenadas. Agora, a variável xco2 passa a ser o atributo associado a cada ponto espacial.

form <- xco2 ~ 1 # "~ 1" modelo de média constante, mas desconhecida (somente intercepto - assume média global, sem covariáveis).

vari_exp <- gstat::variogram(form, data = df_aux,
                      cressie = FALSE, #estimador clássico do semivar.
                      cutoff = 1, # distância máxima = 8
                      width = 0.075) # distancia entre pontos
# vari_exp  |>
#   ggplot(aes(x=dist, y=gamma)) +
#   geom_point() +
#   labs(x="lag (º)",
#        y=expression(paste(gamma,"(h)")))

patamar=1.4
alcance=0.2
epepita=0.5
modelo_1 <- fit.variogram(vari_exp,vgm(patamar,"Sph",alcance,epepita))
modelo_2 <- fit.variogram(vari_exp,vgm(patamar,"Exp",alcance,epepita))
modelo_3 <- fit.variogram(vari_exp,vgm(patamar,"Gau",alcance,epepita))
plot_my_models(modelo_1,modelo_2,modelo_3)
modelo <- modelo_2 ## sempre MODIFICAR

✅ Selecionado o melhor modelo, vamos guardá-lo

Este chunk do código tem por objetivo selecionar o melhor modelo dos semivariogramas experimentais para as diferentes variáveis. Assim, pode ser “reutilizado” para qualquer uma das variáveis analisadas, desde que, sejam as variáveis espaciais estimadas pela krigaem ordinária, rode o novo semivariograma desejado e altere:

1. O modelo - ~ Linha 472 (já será alterado ao rodar o novo semivariograma) 2. A variável - ~ Linha 543 e 547 3. O nome do arquivo salvo - ~ Linha 547

O novo semivariograma a ser rodado deve abranger, por exemplo, desde a linha 784 à 904.

As variáveis krigadas foram: XCO2, XCH4, SIF e as do NASA POWER - A base do AppEEARS não foi krigada pois o sensor abrangeu quase que integralmente todos os municípios

# Extraindo dos parâmetros dos modelos
sqr.f1 <- round(attr(modelo_1, "SSErr"), 4)
c01    <- round(modelo_1$psill[[1]], 4)
c0_c11 <- round(sum(modelo_1$psill), 4)
a1     <- round(modelo_1$range[[2]], 2)

sqr.f2 <- round(attr(modelo_2, "SSErr"), 4)
c02    <- round(modelo_2$psill[[1]], 4)
c0_c12 <- round(sum(modelo_2$psill), 4)
a2     <- round(3 * modelo_2$range[[2]], 2)

sqr.f3 <- round(attr(modelo_3, "SSErr"), 4)
c03    <- round(modelo_3$psill[[1]], 4)
c0_c13 <- round(sum(modelo_3$psill), 4)
a3     <- round(modelo_3$range[[2]] * sqrt(3), 2)

# Cálculo dos R² para cada modelo
df_aux <- vari_exp |>
  mutate(
    gamma_m1 = ifelse(dist <= a1,
                      c01 + (c0_c11 - c01) * (1.5*(dist/a1) - 0.5*(dist/a1)^3),
                      c0_c11),

    gamma_m2 = c02 + (c0_c12 - c02) * (1 - exp(-3 * (dist/a2))),

    gamma_m3 = c03 + (c0_c13 - c03) * (1 - exp(-3 * (dist/a3)^2)),

    residuo_total = (gamma - mean(gamma))^2,
    residuo_mod_1 = (gamma - gamma_m1)^2,
    residuo_mod_2 = (gamma - gamma_m2)^2,
    residuo_mod_3 = (gamma - gamma_m3)^2
  ) |>
  summarise(
    r2_1 = (sum(residuo_total) - sum(residuo_mod_1)) / sum(residuo_total),
    r2_2 = (sum(residuo_total) - sum(residuo_mod_2)) / sum(residuo_total),
    r2_3 = (sum(residuo_total) - sum(residuo_mod_3)) / sum(residuo_total)
  )

r21 <- as.vector(round(df_aux[1], 4))
r22 <- as.vector(round(df_aux[2], 4))
r23 <- as.vector(round(df_aux[3], 4))

# Salvando os parâmetros dos melhores modelo
model <- modelo |> slice(2) |> pull(model)
rss <- round(attr(modelo, "SSErr"),4) # Extrai SSErr (Sum of Squared Errors) do modelo
c0 <- round(modelo$psill[[1]],4) # psill[[1]] = o nugget (C0)
c0_c1 <- round(sum(modelo$psill),4) # Soma patamar + efeito pepita = sill total 
a <- ifelse(model == "Gau", round(modelo$range[[2]]*(3^.5),2),
            ifelse(model == "Exp",round(3*modelo$range[[2]],2),
            round(modelo$range[[2]],2))) # alcance ajustado conforme modelo (alcance efetivo)
r2 <- ifelse(model == "Gau", r23,
            ifelse(model == "Exp",r22,
            r21)) |> pluck(1) # R quadrado
tibble( #criando tibble com todos os resultados
  my_states, my_year, model, c0, c0_c1, a, rss, r2
) |> mutate(gde = c0/c0_c1, .after = "a") |>
  rename(state=my_states,year=my_year) |> 
  write_csv(paste0("../output/",paste(my_states,
                                            collapse = "_"), 
                   "-",my_year,"-xco2.csv")) # MODIFICAR variável

ls_csv <- list.files("../output/",full.names = TRUE,pattern = "-xco2.csv") # MODIFICAR
map_df(ls_csv, read_csv) |> 
  writexl::write_xlsx("../output/semivariogram-models-xco2.xlsx") # MODIFICAR

ko_variavel <- krige(formula=form, df_aux, grid_geral, model=modelo,
                     block=c(0.1,0.1),
                     nsim=0,
                     na.action=na.pass,
                     debug.level=-1
)

ko_variavel |> 
  as_tibble() |> 
    ggplot(aes(x=X, y=Y)) +
  geom_tile(aes(fill = var1.pred)) +
  scale_fill_viridis_c(option = "mako") +
  coord_equal() +
  labs(x="Longitude",
       y="Latitude",
       fill="xco2",
       title = my_year) +
  theme_bw()

df_kgr <- ko_variavel |> 
      as_tibble() |> 
      select(-var1.var) |> 
      rename(longitude=X,latitude=Y,xco2=var1.pred)  |> 
      mutate(city_ref = "Other",
             state = ifelse(def_pol(longitude, latitude, pol_ms),"MS",
                            ifelse(def_pol(longitude, latitude, pol_mt),"MT",
                            ifelse(def_pol(longitude, latitude, pol_go),"GO",
                            "DF"))) 
             )
resul <- vector()
estado <- df_kgr$state
for(i in 1:nrow(df_kgr)){
  if(estado[i]!="Other"){
    my_citys_obj <- municipality %>%
      filter(abbrev_state == estado[i])
    n_citys <- nrow(my_citys_obj)
    my_citys_names <- my_citys_obj %>% pull(name_muni)
    resul[i] <- "Other"
    for(j in 1:n_citys){
      pol_city <- my_citys_obj$geom  %>%
        purrr::pluck(j) %>%
        as.matrix()
      if(def_pol(df_kgr$longitude[i],
                 df_kgr$latitude[i],
                 pol_city)){
        resul[i] <- my_citys_names[j]
      }
    }
  }
}
df_kgr$city_ref <- resul

municipality |> 
  filter(abbrev_state %in% my_states) |> 
  left_join(
    nasa_xco2 |> 
      filter(year == my_year) |> 
      select(longitude,latitude,xco2,state,city_ref) |> 
      rbind(
        df_kgr
      ) |> 
      group_by(city_ref) |> 
      summarise(
        xco2 = mean(xco2,na.rm=TRUE),
        .groups = "drop"
      ) |> 
      rename(  name_muni = city_ref),
    by = c("name_muni")
  ) |> 
  mutate(
    xco2 = ifelse(is.na(xco2),median(xco2,na.rm = TRUE),xco2)) |>
  ggplot()  +
  geom_sf(aes(fill=xco2), color="transparent",
          size=.05, show.legend = TRUE)  +
  # geom_point(data = df_kgr, 
  #            aes(longitude, latitude, #size = emission,
  #                color="red")) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = rel(.9), color = "black"),
    axis.title.x = element_text(size = rel(1.1), color = "black"),
    axis.text.y = element_text(size = rel(.9), color = "black"),
    axis.title.y = element_text(size = rel(1.1), color = "black"),
    legend.text = element_text(size = rel(1), color = "black"),
    legend.title = element_text(face = 'bold', size = rel(1.2))
  ) +
  labs(fill = 'xco2',
       x = 'Longitude',
       y = 'Latitude') +
  scale_fill_viridis_c()

Estimativa de XCO2 para o Brasil Central

nasa_xco2_kgr <- map_df(2015:2023,~{
  set.seed(1235)
  df_aux <- nasa_xco2 |>
    group_by(longitude, latitude) |>
    filter(year == .x) |>
    summarise(
      xco2 = mean(xco2,na.rm=TRUE),
      .groups = "drop"
    ) |> sample_n(8000)
  sp::coordinates(df_aux) = ~ longitude + latitude
  form <- xco2 ~ 1
  vari_exp <- gstat::variogram(form, data = df_aux,
                               cressie = FALSE,
                               cutoff = 1, # distância máxima 8
                               width = 0.075) # distancia entre pontos
  vari_exp  |>
    ggplot(aes(x=dist, y=gamma)) +
    geom_point() +
    labs(x="lag (º)",
         y=expression(paste(gamma,"(h)")))
  patamar=1.4
  alcance=0.2
  epepita=0.5
  modelo <- fit.variogram(vari_exp,vgm(patamar,"Exp",alcance,epepita))
  ko_variavel <- krige(formula=form, df_aux, grid_geral, model=modelo,
                       block=c(0.1,0.1),
                       nsim=0,
                       na.action=na.pass,
                       debug.level=-1
  )
  df_kgr <- ko_variavel |>
      as_tibble() |>
      select(-var1.var) |>
      rename(longitude=X,latitude=Y,xco2=var1.pred)  |>
      mutate(city_ref = "Other",
             state = ifelse(def_pol(longitude, latitude, pol_ms),"MS",
                            ifelse(def_pol(longitude, latitude, pol_mt),"MT",
                            ifelse(def_pol(longitude, latitude, pol_go),"GO",
                            "DF")))
             )
  resul <- vector()
  estado <- df_kgr$state
  for(i in 1:nrow(df_kgr)){
    if(estado[i]!="Other"){
      my_citys_obj <- municipality %>%
        filter(abbrev_state == estado[i])
      n_citys <- nrow(my_citys_obj)
      my_citys_names <- my_citys_obj %>% pull(name_muni)
      resul[i] <- "Other"
      for(j in 1:n_citys){
        pol_city <- my_citys_obj$geom  %>%
          purrr::pluck(j) %>%
          as.matrix()
        if(def_pol(df_kgr$longitude[i],
                   df_kgr$latitude[i],
                   pol_city)){
          resul[i] <- my_citys_names[j]
        }
      }
    }
  }
  df_kgr$city_ref <- resul

  df_final <- df_kgr |>
        mutate(year = .x)
  return(df_final)
})
# write_rds(nasa_xco2_kgr,"../data-raw/nasa-xco2-kgr.rds")

nasa_xco2_kgr <- read_rds("../data-raw/nasa-xco2-kgr.rds")

nasa_xco2_bind <- nasa_xco2 |> 
  select(longitude,latitude,xco2,city_ref,state,year) |> 
  rbind(nasa_xco2_kgr)

# Conferindo o municípios de GO que estavam faltando --> Agora estão certos (~230)
# nasa_xco2_bind |> 
#   filter(state == "GO") |> 
#   pull(city_ref) |> unique()

map(2015:2023,~{
  municipality |> 
    filter(abbrev_state %in% my_states) |> 
    left_join(
      nasa_xco2_bind |> 
        filter(year == .x) |> 
        select(longitude,latitude,xco2,state,city_ref) |> 
        group_by(city_ref) |> 
        summarise(
          xco2 = mean(xco2,na.rm=TRUE),
          .groups = "drop"
        ) |> 
        rename(  name_muni = city_ref),
      by = c("name_muni")
    ) |> 
    mutate(
      xco2 = ifelse(is.na(xco2),median(xco2,na.rm = TRUE),xco2)) |>
    ggplot()  +
    geom_sf(aes(fill=xco2), color="transparent",
            size=.05, show.legend = TRUE)  +
    # geom_point(data = df_kgr, 
    #            aes(longitude, latitude, #size = emission,
    #                color="red")) +
    theme_bw() +
    theme(
      axis.text.x = element_text(size = rel(.9), color = "black"),
      axis.title.x = element_text(size = rel(1.1), color = "black"),
      axis.text.y = element_text(size = rel(.9), color = "black"),
      axis.title.y = element_text(size = rel(1.1), color = "black"),
      legend.text = element_text(size = rel(1), color = "black"),
      legend.title = element_text(face = 'bold', size = rel(1.2))
    ) +
    labs(fill = 'xco2',
         x = 'Longitude',
         y = 'Latitude',
         title = .x) +
    scale_fill_viridis_c()})

# Salvar arquivo interpolado
# write_rds(nasa_xco2_bind,"../data/nasa-xco2-bind.rds")

💨 Entrada com a Base: `gosat-xch4.rds`

gosat_xch4 <- read_rds("../data/gosat_xch4.rds") |> 
  filter(state %in% my_states) |> 
  select(-flag_nordeste, -flag_br)

glimpse(gosat_xch4)

gosat_xch4 |> 
  filter(year == 2021) |> 
  ggplot(aes(x=longitude,y=latitude)) +
  geom_point()

Classificando cada ponto em município

Já feito no arquivo da pasta docs

Gerar mapa

gosat_xch4 <- read_rds("../data/gosat_xch4.rds") # os arquivos com underline foram tratados, classificados por município e deverão ser repassados para a pasta data futuramente

my_year = 2020
municipality |> 
  filter(abbrev_state %in% my_states) |> 
  left_join(
    gosat_xch4 |> 
      group_by(year, city_ref) |> 
      summarise(
        xch4 = mean(xch4,na.rm=TRUE),
        .groups = "drop"
      ) |> 
      rename(name_muni = city_ref),
    by = c("name_muni")
  ) |> 
  filter(year == my_year) |> 
  ggplot()  +
  geom_sf(aes(fill=xch4), color="transparent",
          size=.05, show.legend = TRUE)  +
  # geom_point(data = gosat_xch4 |> 
  #              filter(year==my_year), 
  #              aes(longitude, latitude, #size = emission,
  #                  color="red"))+
  theme_bw() +
  theme(
    axis.text.x = element_text(size = rel(.9), color = "black"),
    axis.title.x = element_text(size = rel(1.1), color = "black"),
    axis.text.y = element_text(size = rel(.9), color = "black"),
    axis.title.y = element_text(size = rel(1.1), color = "black"),
    legend.text = element_text(size = rel(1), color = "black"),
    legend.title = element_text(face = 'bold', size = rel(1.2))
  ) +
  labs(fill = 'xch4',
       x = 'Longitude',
       y = 'Latitude') +
  scale_fill_viridis_c()

# unidade partes por blihão (ppb)

Criando o grid (malha de pontos) para valores não amostrados - gosat-xch4

# vetores para coordenadas x e y selecionadas da base do IBGE1
# Extrair coordenadas da base gosat_xch4 para definir extensão do grid
x<-gosat_xch4$longitude
y<-gosat_xch4$latitude
dis <- 0.5 # distância (do grid) para adensamento de pontos nos estados
grid_geral <- expand.grid( # Criar malha regular
  X=seq(min(x),max(x),dis), # Combinar cada x e y
  Y=seq(min(y),max(y),dis)) |>
  mutate( #teste ponto-no-polígono, TRUE se (X,Y) caem no polígono
    flag_ms = def_pol(X, Y, pol_ms),
    flag_mt = def_pol(X, Y, pol_mt),
    flag_go = def_pol(X, Y, pol_go),
    flag_df = def_pol(X, Y, pol_df)
  ) |>
  filter(flag_ms | flag_go | flag_mt | flag_df) |> # manter pontos correspondentes ao menos 1 dos limites
  select(-c(flag_ms,flag_mt,flag_go,flag_df)) # remover colunas criadas
plot(grid_geral)
sp::gridded(grid_geral) = ~ X + Y

Interpolação por Krigagem Ordinária - gosat-xch4

Estimar valores médios de concentração de CH4 entre 2015 e 2021

df_aux <- gosat_xch4 |> 
  filter(year == my_year) |> 
  group_by(longitude, latitude) |> 
  summarise(
    xch4 = mean(xch4,na.rm=TRUE), # média xch4
    .groups = "drop"
  ) 
sp::coordinates(df_aux) = ~ longitude + latitude # Converter data frame para objeto espacial - atribuir as colunas longitude/latitude como coordenadas. Várias funções de geoestatística do pacote gstat (como variogram() e krige()) não aceitam um data frame, mas sim um objeto com coordenadas. Agora, a variável xch4 passa a ser o atributo associado a cada ponto espacial.

form <- xch4 ~ 1 # "~ 1" modelo de média constante, mas desconhecida (somente intercepto - assume média global, sem covariáveis).

vari_exp <- gstat::variogram(form, data = df_aux,
                      cressie = FALSE, #estimador clássico do semivar.
                      cutoff = 32, # distância  
                      width = 2) # distancia entre pontos
vari_exp  |>
  ggplot(aes(x=dist, y=gamma)) +
  geom_point() +
  labs(x="lag (º)",
       y=expression(paste(gamma,"(h)")))

patamar=1250
alcance=33
epepita=197
modelo_1 <- fit.variogram(vari_exp,vgm(patamar,"Sph",alcance,epepita))
modelo_2 <- fit.variogram(vari_exp,vgm(patamar,"Exp",alcance,epepita))
modelo_3 <- fit.variogram(vari_exp,vgm(patamar,"Gau",alcance,epepita))
plot_my_models(modelo_1,modelo_2,modelo_3)
modelo <- modelo_2 ## sempre modificar

ko_variavel <- krige(formula=form, df_aux, grid_geral, model=modelo,
                     block=c(0.1,0.1),
                     nsim=0,
                     na.action=na.pass,
                     debug.level=-1
)

ko_variavel |> 
  as_tibble() |> 
    ggplot(aes(x=X, y=Y)) +
  geom_tile(aes(fill = var1.pred)) +
  scale_fill_viridis_c(option = "mako") +
  coord_equal() +
  labs(x="Longitude",
       y="Latitude",
       fill="xch4",
       title = my_year) +
  theme_bw()

df_kgr <- ko_variavel |> 
      as_tibble() |> 
      select(-var1.var) |> 
      rename(longitude=X,latitude=Y,xch4=var1.pred)  |> 
      mutate(city_ref = "Other",
             state = ifelse(def_pol(longitude, latitude, pol_ms),"MS",
                            ifelse(def_pol(longitude, latitude, pol_mt),"MT",
                            ifelse(def_pol(longitude, latitude, pol_go),"GO",
                            "DF"))) 
             )
resul <- vector()
estado <- df_kgr$state
for(i in 1:nrow(df_kgr)){
  if(estado[i]!="Other"){
    my_citys_obj <- municipality %>%
      filter(abbrev_state == estado[i])
    n_citys <- nrow(my_citys_obj)
    my_citys_names <- my_citys_obj %>% pull(name_muni)
    resul[i] <- "Other"
    for(j in 1:n_citys){
      pol_city <- my_citys_obj$geom  %>%
        purrr::pluck(j) %>%
        as.matrix()
      if(def_pol(df_kgr$longitude[i],
                 df_kgr$latitude[i],
                 pol_city)){
        resul[i] <- my_citys_names[j]
      }
    }
  }
}
df_kgr$city_ref <- resul

municipality |> 
  filter(abbrev_state %in% my_states) |> 
  left_join(
    gosat_xch4 |> 
      filter(year == my_year) |> 
      select(longitude,latitude,xch4,state,city_ref) |> 
      rbind(
        df_kgr
      ) |> 
      group_by(city_ref) |> 
      summarise(
        xch4 = mean(xch4,na.rm=TRUE),
        .groups = "drop"
      ) |> 
      rename(  name_muni = city_ref),
    by = c("name_muni")
  ) |> 
  mutate(
    xch4 = ifelse(is.na(xch4),median(xch4,na.rm = TRUE),xch4)) |>
  ggplot()  +
  geom_sf(aes(fill=xch4), color="transparent",
          size=.05, show.legend = TRUE)  +
  # geom_point(data = df_kgr, 
  #            aes(longitude, latitude, #size = emission,
  #                color="red")) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = rel(.9), color = "black"),
    axis.title.x = element_text(size = rel(1.1), color = "black"),
    axis.text.y = element_text(size = rel(.9), color = "black"),
    axis.title.y = element_text(size = rel(1.1), color = "black"),
    legend.text = element_text(size = rel(1), color = "black"),
    legend.title = element_text(face = 'bold', size = rel(1.2))
  ) +
  labs(fill = 'xch4',
       x = 'Longitude',
       y = 'Latitude') +
  scale_fill_viridis_c()

Estimativa de XCH4 para o Brasil Central

# Warning haviam sido removidos, mas voltaram aparecer - revisar modelos 
gosat_xch4_kgr <- map_df(2015:2021,~{
  set.seed(1235)
  df_aux <- gosat_xch4 |>
    group_by(longitude, latitude) |>
    filter(year == .x) |>
    summarise(
      xch4 = mean(xch4,na.rm=TRUE),
      .groups = "drop"
    ) |> sample_n(1961)
  sp::coordinates(df_aux) = ~ longitude + latitude
  form <- xch4 ~ 1
  vari_exp <- gstat::variogram(form, data = df_aux,
                               cressie = FALSE,
                               cutoff = 32, # distância máxima 
                               width = 2) # distancia entre pontos
  vari_exp  |>
    ggplot(aes(x=dist, y=gamma)) +
    geom_point() +
    labs(x="lag (º)",
         y=expression(paste(gamma,"(h)")))
  patamar=1250
  alcance=33
  epepita=197
  modelo <- fit.variogram(vari_exp,vgm(patamar,"Exp",alcance,epepita))
  ko_variavel <- krige(formula=form, df_aux, grid_geral, model=modelo,
                       block=c(0.1,0.1),
                       nsim=0,
                       na.action=na.pass,
                       debug.level=-1
  )
  df_kgr <- ko_variavel |>
      as_tibble() |>
      select(-var1.var) |>
      rename(longitude=X,latitude=Y,xch4=var1.pred)  |>
      mutate(city_ref = "Other",
             state = ifelse(def_pol(longitude, latitude, pol_ms),"MS",
                            ifelse(def_pol(longitude, latitude, pol_mt),"MT",
                            ifelse(def_pol(longitude, latitude, pol_go),"GO",
                            "DF")))
             )
  resul <- vector()
  estado <- df_kgr$state
  for(i in 1:nrow(df_kgr)){
    if(estado[i]!="Other"){
      my_citys_obj <- municipality %>%
        filter(abbrev_state == estado[i])
      n_citys <- nrow(my_citys_obj)
      my_citys_names <- my_citys_obj %>% pull(name_muni)
      resul[i] <- "Other"
      for(j in 1:n_citys){
        pol_city <- my_citys_obj$geom  %>%
          purrr::pluck(j) %>%
          as.matrix()
        if(def_pol(df_kgr$longitude[i],
                   df_kgr$latitude[i],
                   pol_city)){
          resul[i] <- my_citys_names[j]
        }
      }
    }
  }
  df_kgr$city_ref <- resul

  df_final <- df_kgr |>
        mutate(year = .x)
  return(df_final)
})
# write_rds(gosat_xch4_kgr,"data-raw/gosat-xch4-kgr.rds")

gosat_xch4_kgr <- read_rds("../data-raw/gosat-xch4-kgr.rds")

gosat_xch4_bind <- gosat_xch4 |> 
  select(longitude,latitude,xch4,city_ref,state,year) |> 
  rbind(gosat_xch4_kgr)

map(2015:2021,~{
  municipality |> 
    filter(abbrev_state %in% my_states) |> 
    left_join(
      gosat_xch4_bind |> 
        filter(year == .x) |> 
        select(longitude,latitude,xch4,state,city_ref) |> 
        group_by(city_ref) |> 
        summarise(
          xch4 = mean(xch4,na.rm=TRUE),
          .groups = "drop"
        ) |> 
        rename(name_muni = city_ref),
      by = c("name_muni")
    ) |> 
    mutate(
      xch4 = ifelse(is.na(xch4),median(xch4,na.rm = TRUE),xch4)) |>
    ggplot()  +
    geom_sf(aes(fill=xch4), color="transparent",
            size=.05, show.legend = TRUE)  +
    # geom_point(data = df_kgr, 
    #            aes(longitude, latitude, #size = emission,
    #                color="red")) +
    theme_bw() +
    theme(
      axis.text.x = element_text(size = rel(.9), color = "black"),
      axis.title.x = element_text(size = rel(1.1), color = "black"),
      axis.text.y = element_text(size = rel(.9), color = "black"),
      axis.title.y = element_text(size = rel(1.1), color = "black"),
      legend.text = element_text(size = rel(1), color = "black"),
      legend.title = element_text(face = 'bold', size = rel(1.2))
    ) +
    labs(fill = 'xch4',
         x = 'Longitude',
         y = 'Latitude',
         title = .x) +
    scale_fill_viridis_c()})

# Salvar arquivo interpolado
# write_rds(gosat_xch4_bind,"../data/gosat_xch4_bind.rds")

🍃 Entrada com a Base: `oco2-sif.rds`

oco2_sif <- read_rds("../data/oco2-sif.rds") |> 
  rename(SIF_757 = daily_sif757)

glimpse(oco2_sif)

oco2_sif |>
  filter(year == 2022,
         state %in% my_states) |>  # Brasil Central
  ggplot(aes(x=longitude,y=latitude)) +
  geom_point()

Gerar mapa

my_year = 2020
municipality |> 
  filter(abbrev_state %in% my_states) |> 
  left_join(
    oco2_sif |> 
      group_by(year, city_ref) |> 
      summarise(
        SIF_757 = mean(SIF_757,na.rm=TRUE),
        .groups = "drop"
      ) |> 
      rename(name_muni = city_ref),
    by = c("name_muni")
  ) |> 
  filter(year == my_year) |> 
  ggplot()  +
  geom_sf(aes(fill=SIF_757), color="transparent",
          size=.05, show.legend = TRUE)  +
  # geom_point(data = oco2_sif |> 
  #              filter(year==my_year), 
  #              aes(longitude, latitude, #size = emission,
  #                  color="red"))+
  theme_bw() +
  theme(
    axis.text.x = element_text(size = rel(.9), color = "black"),
    axis.title.x = element_text(size = rel(1.1), color = "black"),
    axis.text.y = element_text(size = rel(.9), color = "black"),
    axis.title.y = element_text(size = rel(1.1), color = "black"),
    legend.text = element_text(size = rel(1), color = "black"),
    legend.title = element_text(face = 'bold', size = rel(1.2))
  ) +
  labs(fill = 'SIF757',
       x = 'Longitude',
       y = 'Latitude') +
  scale_fill_viridis_c()

Criando o grid (malha de pontos) para valores não amostrados - oco2-sif

# vetores para coordenadas x e y selecionadas da base do IBGE1
# Extrair coordenadas da base oco2_sif para definir extensão do grid
x<-oco2_sif$longitude
y<-oco2_sif$latitude
dis <- 0.5 # distância (do grid) para adensamento de pontos nos estados
grid_geral <- expand.grid( # Criar malha regular
  X=seq(min(x),max(x),dis), # Combinar cada x e y
  Y=seq(min(y),max(y),dis)) |>
  mutate( #teste ponto-no-polígono, TRUE se (X,Y) caem no polígono
    flag_ms = def_pol(X, Y, pol_ms),
    flag_mt = def_pol(X, Y, pol_mt),
    flag_go = def_pol(X, Y, pol_go),
    flag_df = def_pol(X, Y, pol_df)
  ) |>
  filter(flag_ms | flag_go | flag_mt | flag_df) |> # manter pontos correspondentes ao menos 1 dos limites
  select(-c(flag_ms,flag_mt,flag_go,flag_df)) # remover colunas criadas
plot(grid_geral)
sp::gridded(grid_geral) = ~ X + Y

Interpolação por Krigagem Ordinária - oco2-sif

Os parâmetros do semivariograma (alcance, patamar e efeito pepita) da sif foram estabelecidos com base na amostragem total dos pontos da base de dados (sem a parte de “sample_n()”).

Estimar valores médios de concentração de SIF entre 2015 e 2023

df_aux <- oco2_sif |> 
  filter(year == my_year) |> 
  group_by(longitude, latitude) |> 
  summarise(
    SIF_757 = mean(SIF_757,na.rm=TRUE), # média sif
    .groups = "drop"
  ) |> sample_n(10000) # amostra
sp::coordinates(df_aux) = ~ longitude + latitude # Converter data frame para objeto espacial - atribuir as colunas longitude/latitude como coordenadas. Várias funções de geoestatística do pacote gstat (como variogram() e krige()) não aceitam um data frame, mas sim um objeto com coordenadas. Agora, a variável sif passa a ser o atributo associado a cada ponto espacial.

form <- SIF_757 ~ 1 # "~ 1" modelo de média constante, mas desconhecida (somente intercepto - assume média global, sem covariáveis).

vari_exp <- gstat::variogram(form, data = df_aux,
                      cressie = FALSE, # estimador clássico do semivar.
                      cutoff = 1.25, # distância semivar
                      width = 0.08) # distancia entre pontos

# vari_exp <- gstat::variogram(form, data = df_aux,
#                       cressie = FALSE, #estimador clássico do semivar.
#                       cutoff = 1, # distância semivar
#                       width = 0.08) # distancia entre pontos

vari_exp  |>
  ggplot(aes(x=dist, y=gamma)) +
  geom_point() +
  labs(x="lag (º)",
       y=expression(paste(gamma,"(h)")))

patamar=0.039
alcance=0.7
epepita=0.0352
modelo_1 <- fit.variogram(vari_exp,vgm(patamar,"Sph",alcance,epepita))
modelo_2 <- fit.variogram(vari_exp,vgm(patamar,"Exp",alcance,epepita))
modelo_3 <- fit.variogram(vari_exp,vgm(patamar,"Gau",alcance,epepita))
plot_my_models(modelo_1,modelo_2,modelo_3)
modelo <- modelo_2 ## sempre modificar

tictoc::tic()
ko_variavel <- krige(formula=form, df_aux, grid_geral, model=modelo,
                     block=c(0.1,0.1),
                     nsim=0,
                     na.action=na.pass,
                     debug.level=-1
)
tictoc::toc()

ko_variavel |> 
  as_tibble() |> 
    ggplot(aes(x=X, y=Y)) +
  geom_tile(aes(fill = var1.pred)) +
  scale_fill_viridis_c(option = "mako") +
  coord_equal() +
  labs(x="Longitude",
       y="Latitude",
       fill="SIF",
       title = my_year) +
  theme_bw()

df_kgr <- ko_variavel |> 
      as_tibble() |> 
      select(-var1.var) |> 
      rename(longitude=X,latitude=Y,SIF_757=var1.pred)  |> 
      mutate(city_ref = "Other",
             state = ifelse(def_pol(longitude, latitude, pol_ms),"MS",
                            ifelse(def_pol(longitude, latitude, pol_mt),"MT",
                            ifelse(def_pol(longitude, latitude, pol_go),"GO",
                            "DF"))) 
             )
resul <- vector()
estado <- df_kgr$state
for(i in 1:nrow(df_kgr)){
  if(estado[i]!="Other"){
    my_citys_obj <- municipality %>%
      filter(abbrev_state == estado[i])
    n_citys <- nrow(my_citys_obj)
    my_citys_names <- my_citys_obj %>% pull(name_muni)
    resul[i] <- "Other"
    for(j in 1:n_citys){
      pol_city <- my_citys_obj$geom  %>%
        purrr::pluck(j) %>%
        as.matrix()
      if(def_pol(df_kgr$longitude[i],
                 df_kgr$latitude[i],
                 pol_city)){
        resul[i] <- my_citys_names[j]
      }
    }
  }
}
df_kgr$city_ref <- resul

municipality |> 
  filter(abbrev_state %in% my_states) |> 
  left_join(
    oco2_sif |> 
      filter(year == my_year) |> 
      select(longitude,latitude,SIF_757,state,city_ref) |> 
      rbind(
        df_kgr
      ) |> 
      group_by(city_ref) |> 
      summarise(
        SIF_757 = mean(SIF_757,na.rm=TRUE),
        .groups = "drop"
      ) |> 
      rename(name_muni = city_ref),
    by = c("name_muni")
  ) |> 
  mutate(
    SIF_757 = ifelse(is.na(SIF_757),median(SIF_757,na.rm = TRUE),SIF_757)) |>
  ggplot()  +
  geom_sf(aes(fill=SIF_757), color="transparent",
          size=.05, show.legend = TRUE)  +
  # geom_point(data = df_kgr, 
  #            aes(longitude, latitude, #size = emission,
  #                color="red")) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = rel(.9), color = "black"),
    axis.title.x = element_text(size = rel(1.1), color = "black"),
    axis.text.y = element_text(size = rel(.9), color = "black"),
    axis.title.y = element_text(size = rel(1.1), color = "black"),
    legend.text = element_text(size = rel(1), color = "black"),
    legend.title = element_text(face = 'bold', size = rel(1.2))
  ) +
  labs(fill = 'SIF_757',
       x = 'Longitude',
       y = 'Latitude') +
  scale_fill_viridis_c()

Estimativa de SIF para o Brasil Central

oco2_sif_kgr <- map_df(2015:2023,~{
  set.seed(1235)
  df_aux <- oco2_sif |>
    group_by(longitude, latitude) |>
    filter(year == .x) |>
    summarise(
      SIF_757 = mean(SIF_757,na.rm=TRUE),
      .groups = "drop"
    ) |> sample_n(10000)
  sp::coordinates(df_aux) = ~ longitude + latitude
  form <- SIF_757 ~ 1
  vari_exp <- gstat::variogram(form, data = df_aux,
                               cressie = FALSE,
                               cutoff = 1.25, # distância máxima
                               width = 0.08) # distancia entre pontos
  vari_exp  |>
    ggplot(aes(x=dist, y=gamma)) +
    geom_point() +
    labs(x="lag (º)",
         y=expression(paste(gamma,"(h)")))
  patamar=0.039
  alcance=0.7
  epepita=0.0352
  modelo <- fit.variogram(vari_exp,vgm(patamar,"Sph",alcance,epepita))
  ko_variavel <- krige(formula=form, df_aux, grid_geral, model=modelo,
                       block=c(0.1,0.1),
                       nsim=0,
                       na.action=na.pass,
                       debug.level=-1
  )
  df_kgr <- ko_variavel |>
      as_tibble() |>
      select(-var1.var) |>
      rename(longitude=X,latitude=Y,SIF_757=var1.pred)  |>
      mutate(city_ref = "Other",
             state = ifelse(def_pol(longitude, latitude, pol_ms),"MS",
                            ifelse(def_pol(longitude, latitude, pol_mt),"MT",
                            ifelse(def_pol(longitude, latitude, pol_go),"GO",
                            "DF")))
             )
  resul <- vector()
  estado <- df_kgr$state
  for(i in 1:nrow(df_kgr)){
    if(estado[i]!="Other"){
      my_citys_obj <- municipality %>%
        filter(abbrev_state == estado[i])
      n_citys <- nrow(my_citys_obj)
      my_citys_names <- my_citys_obj %>% pull(name_muni)
      resul[i] <- "Other"
      for(j in 1:n_citys){
        pol_city <- my_citys_obj$geom  %>%
          purrr::pluck(j) %>%
          as.matrix()
        if(def_pol(df_kgr$longitude[i],
                   df_kgr$latitude[i],
                   pol_city)){
          resul[i] <- my_citys_names[j]
        }
      }
    }
  }
  df_kgr$city_ref <- resul

  df_final <- df_kgr |>
        mutate(year = .x)
  return(df_final)
})
# write_rds(oco2_sif_kgr,"../data-raw/oco2-sif-kgr.rds")

oco2_sif_kgr <- read_rds("../data-raw/oco2-sif-kgr.rds")

oco2_sif_bind <- oco2_sif |> 
  select(longitude,latitude,SIF_757,city_ref,state,year) |> 
  rbind(oco2_sif_kgr)

map(2015:2023,~{
  municipality |> 
    filter(abbrev_state %in% my_states) |> 
    left_join(
      oco2_sif_bind |> 
        filter(year == .x) |> 
        select(longitude,latitude,SIF_757,state,city_ref) |> 
        group_by(city_ref) |> 
        summarise(
          SIF_757 = mean(SIF_757,na.rm=TRUE),
          .groups = "drop"
        ) |> 
        rename(name_muni = city_ref),
      by = c("name_muni")
    ) |> 
    mutate(
      SIF_757 = ifelse(is.na(SIF_757),median(SIF_757,na.rm = TRUE),SIF_757)) |>
    ggplot()  +
    geom_sf(aes(fill=SIF_757), color="transparent",
            size=.05, show.legend = TRUE)  +
    # geom_point(data = df_kgr, 
    #            aes(longitude, latitude, #size = emission,
    #                color="red")) +
    theme_bw() +
    theme(
      axis.text.x = element_text(size = rel(.9), color = "black"),
      axis.title.x = element_text(size = rel(1.1), color = "black"),
      axis.text.y = element_text(size = rel(.9), color = "black"),
      axis.title.y = element_text(size = rel(1.1), color = "black"),
      legend.text = element_text(size = rel(1), color = "black"),
      legend.title = element_text(face = 'bold', size = rel(1.2))
    ) +
    labs(fill = 'SIF_757',
         x = 'Longitude',
         y = 'Latitude',
         title = .x) +
    scale_fill_viridis_c(limits = c(0,.5))})

# write_rds(oco2_sif_bind,"../data/oco2-sif-bind.rds")

🍃 Entrada com a Base: `appeears-modis.rds`

# original archive "appeears-modis.rds" = 8,7mb
# appeears data requisited by API (project fapesp)

appeears_modis <- read_rds("../data/appeears_modis.rds")

glimpse(appeears_modis)

appeears_modis |> 
  filter(year == 2020) |>
  ggplot(aes(x=lon,y=lat)) +
  geom_point()

Gerar mapas das variáveis - appeears

# my_year = 2020

# Definir a variável
variavel <- "media_lai" # <-- mudar

# Definir os limites esperados para cada variável - padronizar visualização
limits_map <- list(
  media_ndvi = c(-1, 1),
  media_evi  = c(-1, 1),
  media_fpar = c(0, 1),
  media_lai  = c(0, 7), # máximo valor = 6,85 (após remoção de outliers)
  media_et   = NULL     # sem limite fixo
)

# Limite para a variável atual
limites <- limits_map[[variavel]]

# .data[[variavel]] # Para reconhecer o objeto "variavel" como nome de coluna, e não como string

# Definir os anos
anos <- if (variavel == "media_et") 2021:2023 else 2015:2023

# Aplicando função para geração dos mapas de cada ano
# for (my_year in anos) {
# 
#   # Gerando mapas
# p <- municipality |>
#   filter(abbrev_state %in% my_states) |>
#   left_join(
#     appeears_modis |>
#       group_by(year, city_ref) |>
#       reframe( #retornar somente coluna da variavel
#         !!variavel := .data[[variavel]],na.rm=TRUE, # VERIFICAR - municipios em cinza ao cacular a media
#         .groups = "drop"
#       ) |>
#       filter(.data[[variavel]] <= 200) |> # Filtrar outliers absurdos
#       rename(name_muni = city_ref),
#     by = c("name_muni")
#   ) |>
#   filter(year == my_year) |>
#   ggplot()  +
#   geom_sf(aes(fill=.data[[variavel]]), color="transparent",
#           size=.05, show.legend = TRUE)  +
#   # geom_point(data = appeears_modis |>
#   #              filter(year==my_year),
#   #              aes(lon, lat, #size = emission,
#   #                  color="red"))+
#   theme_bw() +
#   theme(
#     axis.text.x = element_text(size = rel(.9), color = "black"),
#     axis.title.x = element_text(size = rel(1.1), color = "black"),
#     axis.text.y = element_text(size = rel(.9), color = "black"),
#     axis.title.y = element_text(size = rel(1.1), color = "black"),
#     legend.text = element_text(size = rel(1), color = "black"),
#     legend.title = element_text(face = 'bold', size = rel(1.2))
#   ) +
#   labs(fill = variavel,
#        x = 'Longitude',
#        y = 'Latitude') +
#   scale_fill_viridis_c(limits = limites)
# 
#   print(p)
# 
#     ggsave(filename = paste0("../img/mapa_", variavel, "_", my_year, ".png"),
#          plot = p)
# } # arquivos gerados excluídos devido memória cheia

Finalizar tratamento para posterior incorporação

# último tratamento da base prodes
# normalizando nomes dos municipios e padronizando nomes de colunas
modis_final <- read_rds("../data/appeears_modis.rds")

modis_final <- modis_final |>
  filter(media_fpar < 250, media_lai < 250) |> 
  mutate(
    city_ref = stri_trans_general(tolower(city_ref), "Latin-ASCII"),
    city_ref = trimws(city_ref)
  ) |> 
  group_by(year, state, city_ref) |> 
  summarise(
    media_fpar = mean(media_fpar, na.rm = TRUE),
    media_lai = mean(media_lai, na.rm = TRUE),
    media_evi = mean(media_evi, na.rm = TRUE),
    media_ndvi = mean(media_ndvi, na.rm = TRUE),
    media_et = mean(media_et, na.rm = TRUE)
  )

# write_rds(modis_final, "../data/appeears_modis_final.rds")

Não foi feita a interpolação para dados do appeears, pois todos os municípios estão presentes. Os dados foram agrupados por média.

💨 Entrada com a Base: `nasa-power.rds`

# original archive "nasa-power.rds" = 174,6mb

nasa_power <- read_rds("../data/nasa_power.rds") |> 
  rename(
    Radiacao = allsky_sfc_sw_dwn,
    Temperatura = t2m,
    Precipitacao = prectotcorr,
    Umidade = rh2m,
    Vento = ws2m,
    Pressao = ps
  )

glimpse(nasa_power)

# Temperatura (t2m), precipitação (prectotcorr), radiação solar (allsky) e umidade relativa a 2 m (rh2m), velocidade do vento a 2 metros (ws2m) e pressão na superfície (ps).

nasa_power |> 
  filter(year == 2023) |> 
  ggplot(aes(x=lon,y=lat)) +
  geom_point()

Classificando cada ponto em município

Já feito no script de faxina “nasa-power.Rmd”

Gerar mapa

my_year <- 2020

# Atribuir ao objeto "variavel" a variável climática da base do nasa power que iremos trabalhar
  # Radiacao 
  # Temperatura
  # Precipitacao
  # Umidade 
  # Vento
  # Pressao

# Definir variável
variavel <- "Temperatura" #mudar

# Estabelecer anos
anos <- 2015:2023

# Aplicando função para geração dos mapas de cada ano
for (my_year in anos) {

n <- municipality |> 
  filter(abbrev_state %in% my_states) |> 
  left_join(
    nasa_power |> 
      group_by(year, city_ref) |>
      summarise(
        !!variavel := mean(.data[[variavel]],na.rm=TRUE), #Criar coluna com mesmo nome da variavel
        .groups = "drop"
      ) |> 
      rename(name_muni = city_ref),
    by = c("name_muni")
  ) |> 
  filter(year == my_year) |> 
  ggplot()  +
  geom_sf(aes(fill=.data[[variavel]]), color="transparent",
          size=.05, show.legend = TRUE)  +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = rel(.9), color = "black"),
    axis.title.x = element_text(size = rel(1.1), color = "black"),
    axis.text.y = element_text(size = rel(.9), color = "black"),
    axis.title.y = element_text(size = rel(1.1), color = "black"),
    legend.text = element_text(size = rel(1), color = "black"),
    legend.title = element_text(face = 'bold', size = rel(1.2))
  ) +
  labs(fill = variavel,
       x = 'Longitude',
       y = 'Latitude') +
  scale_fill_viridis_c()

  print(n)
  
  # ggsave(filename = paste0("img/mapa_nasapower", variavel, "_", my_year, ".png"),
  #        plot = n)
}

Criando o grid (malha de pontos) para valores não amostrados - nasa-power

# vetores para coordenadas x e y selecionadas da base do IBGE1
# Extrair coordenadas da base nasa_power para definir extensão do grid
x<-nasa_power$lon
y<-nasa_power$lat
dis <- 0.5 # distância (do grid) para adensamento de pontos nos estados
grid_geral <- expand.grid( # Criar malha regular
  X=seq(min(x),max(x),dis), # Combinar cada x e y
  Y=seq(min(y),max(y),dis)) |>
  mutate( #teste ponto-no-polígono, TRUE se (X,Y) caem no polígono
    flag_ms = def_pol(X, Y, pol_ms),
    flag_mt = def_pol(X, Y, pol_mt),
    flag_go = def_pol(X, Y, pol_go),
    flag_df = def_pol(X, Y, pol_df)
  ) |>
  filter(flag_ms | flag_go | flag_mt | flag_df) |> # manter pontos correspondentes ao menos 1 dos limites
  select(-c(flag_ms,flag_mt,flag_go,flag_df)) # remover colunas criadas
plot(grid_geral)
sp::gridded(grid_geral) = ~ X + Y

Interpolação por Krigagem Ordinária - nasa-power

Parâmetros testados para contruir o semivariograma, com base na variável:

Temperatura cutoff = 5,
width = 0.09 patamar=2.2 alcance=4.5 epepita=0.2

Radiacao Modelo Gaussiano teve menro R2, mas melhor interpolação cutoff = 12.5 width = 0.5 patamar = 1.4 alcance = 11.25 epepita = 0.1

Precipitacao cutoff = 6,
width = 0.5 patamar= 1.35 alcance= 15 epepita= 0.01

Umidade Muitos avisos de warning Warning in fit.variogram(vari_exp, vgm(patamar, “Sph”, alcance, epepita)) :No convergence after 200 iterations: try different initial values? cutoff = 12,
width = 0.5 patamar= 60 alcance= 11 epepita= 0.1

Vento 3 warning cutoff = 5,
width = 0.08 patamar=.2 alcance=1.5 epepita=0

Pressao cutoff = 7, width = 0.35 patamar=4.2 alcance=5 epepita=.4

df_aux <- nasa_power |> 
  filter(year %in% my_year) |> 
  group_by(lon, lat) |> 
  summarise(
    !!variavel := mean(.data[[variavel]],na.rm=TRUE), # média temperatura
    .groups = "drop"
  ) # |> sample_n(546) # tamanho máximo dos dados
sp::coordinates(df_aux) = ~ lon + lat # Converter data frame para objeto espacial - atribuir as colunas longitude/latitude como coordenadas. Várias funções de geoestatística do pacote gstat (como variogram() e krige()) não aceitam um data frame, mas sim um objeto com coordenadas. Agora, a variável passa a ser o atributo associado a cada ponto espacial.

# Converter string em fórmula para atribuir a variavel preestabelecida
form <- as.formula(paste(variavel, "~ 1")) 

vari_exp <- gstat::variogram(form, data = df_aux,
                      cressie = FALSE, #estimador clássico do semivar.
                      cutoff = 7, # distância  
                      width = 0.35) # distancia entre pontos
vari_exp  |>
  ggplot(aes(x=dist, y=gamma)) +
  geom_point() +
  labs(x="lag (º)",
       y=expression(paste(gamma,"(h)")))

patamar=4.2
alcance=5
epepita=.4
modelo_1 <- fit.variogram(vari_exp,vgm(patamar,"Sph",alcance,epepita))
modelo_2 <- fit.variogram(vari_exp,vgm(patamar,"Exp",alcance,epepita))
modelo_3 <- fit.variogram(vari_exp,vgm(patamar,"Gau",alcance,epepita))
plot_my_models(modelo_1,modelo_2,modelo_3)
modelo <- modelo_1 ## sempre modificar

ko_variavel <- krige(formula=form, df_aux, grid_geral, model=modelo,
                     block=c(0.1,0.1),
                     nsim=0,
                     na.action=na.pass,
                     debug.level=-1
)

ko_variavel |> 
  as_tibble() |> 
    ggplot(aes(x=X, y=Y)) +
  geom_tile(aes(fill = var1.pred)) +
  scale_fill_viridis_c(option = "mako") +
  coord_equal() +
  labs(x="Longitude",
       y="Latitude",
       fill=variavel,
       title = my_year) +
  theme_bw()

df_kgr <- ko_variavel |> 
      as_tibble() |> 
      select(-var1.var) |> 
      rename(lon=X,lat=Y,!!variavel :=var1.pred)  |> 
      mutate(city_ref = "Other",
             state = ifelse(def_pol(lon, lat, pol_ms),"MS",
                            ifelse(def_pol(lon, lat, pol_mt),"MT",
                            ifelse(def_pol(lon, lat, pol_go),"GO",
                            "DF"))) 
             )
resul <- vector()
estado <- df_kgr$state
for(i in 1:nrow(df_kgr)){
  if(estado[i]!="Other"){
    my_citys_obj <- municipality %>%
      filter(abbrev_state == estado[i])
    n_citys <- nrow(my_citys_obj)
    my_citys_names <- my_citys_obj %>% pull(name_muni)
    resul[i] <- "Other"
    for(j in 1:n_citys){
      pol_city <- my_citys_obj$geom  %>%
        purrr::pluck(j) %>%
        as.matrix()
      if(def_pol(df_kgr$lon[i],
                 df_kgr$lat[i],
                 pol_city)){
        resul[i] <- my_citys_names[j]
      }
    }
  }
}
df_kgr$city_ref <- resul

municipality |> 
  filter(abbrev_state %in% my_states) |> 
  left_join(
    nasa_power |> 
      filter(year == my_year) |> 
      select(lon,lat,variavel,state,city_ref) |> 
      rbind(
        df_kgr
      ) |> 
      group_by(city_ref) |> 
      summarise(
        !!variavel := mean(.data[[variavel]],na.rm=TRUE),
        .groups = "drop"
      ) |> 
      rename(  name_muni = city_ref),
    by = c("name_muni")
  ) |> 
  mutate(
    !!variavel := ifelse(is.na(.data[[variavel]]),median(.data[[variavel]],na.rm = TRUE),.data[[variavel]])) |>
  ggplot()  +
  geom_sf(aes(fill=.data[[variavel]]), color="transparent",
          size=.05, show.legend = TRUE)  +
  # geom_point(data = df_kgr, 
  #            aes(lon, lat, #size = emission,
  #                color="red")) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = rel(.9), color = "black"),
    axis.title.x = element_text(size = rel(1.1), color = "black"),
    axis.text.y = element_text(size = rel(.9), color = "black"),
    axis.title.y = element_text(size = rel(1.1), color = "black"),
    legend.text = element_text(size = rel(1), color = "black"),
    legend.title = element_text(face = 'bold', size = rel(1.2))
  ) +
  labs(fill = variavel,
       x = 'Longitude',
       y = 'Latitude') +
  scale_fill_viridis_c()

Estimativa da variável para o Brasil Central - nasa-power

nome_variavel <- paste0("nasa_power_",variavel,"_kgr")

nome_variavel <- map_df(2015:2023,~{
  set.seed(1235)
  df_aux <- nasa_power |>
    group_by(lon, lat) |>
    filter(year == .x) |>
    summarise(
      !!variavel := mean(.data[[variavel]],na.rm=TRUE),
      .groups = "drop"
    ) 
  sp::coordinates(df_aux) = ~ lon + lat
  form <- as.formula(paste(variavel, "~ 1")) 
  vari_exp <- gstat::variogram(form, data = df_aux,
                               cressie = FALSE,
                               cutoff = 7, # !! mudar
                               width = 0.35) # !! mudar
  vari_exp  |>
    ggplot(aes(x=dist, y=gamma)) +
    geom_point() +
    labs(x="lag (º)",
         y=expression(paste(gamma,"(h)")))
    patamar=4.2 #!! mudar
    alcance=5 #!! mudar
    epepita=.4 #!! mudar
  modelo <- fit.variogram(vari_exp,vgm(patamar,"Exp",alcance,epepita))  # !! mudar
  ko_variavel <- krige(formula=form, df_aux, grid_geral, model=modelo,
                       block=c(0.1,0.1),
                       nsim=0,
                       na.action=na.pass,
                       debug.level=-1
  )
  df_kgr <- ko_variavel |>
      as_tibble() |>
      select(-var1.var) |>
      rename(lon=X,lat=Y,!!variavel :=var1.pred)  |>
      mutate(city_ref = "Other",
             state = ifelse(def_pol(lon, lat, pol_ms),"MS",
                            ifelse(def_pol(lon, lat, pol_mt),"MT",
                            ifelse(def_pol(lon, lat, pol_go),"GO",
                            "DF")))
             )
  resul <- vector()
  estado <- df_kgr$state
  for(i in 1:nrow(df_kgr)){
    if(estado[i]!="Other"){
      my_citys_obj <- municipality %>%
        filter(abbrev_state == estado[i])
      n_citys <- nrow(my_citys_obj)
      my_citys_names <- my_citys_obj %>% pull(name_muni)
      resul[i] <- "Other"
      for(j in 1:n_citys){
        pol_city <- my_citys_obj$geom  %>%
          purrr::pluck(j) %>%
          as.matrix()
        if(def_pol(df_kgr$lon[i],
                   df_kgr$lat[i],
                   pol_city)){
          resul[i] <- my_citys_names[j]
        }
      }
    }
  }
  df_kgr$city_ref <- resul

  df_final <- df_kgr |>
        mutate(year = .x)
  return(df_final)
})
# write_rds(nome_variavel,paste0("../data-raw/nasa-power-",variavel,"-kgr.rds"))

# Rodar somente quando necessário!!!!
nasa_power_kgr <- read_rds(paste0("../data-raw/nasa-power-",variavel,"-kgr.rds"))

nasa_power_bind <- nasa_power |> 
  select(lon,lat,.data[[variavel]],city_ref,state,year) |> 
  rbind(nasa_power_kgr)

map(2015:2023,~{
  municipality |> 
    filter(abbrev_state %in% my_states) |> 
    left_join(
      nasa_power_bind |> 
        filter(year == .x) |> 
        select(lon,lat,.data[[variavel]],state,city_ref) |> 
        group_by(city_ref) |> 
        summarise(
          !!variavel := mean(.data[[variavel]],na.rm=TRUE),
          .groups = "drop"
        ) |> 
        rename(name_muni = city_ref),
      by = c("name_muni")
    ) |> 
    mutate(
      !!variavel := ifelse(is.na(.data[[variavel]]),median(.data[[variavel]],na.rm = TRUE),.data[[variavel]])) |>
    ggplot()  +
    geom_sf(aes(fill=.data[[variavel]]), color="transparent",
            size=.05, show.legend = TRUE)  +
    # geom_point(data = df_kgr, 
    #            aes(longitude, latitude, #size = emission,
    #                color="red")) +
    theme_bw() +
    theme(
      axis.text.x = element_text(size = rel(.9), color = "black"),
      axis.title.x = element_text(size = rel(1.1), color = "black"),
      axis.text.y = element_text(size = rel(.9), color = "black"),
      axis.title.y = element_text(size = rel(1.1), color = "black"),
      legend.text = element_text(size = rel(1), color = "black"),
      legend.title = element_text(face = 'bold', size = rel(1.2))
    ) +
    labs(fill = variavel,
         x = 'Longitude',
         y = 'Latitude',
         title = .x) +
    scale_fill_viridis_c()})

# write_rds(nasa_power_bind,paste0("../data/nasa-power-",variavel,"-bind.rds"))

🪓 Entrada com a Base: `prodes-deforestation.rds`

# original archive "prodes-deforestation.rds" = 2.2gb

prodes_deforestation <- read_rds("../data/prodes_deforestation.rds")
  # sample_n(1000000)
  # "state" column previously filtred 

glimpse(prodes_deforestation)

prodes_deforestation |> 
  sample_n(100000) |> 
  filter(categorie == 15) |> 
  ggplot(aes(x=x,y=y)) +
  geom_point()

Gerar mapa

# Alocando intensidade com base nas observações de desmatamento no ano
df_prodes <- prodes_deforestation |> 
  group_by(categorie, state, muni) |> 
  count() |> 
  group_by(categorie) |> 
  mutate(
    percent = n/sum(n)*100
  )

# df_prodes |> 
#   filter(state == 'MS')

# Gerando o mapa de intensidade de desmatamento
map(15:22,~{
municipality |> 
    filter(abbrev_state %in% my_states) |> 
    left_join(
      df_prodes |> 
        filter(categorie == .x) |>  # mudar anos
        rename(name_muni = muni),
      by = c("name_muni")
    ) |> 
  ggplot() +
   geom_sf(aes(fill=percent), color="transparent",
            size=.05, show.legend = TRUE)  +
    # geom_point(data = df_kgr, 
    #            aes(longitude, latitude, #size = emission,
    #                color="red")) +
    theme_bw() +
    theme(
      axis.text.x = element_text(size = rel(.9), color = "black"),
      axis.title.x = element_text(size = rel(1.1), color = "black"),
      axis.text.y = element_text(size = rel(.9), color = "black"),
      axis.title.y = element_text(size = rel(1.1), color = "black"),
      legend.text = element_text(size = rel(1), color = "black"),
      legend.title = element_text(face = 'bold', size = rel(1.2))
    ) +
    labs(fill = 'Intensidade',
         x = 'Longitude',
         y = 'Latitude',
         title = 'Desmatamento') +
    scale_fill_viridis_c()})

# Município em destaque no mapa é Corumbá (para anos como 2022)
 # municipality |> 
 #    filter(abbrev_state %in% my_states) |> 
 #   filter(
 #     name_muni == "Corumbá"
 #   ) |> 
 #    left_join(
 #      df_prodes |> 
 #        filter(categorie == 16) |> 
 #        rename(name_muni = muni),
 #      by = c("name_muni")
 #    )

Finalizar tratamento para posterior incorporação

# último tratamento da base prodes
# normalizando nomes dos municipios e padronizando nomes de colunas
prodes_final <- read_rds("../data/prodes_deforestation.rds")
prodes_final <- prodes_final |> 
  filter(categorie %in% 15:22) |> 
  rename(year = categorie,
         city_ref = muni) |> 
  mutate(
    year = 2000 + as.numeric(year)
    ) |>
  group_by(state, city_ref, year) |> 
  count() |> 
  group_by(year) |> 
  mutate(
    percent = n/sum(n)*100
  ) |> 
  rename(desmatamento = percent) |> # destamento proporcional ao número de ocorrências em relação ao total de registros daquele ano.
  mutate(
    city_ref = stri_trans_general(tolower(city_ref), "Latin-ASCII"),
    city_ref = trimws(city_ref)
  ) |> 
  group_by(year, city_ref, state) |> 
  summarise(desmatamento = mean(desmatamento, na.rm = TRUE))
  
# write_rds(prodes_final, "../data/prodes_final.rds")

🔥 Entrada com a Base: `deter_queimadas.rds`

É importante ressaltar que, embora diversos estados estejam contemplados nas observações, estes não apresentam informações referentes a queimadas, somente desmatamento ou outras.

Apenas parte do estado do Mato Grosso apresenta observações de queimadas

Disponibilidade de dados de cicatriz de queimada para biomas como Pantanal(só 2023) e Amazônia (https://terrabrasilis.dpi.inpe.br/downloads)

# original archive "deter-queimadas.rds" = 122,1mb

deter_queimadas <- read_rds("../data/deter_queimadas.rds")

# Somente MT apresenta dados de queimadas
# deter_queimadas |>
#   filter(UF == "MT") |> # alternar estados
#   pull(CLASSNAME) |> unique() |>
#   glimpse()

deter_queimadas |> 
  # pull(ANO) |> unique()
  filter(ANO == 2020,
         CLASSNAME == "CICATRIZ_DE_QUEIMADA") |> 
  ggplot(aes(x=longitude,y=latitude)) +
  geom_point()

Padronizando nomes dos municípios da base DETER com os do pacote geobr

Este chunk tem como objetivo padronizar para posterior incorporação, necessária para evitar perda de informação durante a geração dos mapas

library(stringi)  # para normalização de strings

# Alocando área de queimadas
df_deter <- deter_queimadas |> 
  filter(UF %in% my_states) |> 
  rename(name_muni = MUNICIPALI) |> 
  select(name_muni, ANO, AREAMUNKM) |> 
  arrange(name_muni) |> 
  group_by(ANO, name_muni) |> 
  summarise(
    area_queimada_muni_km = sum(AREAMUNKM, na.rm = TRUE),
    .groups = "drop"
  ) |> 
  mutate(
    name_muni = stri_trans_general(tolower(name_muni), "Latin-ASCII"),
    name_muni = trimws(name_muni)
  )

# glimpse(df_deter)

# Pacote geobr 
munic_geobr <- municipality |> 
  filter(abbrev_state %in% my_states) |> 
  select(name_muni, abbrev_state) |> 
  arrange(name_muni) |>  
  mutate(
    name_muni = stri_trans_general(tolower(name_muni), "Latin-ASCII"),
    #stri_trans_general(..., "Latin-ASCII") → remove acentos e normaliza para comparação;
    #tolower() → padroniza para minúsculas;
    #trimws() → remove espaços extras.
    name_muni = trimws(name_muni)
  )
# munic_geobr <- as.data.frame(munic_geobr)
# glimpse(munic_geobr)

Finalizar tratamento para posterior incorporação

# último tratamento da base deter
# normlizando nomes dos municipios e padronizando nomes de colunas
deter_final <- read_rds("../data/deter_queimadas.rds")
deter_final <- deter_final |> 
  rename(state = UF,
         year = ANO,
         city_ref = MUNICIPALI,
         area_queimada = AREAMUNKM) |> 
  select(state, city_ref, year, area_queimada) |> 
  filter(state %in% my_states,
         year %in% 2015:2023) |> 
  group_by(state, city_ref, year, area_queimada) |> 
  mutate(
    city_ref = stri_trans_general(tolower(city_ref), "Latin-ASCII"),
    city_ref = trimws(city_ref)
  ) |> 
  group_by(year, city_ref, state) |> 
  summarise(area_queimada = mean(area_queimada, na.rm = TRUE))

# write_rds(deter_final, "../data/deter_final.rds")

Gerar mapa - deter

my_year = 2020 # mudar anos (detalhe: alguns anos, como 2016, só tem para 1 estado)

# Fazer filtro por ano antes do join para passar NA = 0 nos municípios faltantes
df_deter_ano <- df_deter |> 
  filter(ANO %in% my_year)

map(2016:2022, ~{
munic_geobr |> 
  filter(abbrev_state %in% my_states) |> 
  left_join(
    df_deter |> 
      filter(ANO == .x) |> 
      group_by(ANO, name_muni) |> 
      summarise(area_queimada_muni_km = sum(area_queimada_muni_km, na.rm = TRUE)),
    by = "name_muni"
  ) |> 
  mutate(area_queimada_muni_km = replace_na(area_queimada_muni_km, 0)) |>
  ggplot() +
  geom_sf(aes(fill=area_queimada_muni_km), color="transparent",
          size=.05, show.legend = TRUE) +
  theme_bw() +
  theme(
    axis.text.x = element_text(size = rel(.9), color = "black"),
    axis.title.x = element_text(size = rel(1.1), color = "black"),
    axis.text.y = element_text(size = rel(.9), color = "black"),
    axis.title.y = element_text(size = rel(1.1), color = "black"),
    legend.text = element_text(size = rel(1), color = "black"),
    legend.title = element_text(face = 'bold', size = rel(1.2))
  ) +
  labs(fill = 'Queimadas (área km)',
       x = 'Longitude',
       y = 'Latitude') +
  scale_fill_viridis_c(
    option = "inferno"
  )})

Gerar mapa - deter (focos de queimadas)

Essa é a outra base do deter, que contém somente os focos de queimadas por estados

Análise pouca representativa (opcionalmente comentar no relatório)

focos_queimadas <- read_rds(("../data/focos_queimadas.rds")
# focos_queimadas <- focos_queimadas |>
#   rename(abbrev_state = uf)


focos_queimadas |> 
  filter(ano %in% 2018:2023) |> 
  group_by(ano, abbrev_state) |> 
  summarise(focos_soma = sum(focos, na.rm = TRUE))
  

map(2018:2023,~{
munic_geobr |>
  left_join(
    focos_queimadas |> 
      filter(ano == .x) |>
      group_by(ano, abbrev_state) |> 
      summarise(focos_soma = sum(focos, na.rm = TRUE)),
    by = "abbrev_state"
    )|> 
      ggplot() +
      geom_sf(aes(fill=focos_soma), color="transparent",
              size=.05, show.legend = TRUE) +
      theme_bw() +
      theme(
        axis.text.x = element_text(size = rel(.9), color = "black"),
        axis.title.x = element_text(size = rel(1.1), color = "black"),
        axis.text.y = element_text(size = rel(.9), color = "black"),
        axis.title.y = element_text(size = rel(1.1), color = "black"),
        legend.text = element_text(size = rel(1), color = "black"),
        legend.title = element_text(face = 'bold', size = rel(1.2))
      ) +
      labs(fill = 'Focos de queimadas',
           x = 'Longitude',
           y = 'Latitude') +
      scale_fill_viridis_c(
        option = "inferno"
      )})

Pré-Processamento

Oliveria, PHM; Panosso, AR

2025-10-02

🛠️ Pré-processameto script

🛠 Preparação dos dados para análise.

Definindo estados

💨 Entrada com a Base: `emissions-sources.rds`

💨 Entrada com a Base: `climate-trace-br.xlsx`

💨 Entrada com a Base: `nasa-xco2.rds`

Filtrando os polígonos do municípios

Gerar mapa

Criando o grid (malha de pontos) para valores não amostrados - nasa-xco2

Interpolação por Krigagem Ordinária

✅ Selecionado o melhor modelo, vamos guardá-lo

Estimativa de XCO2 para o Brasil Central

💨 Entrada com a Base: `gosat-xch4.rds`

Classificando cada ponto em município

Gerar mapa

Criando o grid (malha de pontos) para valores não amostrados - gosat-xch4

Interpolação por Krigagem Ordinária - gosat-xch4

Estimativa de XCH4 para o Brasil Central

🍃 Entrada com a Base: `oco2-sif.rds`

Gerar mapa

Criando o grid (malha de pontos) para valores não amostrados - oco2-sif

Interpolação por Krigagem Ordinária - oco2-sif

Estimativa de SIF para o Brasil Central

🍃 Entrada com a Base: `appeears-modis.rds`

Gerar mapas das variáveis - appeears

💨 Entrada com a Base: `nasa-power.rds`

Classificando cada ponto em município

Gerar mapa

Criando o grid (malha de pontos) para valores não amostrados - nasa-power

Interpolação por Krigagem Ordinária - nasa-power

Estimativa da variável para o Brasil Central - nasa-power

🪓 Entrada com a Base: `prodes-deforestation.rds`

Gerar mapa

🔥 Entrada com a Base: `deter_queimadas.rds`

Padronizando nomes dos municípios da base DETER com os do pacote geobr

Gerar mapa - deter

Gerar mapa - deter (focos de queimadas)

Pré-Processamento

Oliveria, PHM; Panosso, AR

2025-10-02

🛠️ Pré-processameto script

🛠 Preparação dos dados para análise.

Definindo estados

💨 Entrada com a Base: emissions-sources.rds

💨 Entrada com a Base: climate-trace-br.xlsx

💨 Entrada com a Base: nasa-xco2.rds

Filtrando os polígonos do municípios

Gerar mapa

Criando o grid (malha de pontos) para valores não amostrados - nasa-xco2

Interpolação por Krigagem Ordinária

✅ Selecionado o melhor modelo, vamos guardá-lo

Estimativa de XCO2 para o Brasil Central

💨 Entrada com a Base: gosat-xch4.rds

Classificando cada ponto em município

Gerar mapa

Criando o grid (malha de pontos) para valores não amostrados - gosat-xch4

Interpolação por Krigagem Ordinária - gosat-xch4

Estimativa de XCH4 para o Brasil Central

🍃 Entrada com a Base: oco2-sif.rds

Gerar mapa

Criando o grid (malha de pontos) para valores não amostrados - oco2-sif

Interpolação por Krigagem Ordinária - oco2-sif

Estimativa de SIF para o Brasil Central

🍃 Entrada com a Base: appeears-modis.rds

Gerar mapas das variáveis - appeears

💨 Entrada com a Base: nasa-power.rds

Classificando cada ponto em município

Gerar mapa

Criando o grid (malha de pontos) para valores não amostrados - nasa-power

Interpolação por Krigagem Ordinária - nasa-power

Estimativa da variável para o Brasil Central - nasa-power

🪓 Entrada com a Base: prodes-deforestation.rds

Gerar mapa

🔥 Entrada com a Base: deter_queimadas.rds

Padronizando nomes dos municípios da base DETER com os do pacote geobr

Gerar mapa - deter

Gerar mapa - deter (focos de queimadas)

💨 Entrada com a Base: `emissions-sources.rds`

💨 Entrada com a Base: `climate-trace-br.xlsx`

💨 Entrada com a Base: `nasa-xco2.rds`

💨 Entrada com a Base: `gosat-xch4.rds`

🍃 Entrada com a Base: `oco2-sif.rds`

🍃 Entrada com a Base: `appeears-modis.rds`

💨 Entrada com a Base: `nasa-power.rds`

🪓 Entrada com a Base: `prodes-deforestation.rds`

🔥 Entrada com a Base: `deter_queimadas.rds`