Skip to content

Latest commit

 

History

History
833 lines (723 loc) · 34.3 KB

README.md

File metadata and controls

833 lines (723 loc) · 34.3 KB
Raw

Determinação de fontes e sumidouros de carbono atmosférico no Brasil por meio da análise de regressão linear

Aluno: Alan R. Panosso

Data: 31/08/2021

Projeto final apresentado os instrutores Athos Damiani e Fernando Correa da curso-R como parte das exigências para a finalização do curso de Regressão Linear (Trilha de Machine Learning - Turma Agosto 2021).

Introdução

As concentrações atmosféricas de gases do efeito estufa (GEE) têm aumentado a níveis preocupantes. De acordo com a Organização Meteorológica Mundial (WMO) as concentrações atmosféricas de dióxido de carbono (CO2), metano (CH4) e óxido nitroso (N2O) atingiram novas máximas no ano de 2015 com CO2 a 400,0 ± 0,1 ppm, CH4 em 1845 ± 2 ppb e N2O em 328,0 ± 0,1 ppb, valores que representam, respectivamente, 144%, 256% e 121% dos níveis pré-industriais (WMO, 2016).

O Observatório Orbital de Carbono-2 (OCO-2) foi projetado pela Agência Espacial Americana (National Aeronautics and Space Administration - NASA) para apoiar a quantificação e o monitoramento das emissões antropogênicas de CO2.

O satélite OCO-2 foi lançado em julho de 2014 e desde então mede a concentração de CO2 atmosférico indiretamente por meio da intensidade da radiação solar refletida em função da presença de dióxido de carbono em uma coluna de ar. Essas medições resultam em mapas espaciais densos e em escala fina de frações molares médias de coluna de ar seco de dióxido de carbono (XCO2).

Nesse contexto, a variação de XCO2 opode ser modelada por meio da análise de regressão linear simples, uma vez que se espera um aumento dessas concentrações com o passar do tempo. Em adição, as estimativas do coeficiente angular β1 fornece informações importantes para uma determinada região, haja visto que se β1’ (observado para essa região) for significativamente maior ao β1 padrão (observado para uma macro-região), tal área poderá ser considerada uma potencial fonte de carbono para a atmosfera, caso contrário (β1’ < β1), a área em questão poderá ser considerada um sumidouro de CO2 atmosférico, mitigando o efeito estufa adicional e, consequentemente, as mudanças climáticas globais.

Hipótese e objetivo

A hipótese do projeto é que essa tendência de aumento da concentração atmosférica de CO2 pode ser utilizada como um indicativo para a classificação de áreas como fontes ou sumidouros de CO2 utilizando as estimativas de β1 provenientes da análise de regressão linear simples.

Material e Métodos

Para a aquisição de dados será utilizado metodologia desenvolvida e apresentada anteriormente no curso de R para Ciências de dados 2, ministrados pela Curso-r (Projeto Final r4ds2), onde foram utilizadas técnicas de web scraping e faxina de dados para obtenção dos valores de XCO2.

Breve descrição das variáveis da base:

longitude: coordenada geográfica que especifica a posição leste-oeste de um ponto na superfície da Terra;

longitude_bnds: são, respectivamente, os limites superior e inferior da coordenada, onde a longitude para um ponto foi dada pela média desses limites;

latitude: é uma coordenada geográfica que especifica a posição norte-sul de um ponto na superfície da Terra;

latitude_bnds: são, respectivamente, os limites superior e inferior da coordenada, onde a latitude para um ponto foi dada pela média desses limites;

time_yyyymmddhhmmss: data de leitura, em ano, mês, dia, horas minutos e segundos;

time_bnds_yyyymmddhhmmss: limites de tempo utilizados para o cálculo da data de leitura;

altitude_km: altitude média em km;

alt_bnds_km: limites da altitude, 0 (nível do mar) e altitude do satélite no momento de leitura;

fluorescence_offset_relative_771nm_idp: Fração de radiância de nível contínuo explicada por um termo de deslocamento aditivo na janela espectral de 757 nm (sem unidade);

fluorescence_offset_relative_757nm_idp: Fração da radiância de nível contínuo explicada por um termo de deslocamento aditivo na janela espectral de 771 nm (sem unidade);

xco2_moles_mole_1: Fração molar de ar seco de CO2 em média da coluna.

oco2 <- readr::read_rds("data/oco2.rds")
dplyr::glimpse(oco2)
#> Rows: 361,615
#> Columns: 11
#> $ longitude                              <dbl> -74.58225, -74.58225, -74.58225…
#> $ longitude_bnds                         <chr> "-74.70703125:-74.4574652778", …
#> $ latitude                               <dbl> -30.22572489, -29.97654828, -29…
#> $ latitude_bnds                          <chr> "-30.3503131952:-30.1011365845"…
#> $ time_yyyymmddhhmmss                    <dbl> 2.014091e+13, 2.014091e+13, 2.0…
#> $ time_bnds_yyyymmddhhmmss               <chr> "20140909000000:20140910000000"…
#> $ altitude_km                            <dbl> 3307.8, 3307.8, 3307.8, 3307.8,…
#> $ alt_bnds_km                            <chr> "0.0:6615.59960938", "0.0:6615.…
#> $ fluorescence_offset_relative_771nm_idp <dbl> 0.012406800, 0.010696600, -0.00…
#> $ fluorescence_offset_relative_757nm_idp <dbl> -3.58630e+00, 8.81219e-02, -3.6…
#> $ xco2_moles_mole_1                      <dbl> 0.000394333, 0.000395080, 0.000…

Será necessário transformar os dados de XCO2, para ppm em seguida criar as variáveis de data a partir da variável time_yyyymmddhhmmss.

oco2<-oco2 |>
         dplyr::mutate(
           xco2 = xco2_moles_mole_1*1e06,
           data = lubridate::ymd_hms(time_yyyymmddhhmmss),
           year = lubridate::year(data),
           month = lubridate::month(data),
           day = lubridate::day(data),
           dia = difftime(data,"2014-01-09", units = "days"),
           day_week = lubridate::wday(data),
           month_year = lubridate::make_date(year, month, 1) )

Existe uma tendência de aumento monotônica mundial da concentração de CO2 na atmosfera, assim, ela deve ser modelada para obtermos β1 para ser considerado o padrão para comparação às tendências de uma determinada região. Devido à periodicidade de retorno do satélite em um ponto (ao redor de 16 dias) os dados devem ser agrupados pelo mês dentro de um determinado ano.

oco2 |> 
  dplyr::filter(year %in% 2015:2020) |> 
  dplyr::group_by(dia) |> 
  dplyr::summarise(xco2_mean = mean(xco2, na.rm =TRUE)) |> 
  ggplot2::ggplot(ggplot2::aes(x=dia,y=xco2_mean )) +
  ggplot2::geom_point(shape=21,color="black",fill="gray") +
  ggplot2::geom_line(color="red") +
  ggplot2::geom_smooth(method = "lm") +
  ggpubr::stat_regline_equation(ggplot2::aes(
  label =  paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~~"))) +
  ggplot2::theme_bw()

oco2 |>
  dplyr::filter(year %in% 2015:2020) |>
  dplyr::group_by(year, dia) |>
  dplyr::summarise(xco2_mean = mean(xco2, na.rm =TRUE)) |>
  ggplot2::ggplot(ggplot2::aes(x=dia,y=xco2_mean,
                               fill=forcats::as_factor(year))) +
  ggplot2::geom_point(shape=21,color="black") +
  #ggplot2::geom_line(color="red") +
  ggplot2::geom_smooth(method = "lm") +
  ggplot2::facet_wrap(~year,scales = "free")+
  ggpubr::stat_regline_equation(ggplot2::aes(
  label =  paste(..eq.label.., ..rr.label.., sep = "*plain(\",\")~~"))) +
  ggplot2::theme_bw()
#> `summarise()` has grouped output by 'year'. You can override using the
#> `.groups` argument.
#> Don't know how to automatically pick scale for object of type <difftime>.
#> Defaulting to continuous.
#> `geom_smooth()` using formula = 'y ~ x'

Para ajustar o modelo linear, usamos lm().

oco2_aux <- oco2 |> 
  dplyr::filter(year == 2015) |>
  dplyr::group_by(dia) |> 
  dplyr::summarise(xco2_mean = mean(xco2, na.rm =TRUE))
mod <- lm(xco2_mean ~ dia, data = oco2_aux)
br <- mod$coefficients[2] # beta regional
ep <- summary(mod)$coefficients[2,2] # erro padrão

# limite inferior do Beta regional
limite_inferior_beta_regional <- br - ep

# limite superior do Beta regional
limite_superior_beta_regional <- br + ep

Vamos olhar o diagnóstico da análise.

broom::augment(mod, interval="confidence")
#> # A tibble: 317 × 10
#>    xco2_mean dia   .fitted .lower .upper .resid   .hat .sigma .cooksd .std.resid
#>        <dbl> <drt>   <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>   <dbl>      <dbl>
#>  1      391. 357.…    392.   392.   392. -0.812 0.0134   1.22 3.06e-3     -0.672
#>  2      393. 358.…    392.   392.   392.  1.29  0.0133   1.22 7.60e-3      1.06 
#>  3      393. 359.…    392.   392.   392.  1.00  0.0131   1.22 4.58e-3      0.829
#>  4      393. 360.…    392.   392.   392.  0.854 0.0130   1.22 3.29e-3      0.706
#>  5      392. 361.…    392.   392.   392.  0.204 0.0129   1.22 1.86e-4      0.168
#>  6      392. 366.…    392.   392.   392. -0.389 0.0124   1.22 6.50e-4     -0.322
#>  7      393. 367.…    392.   392.   392.  0.573 0.0123   1.22 1.40e-3      0.473
#>  8      393. 368.…    392.   392.   392.  0.879 0.0122   1.22 3.26e-3      0.726
#>  9      392. 369.…    392.   392.   392.  0.352 0.0121   1.22 5.19e-4      0.291
#> 10      393. 370.…    392.   392.   392.  1.19  0.0120   1.22 5.89e-3      0.985
#> # ℹ 307 more rows
plot(mod)

cooks.distance(mod)
#>            1            2            3            4            5            6 
#> 3.055247e-03 7.597057e-03 4.580254e-03 3.292568e-03 1.859488e-04 6.499081e-04 
#>            7            8            9           10           11           12 
#> 1.395253e-03 3.256227e-03 5.185250e-04 5.886842e-03 3.603522e-06 5.516671e-04 
#>           13           14           15           16           17           18 
#> 8.850543e-04 2.806128e-03 1.350571e-03 1.352814e-06 6.925387e-04 3.456922e-03 
#>           19           20           21           22           23           24 
#> 2.092351e-03 7.621595e-04 1.395225e-02 2.648997e-03 2.675939e-04 2.984946e-06 
#>           25           26           27           28           29           30 
#> 1.563528e-03 4.901811e-04 1.995668e-03 9.084349e-03 5.843687e-03 2.077896e-03 
#>           31           32           33           34           35           36 
#> 1.760198e-03 3.571221e-03 9.101033e-03 3.812346e-04 8.390584e-03 2.005749e-04 
#>           37           38           39           40           41           42 
#> 3.926647e-03 5.876160e-03 1.352124e-02 1.498654e-02 7.068567e-03 1.146305e-03 
#>           43           44           45           46           47           48 
#> 1.888846e-02 2.555641e-04 5.158430e-03 1.940192e-04 2.861215e-02 3.983845e-02 
#>           49           50           51           52           53           54 
#> 6.439778e-04 7.267364e-03 1.766793e-03 2.623932e-05 5.367644e-06 4.639781e-04 
#>           55           56           57           58           59           60 
#> 1.132913e-03 1.251678e-02 1.424559e-02 3.885648e-06 2.679798e-04 6.501314e-04 
#>           61           62           63           64           65           66 
#> 2.208402e-03 2.121131e-04 3.573741e-03 8.706761e-04 9.375163e-04 4.898626e-02 
#>           67           68           69           70           71           72 
#> 1.222203e-03 2.448546e-02 7.243650e-07 5.476660e-02 8.777554e-03 1.047693e-03 
#>           73           74           75           76           77           78 
#> 9.163364e-04 1.213957e-02 1.756302e-02 7.231332e-04 4.357834e-04 6.094724e-06 
#>           79           80           81           82           83           84 
#> 1.932321e-02 7.714448e-04 6.727077e-05 4.733975e-03 8.841795e-07 6.233330e-04 
#>           85           86           87           88           89           90 
#> 4.719103e-05 8.251322e-04 4.749394e-03 4.614332e-03 6.587137e-04 8.577621e-06 
#>           91           92           93           94           95           96 
#> 9.731794e-04 7.944171e-05 8.122349e-05 2.038396e-03 2.016680e-03 2.129724e-04 
#>           97           98           99          100          101          102 
#> 1.119951e-06 1.483962e-04 7.735784e-05 1.010943e-04 2.697623e-03 8.782719e-05 
#>          103          104          105          106          107          108 
#> 1.094433e-04 1.803959e-04 1.683447e-03 1.075276e-04 3.332941e-03 2.553531e-05 
#>          109          110          111          112          113          114 
#> 1.666449e-05 1.734787e-03 2.074535e-03 5.185004e-03 1.804371e-03 5.281733e-03 
#>          115          116          117          118          119          120 
#> 1.827251e-03 2.530600e-03 3.833097e-04 8.532830e-04 8.289624e-04 1.800869e-03 
#>          121          122          123          124          125          126 
#> 8.261726e-04 5.410450e-04 3.386293e-03 7.617568e-03 1.339824e-03 1.323014e-03 
#>          127          128          129          130          131          132 
#> 3.150965e-04 6.371985e-03 9.730021e-04 1.999290e-03 1.383460e-03 2.387345e-03 
#>          133          134          135          136          137          138 
#> 1.676880e-03 6.818059e-04 2.915351e-03 4.782992e-04 3.015790e-03 2.435112e-04 
#>          139          140          141          142          143          144 
#> 1.173490e-03 7.952517e-04 2.726686e-03 4.292316e-03 2.510919e-03 3.847819e-03 
#>          145          146          147          148          149          150 
#> 3.328335e-03 3.942858e-03 2.327937e-03 3.486587e-03 3.052312e-03 3.846841e-03 
#>          151          152          153          154          155          156 
#> 6.136751e-04 2.179549e-03 1.545625e-03 2.048557e-06 3.243162e-03 7.573242e-05 
#>          157          158          159          160          161          162 
#> 2.610463e-03 3.099234e-03 1.494351e-03 3.275413e-03 1.392188e-04 3.305571e-03 
#>          163          164          165          166          167          168 
#> 2.411267e-03 5.202507e-04 4.023833e-03 2.633894e-07 3.176671e-03 3.406761e-04 
#>          169          170          171          172          173          174 
#> 9.798096e-04 7.062379e-04 9.203146e-06 1.689173e-03 2.030393e-04 3.123383e-03 
#>          175          176          177          178          179          180 
#> 3.266551e-04 1.478282e-03 6.243266e-04 8.765316e-04 1.588683e-04 2.017707e-03 
#>          181          182          183          184          185          186 
#> 1.545060e-03 4.403041e-05 4.057039e-03 6.459277e-04 4.231227e-04 7.805765e-06 
#>          187          188          189          190          191          192 
#> 1.072661e-03 8.872802e-04 1.653733e-04 3.329200e-03 6.823652e-04 1.051280e-03 
#>          193          194          195          196          197          198 
#> 1.980442e-05 7.718141e-06 4.839467e-04 1.756631e-04 9.177032e-05 2.073792e-05 
#>          199          200          201          202          203          204 
#> 7.124094e-03 4.066240e-04 6.172188e-06 1.078655e-04 1.835665e-03 4.935996e-06 
#>          205          206          207          208          209          210 
#> 1.280712e-05 2.268246e-05 1.293692e-03 1.504142e-04 2.547836e-05 2.620210e-05 
#>          211          212          213          214          215          216 
#> 1.095968e-07 1.573884e-04 2.438335e-04 2.471376e-06 1.291552e-04 8.870286e-04 
#>          217          218          219          220          221          222 
#> 1.961707e-03 2.268347e-04 1.732516e-06 8.410656e-05 2.009838e-04 4.866213e-04 
#>          223          224          225          226          227          228 
#> 4.516526e-05 8.259732e-04 1.410293e-03 6.136655e-03 3.234320e-04 7.513052e-05 
#>          229          230          231          232          233          234 
#> 4.082083e-04 1.294844e-04 3.179133e-03 2.507150e-03 1.055601e-03 2.887959e-05 
#>          235          236          237          238          239          240 
#> 1.996825e-03 7.461275e-05 1.461838e-05 2.616779e-04 3.717024e-03 2.133591e-03 
#>          241          242          243          244          245          246 
#> 7.753288e-03 2.918223e-03 8.795129e-07 4.355446e-04 4.210778e-04 1.175078e-04 
#>          247          248          249          250          251          252 
#> 5.395151e-04 4.123921e-05 5.564857e-03 9.889962e-05 3.070257e-05 2.451659e-05 
#>          253          254          255          256          257          258 
#> 8.146741e-04 1.197111e-04 7.264636e-06 2.163958e-03 6.435351e-05 9.541303e-03 
#>          259          260          261          262          263          264 
#> 1.599980e-02 2.724419e-03 1.854510e-03 7.370544e-03 1.287532e-03 3.584743e-03 
#>          265          266          267          268          269          270 
#> 2.191274e-03 6.155308e-04 4.255415e-03 6.615302e-04 3.611688e-04 3.615010e-04 
#>          271          272          273          274          275          276 
#> 4.825631e-03 2.217225e-02 3.402189e-06 4.926265e-06 1.812464e-05 1.914446e-04 
#>          277          278          279          280          281          282 
#> 2.051160e-03 6.859079e-04 3.809780e-03 1.522015e-03 1.664252e-02 7.383826e-03 
#>          283          284          285          286          287          288 
#> 2.327033e-03 9.724409e-07 1.515583e-04 2.439975e-03 3.263679e-03 5.762377e-04 
#>          289          290          291          292          293          294 
#> 1.220151e-04 5.633927e-04 7.224024e-05 7.610696e-04 2.161567e-05 7.820580e-03 
#>          295          296          297          298          299          300 
#> 6.254593e-04 7.875021e-02 1.184216e-03 8.799616e-04 1.085939e-03 2.896817e-04 
#>          301          302          303          304          305          306 
#> 5.132459e-03 2.182456e-03 1.423779e-02 9.297009e-04 1.598796e-04 3.991303e-06 
#>          307          308          309          310          311          312 
#> 2.487110e-03 3.138805e-03 3.225408e-03 3.021315e-04 2.630219e-03 1.591272e-02 
#>          313          314          315          316          317 
#> 1.351446e-04 5.807543e-03 1.264129e-02 7.193662e-05 3.172273e-03

A próxima operação é selecionar na base de dados somente os pontos pertencentes ao território brasileiro. Assim vamos utilizar o pacote {geobr} para criarmos os filtros a partir dos polígonos das diferentes regiões do Brasil.

regiao <- geobr::read_region(showProgress = FALSE)
#> Loading required namespace: sf
#> Using year/date 2010
br <- geobr::read_country(showProgress = FALSE)
#> Using year/date 2010

Agora podemos extrair os polígonos.

### Polígono Brasil
pol_br <- br$geom |> purrr::pluck(1) |> as.matrix()
### Polígonos das Regiões
pol_norte <- regiao$geom |> purrr::pluck(1) |> as.matrix()
pol_nordeste <- regiao$geom |> purrr::pluck(2) |> as.matrix()
pol_sudeste <- regiao$geom |> purrr::pluck(3) |> as.matrix()
pol_sul <- regiao$geom |> purrr::pluck(4) |> as.matrix()
pol_centroeste<- regiao$geom |> purrr::pluck(5) |> as.matrix()
# Retirando alguns pontos
pol_br <- pol_br[pol_br[,1]<=-34,]
pol_br <- pol_br[!((pol_br[,1]>=-38.8 & pol_br[,1]<=-38.6) &
                              (pol_br[,2]>= -19 & pol_br[,2]<= -16)),]

# Arrumando alguns pontos
pol_nordeste <- pol_nordeste[pol_nordeste[,1]<=-34,]
pol_nordeste <- pol_nordeste[!((pol_nordeste[,1]>=-38.7 & pol_nordeste[,1]<=-38.6) & pol_nordeste[,2]<= -15),]

# retirando pontos do sudeste
pol_sudeste <- pol_sudeste[pol_sudeste[,1]<=-30,]

Plot de todos os pontos.

 br |>
   ggplot2::ggplot() +
   ggplot2::geom_sf(fill="#2D3E50", color="#FEBF57",
           size=.15, show.legend = FALSE) +
   ggplot2::geom_point(data=oco2 |> dplyr::sample_n(20000) |> 
                         dplyr::filter(year == 2014) ,
              ggplot2::aes(x=longitude, y=latitude),
              shape=17,
              col="red",
              alpha=0.2)

Definindo uma função para criar as flags das diferentes regiões do país, a partir da função point.in.polygon do pacote {sp}.

def_pol <- function(x, y, pol){
  as.logical(sp::point.in.polygon(point.x = x,
                                  point.y = y,
                                  pol.x = pol[,1],
                                  pol.y = pol[,2]))
}

Vamos criar o filtro para os pontos pertencentes ao polígono do Brasill. Devemos salientar que esses dados estão com a tendência de aumento ao longo do tempo, pois, vamos salvar o arquivo para posterior disponibilização.

oco2 <- oco2 |>
          dplyr::mutate(
            flag_br = def_pol(longitude, latitude, pol_br),
            flag_norte = def_pol(longitude, latitude, pol_norte),
            flag_nordeste = def_pol(longitude, latitude, pol_nordeste),
            flag_sul = def_pol(longitude, latitude, pol_sul),
            flag_sudeste = def_pol(longitude, latitude, pol_sudeste),
            flag_centroeste = def_pol(longitude, latitude, pol_centroeste)
          )

Verificação dos pontos dentro do território brasileiro.

 br |>
   ggplot2::ggplot() +
   ggplot2::geom_sf(fill="#2D3E50", color="#FEBF57",
           size=.15, show.legend = FALSE) +
   ggplot2::geom_point(data=oco2 |> 
                         dplyr::filter(flag_norte |
                                         flag_sul |
                                         flag_sudeste|
                                         flag_centroeste|
                                         flag_nordeste) |> 
                         dplyr::sample_n(20000) |> 
                         dplyr::filter(year == 2014) ,
              ggplot2::aes(x=longitude, y=latitude),
              shape=17,
              col="red",
              alpha=0.2)

Pra garantir a reprodutibilidade desse material, vamos salvar/disponibilizar os dados na base oco2_br_trend.rds, somente com os pontos dentro do território nacional.

  readr::write_rds(oco2 |>
                     dplyr::filter(flag_norte |
                                         flag_sul |
                                         flag_sudeste|
                                         flag_centroeste|
                                         flag_nordeste), 
                     "data/oco2_br_trend.rds")

Vamos ler o banco de dados [com a tendência].

oco2_br_trend <- readr::read_rds("data/oco2_br_trend.rds")

Devemos criar a análise para um ano específico, pois a comparação do artigo será ano a ano.

oco2_nest <- oco2_br_trend |>
  dplyr::filter(year == 2015) |> 
  tibble::as_tibble() |> 
  dplyr::mutate(quarter = lubridate::quarter(data),
                quarter_year = lubridate::make_date(year, quarter, 1)) |>   tidyr::pivot_longer(
    starts_with("flag"),
    names_to = "region",
    values_to = "flag",
  ) |> 
  dplyr::filter(flag) |> 
  dplyr::mutate(region = stringr::str_remove(region,"flag_")) |> 
  dplyr::group_by(region, longitude, latitude, dia) |> 
  dplyr::summarise(xco2_mean = mean(xco2, na.rm=TRUE)) |> 
  dplyr::mutate(
    regi = region,
    id_time = dia
  ) |> 
  dplyr::group_by(region, latitude, longitude) |> 
  tidyr::nest()

Função para construção da análise de regressão linear para cada pixel, e diagnósticos.

linear_reg <- function(df, output="beta1"){
  # Modelo para cada pixel
  modelo <- lm(xco2_mean ~ dia, data=df)
  beta_1 <- c(summary(modelo)$coefficients[2])
  
  # Definindo o modelo
  if(output=="beta1"){ 
    return(beta_1)
  }
  
  # Salvando o valor P
  if(output=="p_value"){
    if(is.nan(beta_1)){
      beta_1 <- 0
      p <- 1
    }else{
      p <- summary(modelo)$coefficients[2,4]
      if(is.nan(p)) p <- 1
    }
    return(p)
  }
  
  # Criando gráfico
  if(output=="plot"){
    plot <- df |> 
      ggplot2::ggplot(ggplot2::aes(x=dia,y=xco2_mean)) +
      ggplot2::geom_point() +
      ggplot2::theme_bw()
    return(plot)
  }
  if(output=="hist"){
    hist <- df |> 
      ggplot2::ggplot(ggplot2::aes(x=xco2_mean, y=..density..)) +
      ggplot2::geom_histogram(bins=10, color="black", fill="lightgray") +
      ggplot2::geom_density()+
      ggplot2::theme_bw()
    return(hist)
  }
  
  # Anomalia é o Xco2 do regional menos o Xco2 do pixel, melhor é o contrário.
  if(output == "partial"){
     partial <- df |> 
       dplyr::summarise(xco2 = mean(xco2_mean), na.mr=TRUE) |> 
       dplyr::pull(xco2)
     return(partial)
  }
  
    if(output == "n"){
     return(nrow(df))
   }
}

Vamos aplicar a função para cada ponto de amostragem do satélite (pixel).

oco2_nest <- oco2_nest |> 
  dplyr::mutate( 
    beta_line = purrr::map(data,linear_reg, output="beta1"),
    p_value = purrr::map(data,linear_reg, output="p_value"),
    partial = purrr::map(data,linear_reg, output="partial"),
    n_obs = purrr::map(data,linear_reg, output="n")
    #plot = purrr::map(data,linear_reg, output="plot"),
    #hist = purrr::map(data,linear_reg, output="hist")
  )
oco2_nest$data[[1]]
#> # A tibble: 6 × 4
#>   dia           xco2_mean regi       id_time      
#>   <drtn>            <dbl> <chr>      <drtn>       
#> 1 454.4167 days      394. centroeste 454.4167 days
#> 2 525.4167 days      392. centroeste 525.4167 days
#> 3 557.4167 days      394. centroeste 557.4167 days
#> 4 589.4167 days      397. centroeste 589.4167 days
#> 5 598.4167 days      398. centroeste 598.4167 days
#> 6 621.4167 days      394. centroeste 621.4167 days
# oco2_nest |> 
#   dplyr::filter(region == "norte") |> 
#   dplyr::filter(p_value < 0.05, beta_line < 0) |> 
#   dplyr::pull(plot)
oco2_nest |> 
  # dplyr::filter(p_value < 0.05) |>
  dplyr::filter(n_obs > 5) |> 
  # dplyr::mutate(class = ifelse(beta_line > limite_inferior_beta_regional,
  #                              1,ifelse(beta_line < limite_inferior_beta_regional, -1, 0))
  #               ) |> 
  dplyr::select(longitude, latitude, n_obs) |> 
  ggplot2::ggplot(ggplot2::aes(x=longitude, y=latitude, color = n_obs)) +
  ggplot2::geom_point()

oco2_aux <- oco2_nest |> 
  # dplyr::filter(region == "norte") |> 
  # dplyr::filter(p_value < 0.05) |> 
  dplyr::filter(n_obs > 7) |> 
  tidyr::unnest(cols = c(beta_line, partial)) |>
  dplyr::ungroup() |>
  dplyr::select(longitude, latitude, beta_line, partial)

q3_oco2 <- oco2_aux |> dplyr::pull(beta_line) |> quantile(.75)

oco2_aux <- oco2_aux |> 
  dplyr::mutate(
    anomaly =  partial - oco2_aux |> 
      dplyr::pull(partial) |>  
      mean(),
    Dbeta = beta_line - oco2_aux |> 
      dplyr::pull(beta_line) |> mean(na.rm=TRUE)
  )
q3_anom <- oco2_aux |> dplyr::pull(anomaly) |> quantile(.75)


oco2_aux <- oco2_aux |> 
  dplyr::mutate(
    beta_index =  ifelse(beta_line <=q3_oco2, 0, 1)
  )

# Mapeando
oco2_aux |> 
  ggplot2::ggplot(ggplot2::aes(x=longitude, y=latitude) ) + 
  ggplot2::geom_point()

oco2_aux |> 
  ggplot2::ggplot(ggplot2::aes(x=beta_line)) +
  ggplot2::geom_histogram(bins=30,
                          fill="orange",
                          color="black") +
  ggplot2::labs(x="βpixel",y="Count") +
  ggplot2::geom_vline(xintercept = q3_oco2,
                      color = "red",
                      lty=2) +
  gghighlight::gghighlight(beta_line > q3_oco2,
                           unhighlighted_params = list(
                               color = "darkgray",
                               fill = "lightgray")) +
  ggplot2::theme_minimal()
#> Warning: Could not calculate the predicate for layer 2; ignored

oco2_aux |>
  ggplot2::ggplot(ggplot2::aes(x=anomaly)) +
  ggplot2::geom_histogram(bins=30,
                          fill="lightblue",
                          color="black") +
  ggplot2::labs(x="Anomaly",y="Count") +
  ggplot2::geom_vline(xintercept = q3_anom,
                      color = "red",
                      lty=2) +
  gghighlight::gghighlight(anomaly > q3_anom,
                           unhighlighted_params = list(
                               color = "darkgray",
                               fill = "lightgray")) +
  ggplot2::theme_minimal()
#> Warning: Could not calculate the predicate for layer 2; ignored

sp::coordinates(oco2_aux)=~ longitude+latitude  
form_beta<-beta_line~1
form_anom<-anomaly~1
form_index<-beta_index~1

Semivariograma para beta

vari_beta <- gstat::variogram(form_beta, data=oco2_aux)
m_beta <- gstat::fit.variogram(vari_beta,fit.method = 7,
                               gstat::vgm(1, "Sph", 8, 1))
plot(vari_beta,model=m_beta, col=1,pl=F,pch=16)

Semivariograma para anomalia

vari_anom<-gstat::variogram(form_anom, data=oco2_aux)
m_anom <- gstat::fit.variogram(vari_anom,gstat::vgm(.8,"Sph",9,.2))
plot(vari_anom, model=m_anom, col=1,pl=F,pch=16)

Semivariograma para beta_index

vari_index <- gstat::variogram(form_index, data=oco2_aux,
                               cutoff = 5, 
                               width = 5/15)
m_index <- gstat::fit.variogram(vari_index,fit.method = 7,
                               gstat::vgm(1, "Sph", 4, 1))
plot(vari_index,model=m_index, col=1,pl=F,pch=16)

x<-oco2_aux$longitude
y<-oco2_aux$latitude
dis <- .1 #Distância entre pontos
grid <- expand.grid(X=seq(min(x),max(x),dis), Y=seq(min(y),max(y),dis))
sp::gridded(grid) = ~ X + Y
ko_beta<-gstat::krige(formula=form_beta, oco2_aux, grid, model=m_beta, 
    block=c(1,1),
    nsim=0,
    na.action=na.pass,
    debug.level=-1,  
    )
#> [using ordinary kriging]
#>   1% done  3% done  5% done  6% done  7% done  8% done  9% done 10% done 11% done 12% done 13% done 14% done 15% done 16% done 17% done 18% done 19% done 20% done 21% done 22% done 23% done 24% done 25% done 26% done 27% done 28% done 29% done 30% done 31% done 32% done 33% done 34% done 35% done 36% done 37% done 38% done 39% done 40% done 41% done 42% done 43% done 44% done 45% done 46% done 47% done 48% done 49% done 50% done 51% done 52% done 53% done 55% done 57% done 59% done 60% done 62% done 64% done 66% done 68% done 69% done 70% done 71% done 72% done 73% done 74% done 75% done 76% done 77% done 78% done 79% done 80% done 81% done 82% done 83% done 84% done 85% done 86% done 87% done 88% done 89% done 90% done 91% done 92% done 93% done 94% done 95% done 96% done 97% done 98% done 99% done100% done

ko_anom<-gstat::krige(formula=form_anom, oco2_aux, grid, model=m_anom, 
    block=c(0,0),
    nsim=0,
    na.action=na.pass,
    debug.level=-1,  
    )
#> [using ordinary kriging]
#>   0% done  1% done  2% done  3% done  4% done  5% done  6% done  7% done  8% done  9% done 10% done 11% done 12% done 13% done 14% done 15% done 16% done 17% done 18% done 19% done 20% done 21% done 22% done 23% done 24% done 25% done 26% done 27% done 28% done 29% done 30% done 31% done 32% done 33% done 34% done 35% done 36% done 37% done 38% done 39% done 40% done 41% done 42% done 44% done 46% done 48% done 51% done 54% done 56% done 58% done 61% done 63% done 66% done 68% done 70% done 73% done 75% done 78% done 80% done 83% done 85% done 88% done 90% done 93% done 95% done 98% done100% done

ko_index<-gstat::krige(formula=form_index, oco2_aux, grid, model=m_index, 
    block=c(0,0),
    nsim=0,
    na.action=na.pass,
    debug.level=-1,  
    )
#> [using ordinary kriging]
#>   2% done  4% done  7% done  9% done 11% done 14% done 16% done 19% done 22% done 24% done 27% done 29% done 32% done 34% done 37% done 39% done 42% done 45% done 47% done 50% done 52% done 55% done 58% done 60% done 63% done 65% done 68% done 70% done 73% done 75% done 78% done 81% done 83% done 86% done 88% done 91% done 93% done 96% done 98% done100% done
mapa <- geobr::read_state(showProgress = FALSE)
#> Using year/date 2010
get_pol_in_pol <- function(indice, lista, gradeado){
  poligono <- lista |> purrr::pluck(indice) |> as.matrix()
  flag <- def_pol(gradeado$X, gradeado$Y, poligono)
  return(flag)
}
flag <- purrr::map_dfc(1:27, get_pol_in_pol, lista=mapa$geom, gradeado = grid)
#> New names:
#> • `` -> `...1`
#> • `` -> `...2`
#> • `` -> `...3`
#> • `` -> `...4`
#> • `` -> `...5`
#> • `` -> `...6`
#> • `` -> `...7`
#> • `` -> `...8`
#> • `` -> `...9`
#> • `` -> `...10`
#> • `` -> `...11`
#> • `` -> `...12`
#> • `` -> `...13`
#> • `` -> `...14`
#> • `` -> `...15`
#> • `` -> `...16`
#> • `` -> `...17`
#> • `` -> `...18`
#> • `` -> `...19`
#> • `` -> `...20`
#> • `` -> `...21`
#> • `` -> `...22`
#> • `` -> `...23`
#> • `` -> `...24`
#> • `` -> `...25`
#> • `` -> `...26`
#> • `` -> `...27`
flag_br <- apply(flag, 1, sum) != 0
tibble::as_tibble(ko_beta) |> 
  tibble::add_column(flag_br) |>
  dplyr::filter(flag_br) |>
  ggplot2::ggplot(ggplot2::aes(x=X, y=Y),color="black") + 
  ggplot2::geom_tile(ggplot2::aes(fill = var1.pred)) +
  ggplot2::scale_fill_gradient(low = "yellow", high = "blue") + 
  ggplot2::coord_equal()+
  ggplot2::labs(fill="βpixel") +
  ggplot2::theme_bw()

tibble::as_tibble(ko_anom) |> 
  tibble::add_column(flag_br) |> 
  dplyr::filter(flag_br) |> 
  ggplot2::ggplot(ggplot2::aes(x=X, y=Y),color="black") + 
  ggplot2::geom_tile(ggplot2::aes(fill = var1.pred)) +
  ggplot2::scale_fill_gradient(low = "yellow", high = "blue") + 
  ggplot2::coord_equal()+
  ggplot2::labs(fill="Anomaly") +
  ggplot2::theme_bw()

tibble::as_tibble(ko_index) |> 
  tibble::add_column(flag_br) |> 
  dplyr::filter(flag_br) |> 
  ggplot2::ggplot(ggplot2::aes(x=X, y=Y),color="black") + 
  ggplot2::geom_tile(ggplot2::aes(fill = var1.pred)) +
  ggplot2::scale_fill_gradient(low = "yellow", high = "blue") + 
  ggplot2::coord_equal()+
  ggplot2::labs(fill="P(x = source)") +
  ggplot2::theme_bw()