Como fazer group_by() + mean() de datasets do tipo "sf"?

itamar_cordeiro · Agosto 12, 2022, 12:59am

Estou trabalhando com o seguinte dataset

As latitudes e longitudes representam estações que coletam dados. Neste caso, estamos trabalhando apenas com a Variável “v1”.

Algumas estações coletam mais de um registro. Então, nesse caso, o que eu gostaria era de trabalhar com a médias desses dados coletados.

Ou seja, o que preciso é agrupar por coordenadas geográficas, para em seguida tirar a média da variável “v1”

Importar dataset

Então, primeiro importo o dataset…
Depois o converto para o formato “sf” (com a função st_as_sf())
Depois eu crio uma coluna (com a st_as_text()) que nada mais é do que as informações da coluna “geometry” no formato character (o motivo eu explico logo mais)
Em seguida, eu crio um identificador (coluna “id”) para cada linha (o motivo eu também explico mais adiante)
Por fim, seleciono apenas as minhas colunas de interesse

dataset <- import("dataset.xlsx", setclass = "tibble")

Converter para formato “sf” e acrescentar coluna de identificação

dataset_sf <- dataset %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4674) %>% 
  mutate(geometry_text = st_as_text(geometry),
         id = row_number()) %>% 
  select(id, v1, geometry_text)

Resultado:

dataset_sf


Simple feature collection with 2699 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -44.82357 ymin: -15.31984 xmax: 0 ymax: 0
Geodetic CRS:  SIRGAS 2000
# A tibble: 2,699 x 4
      id    v1 geometry_text                            geometry
   <int> <dbl> <chr>                                 <POINT [°]>
 1     1    40 POINT (-44.82357 -15.28918) (-44.82357 -15.28918)
 2     2    66 POINT (-44.82351 -15.28921) (-44.82351 -15.28921)
 3     3    56 POINT (-44.82351 -15.28921) (-44.82351 -15.28921)
 4     4    33 POINT (-44.82351 -15.28921) (-44.82351 -15.28921)
 5     5    53 POINT (-44.82351 -15.28921) (-44.82351 -15.28921)
 6     6    90 POINT (-44.82344 -15.28923) (-44.82344 -15.28923)
 7     7    60 POINT (-44.82344 -15.28923) (-44.82344 -15.28923)
 8     8    56 POINT (-44.82344 -15.28923) (-44.82344 -15.28923)
 9     9    76 POINT (-44.82338 -15.28925) (-44.82338 -15.28925)
10    10    70 POINT (-44.82338 -15.28925) (-44.82338 -15.28925)
# ... with 2,689 more rows

Definir os grupos

Como meu objetivo é tirar uma média para os mesmos pontos de coleta, eu:

usei a função janitor::tabyl() para obter uma relação única de cada um dos pontos.
Como não consigo aplicar a função tabyl () à coluna de geometria, apliquei-a à coluna geometry_text. Por isso que precisei converter as coordenadas para o formato texto;
Na sequência, criei uma coluna chamada “grupo” e a preenchi com o número da linha.
Por fim, salvei como um novo dataset

grupos_dataset_sf <- dataset_sf %>% 
  tabyl(geometry_text) %>% 
  mutate(grupo = row_number()) %>% 
  select(geometry_textgeometry, grupo)

grupos_dataset_sf %>% 
  head()

     geometry_textgeometry grupo
 POINT (-44.76423 -15.31983)     1
 POINT (-44.76424 -15.31983)     2
 POINT (-44.76425 -15.31982)     3
 POINT (-44.76425 -15.31983)     4
 POINT (-44.76425 -15.31984)     5
 POINT (-44.76426 -15.31981)     6

Adicionar os grupos ao dataset “sf”

Com cada ponto associado à um grupo, eu faço um join usando como referência a coluna do “geometry_text”.

dataset_sf_grp <- dataset_sf %>% 
  left_join(grupos_dataset_sf, by = c("geometry_text"="geometry_textgeometry"))

Como resultado tenho cada ponto do dataset original associado a um grupo

Desta forma, consigo agora agrupar pela variável “grupo” e aplicar média

Para efeito da checagem que pretendo fazer ao final, opto por excluir a coluna das geometrias com st_drop_geometry():

dataset_sf_grp_ITAMAR <- 
dataset_sf_grp %>% 
  group_by(grupo) %>% 
  mutate(media_ITAMAR = mean(v1)) %>% 
  ungroup() %>% 
  st_drop_geometry()

Consegui então um dataset com as médias para cada ponto

dataset_sf_grp_ITAMAR


# A tibble: 2,699 x 5
      id    v1 geometry_text               grupo media_ITAMAR
 * <int> <dbl> <chr>                       <int>        <dbl>
 1     1    40 POINT (-44.82357 -15.28918)   577         40  
 2     2    66 POINT (-44.82351 -15.28921)   576         52  
 3     3    56 POINT (-44.82351 -15.28921)   576         52  
 4     4    33 POINT (-44.82351 -15.28921)   576         52  
 5     5    53 POINT (-44.82351 -15.28921)   576         52  
 6     6    90 POINT (-44.82344 -15.28923)   575         68.7
 7     7    60 POINT (-44.82344 -15.28923)   575         68.7
 8     8    56 POINT (-44.82344 -15.28923)   575         68.7
 9     9    76 POINT (-44.82338 -15.28925)   574         83  
10    10    70 POINT (-44.82338 -15.28925)   574         83  
# ... with 2,689 more rows

-------------------------------------------------------------------------

Acontece que vi esse post do Pebesma onde ele sugere uma forma de agrupar e fazer manipulações a partir de datasets do tipo “sf”

Aplicando a sugestão dele ao dataset em questão:

dataset_sf_grp_PEBESMA <- 
dataset_sf %>% 
 mutate(grupo = st_equals(.) %>% 
        map(min) %>% 
        unlist()) %>% 
 group_by(grupo) %>% 
 mutate(media_PEBESMA = mean(v1)) %>% 
 st_drop_geometry()

Observe que também removi a coluna das geometrias para poder fazer o join com o dataset anterior.
Então tenho um dataset com as médias. Só que dessa vez seguindo outro método.

dataset_sf_grp_PEBESMA

# A tibble: 2,699 x 5
# Groups:   grupo [694]
      id    v1 geometry_text               grupo media_PEBESMA
 * <int> <dbl> <chr>                       <int>         <dbl>
 1     1    40 POINT (-44.82357 -15.28918)     1          40  
 2     2    66 POINT (-44.82351 -15.28921)     2          52  
 3     3    56 POINT (-44.82351 -15.28921)     2          52  
 4     4    33 POINT (-44.82351 -15.28921)     2          52  
 5     5    53 POINT (-44.82351 -15.28921)     2          52  
 6     6    90 POINT (-44.82344 -15.28923)     6          68.7
 7     7    60 POINT (-44.82344 -15.28923)     6          68.7
 8     8    56 POINT (-44.82344 -15.28923)     6          68.7
 9     9    76 POINT (-44.82338 -15.28925)     9          83  
10    10    70 POINT (-44.82338 -15.28925)     9          83  
# ... with 2,689 more rows

Checagem

O que acontece no entanto, é que as médias estão diferentes e não entendo o porquê.

Primeiro: faço um join entre o meu dataset com as médias, e o dataset utilizando o método do Pebesma.
Depois crio um coluna para verificar se as médias seguindo meu método e as médias seguindo o método do Pebesma são iguais.

checagem <- dataset_sf_grp_ITAMAR %>% 
  inner_join(dataset_sf_grp_PEBESMA, by = "id") %>% 
  select(id, v1.x, media_ITAMAR, v1.y, media_PEBESMA) %>% 
  mutate(igual = if_else(condition = media_ITAMAR == media_PEBESMA,
                         true = "SIM",
                         false = "NAO"))

# A tibble: 2,699 x 6
      id  v1.x media_ITAMAR  v1.y media_PEBESMA igual
   <int> <dbl>        <dbl> <dbl>         <dbl> <chr>
 1     1    40         40      40          40   SIM  
 2     2    66         52      66          52   SIM  
 3     3    56         52      56          52   SIM  
 4     4    33         52      33          52   SIM  
 5     5    53         52      53          52   SIM  
 6     6    90         68.7    90          68.7 SIM  
 7     7    60         68.7    60          68.7 SIM  
 8     8    56         68.7    56          68.7 SIM  
 9     9    76         83      76          83   SIM  
10    10    70         83      70          83   SIM  
# ... with 2,689 more rows

Quando faço a contagem verifico que há 617 ocorrências para o “não”

checagem %>% 
  count(igual)

# A tibble: 2 x 2
  igual     n
  <chr> <int>
1 NAO     617
2 SIM    2082

Analisando mais de perto

checagem %>% 
  filter(igual == "NAO")

# A tibble: 617 x 6
      id  v1.x media_ITAMAR  v1.y media_PEBESMA igual
   <int> <dbl>        <dbl> <dbl>         <dbl> <chr>
 1   944    16         29.7    16          26   NAO  
 2   945    36         29.7    36          26   NAO  
 3   946    26         29.7    26          26   NAO  
 4   947    63         35.3    63          46.3 NAO  
 5   948    43         35.3    43          46.3 NAO  
 6   949    33         35.3    33          46.3 NAO  
 7   950    46         37.8    46          39.2 NAO  
 8   951    43         37.8    43          39.2 NAO  
 9   952    40         37.8    40          39.2 NAO  
10   953    33         37.8    33          36.6 NAO  
# ... with 607 more rows

Resolvi analisar os registros de 944:946, onde minha média é 29.7 e a do Pebesma é 26.

Vejo que esses registros estão associados ao grupo “334”

dataset_sf_grp_ITAMAR %>% 
  filter(id %in% 944:946)

# A tibble: 3 x 5
     id    v1 geometry_text               grupo media_ITAMAR
  <int> <dbl> <chr>                       <int>        <dbl>
1   944    16 POINT (-44.79965 -15.30073)   334         29.7
2   945    36 POINT (-44.79965 -15.30073)   334         29.7
3   946    26 POINT (-44.79965 -15.30073)   334         29.7
>

Se eu filtrar o grupo “334”:

dataset_sf_grp_ITAMAR %>% 
 filter(grupo == 334)

# A tibble: 17 x 5
     id    v1 geometry_text               grupo media_ITAMAR
  <int> <dbl> <chr>                       <int>        <dbl>
1   944    16 POINT (-44.79965 -15.30073)   334         29.7
2   945    36 POINT (-44.79965 -15.30073)   334         29.7
3   946    26 POINT (-44.79965 -15.30073)   334         29.7
4   973    33 POINT (-44.79965 -15.30073)   334         29.7
5   974    16 POINT (-44.79965 -15.30073)   334         29.7
6   975    23 POINT (-44.79965 -15.30073)   334         29.7
7   976    23 POINT (-44.79965 -15.30073)   334         29.7
8   977    36 POINT (-44.79965 -15.30073)   334         29.7
9   978    46 POINT (-44.79965 -15.30073)   334         29.7
10   979    36 POINT (-44.79965 -15.30073)   334         29.7
11   990    26 POINT (-44.79965 -15.30073)   334         29.7
12  1025    36 POINT (-44.79965 -15.30073)   334         29.7
13  1026    26 POINT (-44.79965 -15.30073)   334         29.7
14  1027    30 POINT (-44.79965 -15.30073)   334         29.7
15  1028    36 POINT (-44.79965 -15.30073)   334         29.7
16  1029    30 POINT (-44.79965 -15.30073)   334         29.7
17  1030    30 POINT (-44.79965 -15.30073)   334         29.7

Tenho o seguinte ponto: “POINT (-44.79965 -15.30073)”

Se eu fizer o mesmo para o dataset das médias do Pebesma:

Os registros estão associados ao grupo “944”, cujo ponto também é: “POINT (-44.79965 -15.30073)”

dataset_sf_grp_PEBESMA %>% 
  filter(id %in% 944:946)

# A tibble: 3 x 5
# Groups:   grupo [1]
     id    v1 geometry_text               grupo media_PEBESMA
  <int> <dbl> <chr>                       <int>         <dbl>
1   944    16 POINT (-44.79965 -15.30073)   944            26
2   945    36 POINT (-44.79965 -15.30073)   944            26
3   946    26 POINT (-44.79965 -15.30073)   944            26

No entanto, se filtro por esse ponto em específico:

dataset_sf_grp_PEBESMA %>% 
  filter(geometry_text == "POINT (-44.79965 -15.30073)")

# A tibble: 17 x 5
# Groups:   grupo [5]
      id    v1 geometry_text               grupo media_PEBESMA
   <int> <dbl> <chr>                       <int>         <dbl>
 1   944    16 POINT (-44.79965 -15.30073)   944          26  
 2   945    36 POINT (-44.79965 -15.30073)   944          26  
 3   946    26 POINT (-44.79965 -15.30073)   944          26  
 4   973    33 POINT (-44.79965 -15.30073)   973          30.4
 5   974    16 POINT (-44.79965 -15.30073)   973          30.4
 6   975    23 POINT (-44.79965 -15.30073)   973          30.4
 7   976    23 POINT (-44.79965 -15.30073)   973          30.4
 8   977    36 POINT (-44.79965 -15.30073)   973          30.4
 9   978    46 POINT (-44.79965 -15.30073)   973          30.4
10   979    36 POINT (-44.79965 -15.30073)   973          30.4
11   990    26 POINT (-44.79965 -15.30073)   990          26  
12  1025    36 POINT (-44.79965 -15.30073)  1025          30.7
13  1026    26 POINT (-44.79965 -15.30073)  1025          30.7
14  1027    30 POINT (-44.79965 -15.30073)  1025          30.7
15  1028    36 POINT (-44.79965 -15.30073)  1028          32  
16  1029    30 POINT (-44.79965 -15.30073)  1028          32  
17  1030    30 POINT (-44.79965 -15.30073)  1028          32

O resultado gera 5 grupos (“944”, “973”, “990”, “1025”, “1028”) para um mesmo ponto.

Fiquei sem entender.

Em suma, o que eu queria saber é como agrupar e tirar médias de datasets do tipo “sf”.
Alías, creio, inclusive, que o subterfúgio que criei (de criar um coluna com as coordenadas no formato de character) não seja nem mesmo a mais adequada. Talvez se possa fazer isso a partir do própria coluna de geometrias.

Se alguém puder ajudar, agradeceria bastante.

denis-or · Agosto 17, 2022, 3:05pm

Oi Itamar.

você trouxe muitos pontos na sua dúvida e eu quero fazer alguns apontamentos, apenas.

será que as médias são diferentes porque os grupos são diferentes? Veja abaixo:
o “grupo” 1 do df dataset_sf_grp_PEBESMA é o ponto (-44.82357 -15.28918). O mesmo ponto no df dataset_sf_grp_ITAMAR é o “grupo” 577. O que pode ser visto no df grupos_dataset_sf também. Se usar a função tabyl() do janitor ou a count() do dplyr, ambas vão organizar o resultado do menor pro maior e, dessa forma, que você definiu o grupo.
Já no dataset_sf_grp_PEBESMA você usou um predicado de relacionamento espacial dado pela função st_equals(), que significa que foi verificado qual geometria espacial estava totalmente contida por outra. Essa função identifica as geometrias pelo id do registro. Assim, neste momento, você escolheu pra identificar o grupo o menor id (função min()). Ou seja, no dataset_sf_grp_ITAMAR você identificou o grupo pelo geometry_text e no dataset_sf_grp_PEBESMA pela posição da geometria no banco (id da linha). Se agrupar em ambas as tabela pelo geometry_text talvez funcione.
Será que isso que está fazendo está certo?
Bom, como geógrafo, gostaria de lembrá-lo que há uma questão de precisão nas coordenadas espaciais nesta tabela. Quando você identifica o ponto (-44.82357 -15.28918) como sendo um “grupo” diferente do ponto (-44.82351 -15.28921), tenho a seguinte dúvida: Tem certeza que não é o mesmo ponto? Estão muito próximos. Ao dizer “As latitudes e longitudes representam estações que coletam dados”, imagino que os pontos representam a localização da estação. Não acho que teriam duas estações tão próximas. Experimente jogar este ponto no google maps pra verificar. Quando diz estações, seriam meteorológicas? A geometria parece que foi uma medição ao longo de uma estrada:
plot(st_geometry(filter(dataset_sf, geometry_text != "POINT (0 0)")))

image694×428 1.71 KB

Penso que o ideal (pra facilitar apenas) seria fazer os grupos pela “estação” mesmo (o nome ou identificação dela). Caso não consiga, tente utilizar um raio pra cada ponto e verificar quais estão próximas o suficiente para classificar que é uma só. Caso seja uma medição ao longo de uma estrada, talvez tirar a média a cada x metros (ou quilômetros) também faça algum sentido.

Bom. Espero ter ajudado de alguma forma.

Valeu!