R sf calculate average value of each point-CodePudding

I have a forty year monthly timeseries dataset. I would like to create 3 new columns (AvgTMean, AvgTMin, AvgTMax) that would be the average of tmean, Tmin and TMax for each point in the data (each point will have its own unique average value) respectively. I will then map these average values from any of the 3 average value columns on an interactive.

The purpose is to create a map that shows the 40 year average temperature values.

How can I calculate the average of each point?

Sample data (sf):

structure(list(Info = c(NA_character_, NA_character_, NA_character_, 
NA_character_, NA_character_, NA_character_, NA_character_, NA_character_, 
NA_character_), tmean = c(22.2395992279053, 22.7657985687256, 
24.4260005950928, 19.601001739502, 21, 24, 26, 21.45, 27.6), Variable = c("tmean", "tmean", 
    "tmean", "tmean", "tmax", "tmax", "tmax", "tmax", "tmax"), 
     year = c(2021L, 2021L, 1980L, 1980L, 2021L, 2021L, 
    2021L, 2021L, 2021L), month = c(11L, 12L, 0L, 1L, 6L, 7L, 
    8L, 9L, 10L), TMin = c(15, 15.23, 16.12, 13.45, 16.46, 12.11, 11.55, 9.78, 10.56), TMax = c(0, 
    39.69, 40.001, 43.2, 40.6976985931396, 41.7550983428955, 42.1988983154297, 
    41.6512985229492, 40.2621994018555), geometry = structure(list(
        structure(c(-80.2083333327448, 26.2083333333333), class = c("XY", 
        "POINT", "sfg")), structure(c(-80.2083333327448, 26.2083333333333
        ), class = c("XY", "POINT", "sfg")), structure(c(-80.2083333327448, 
        26.2083333333333), class = c("XY", "POINT", "sfg")), 
        structure(c(-80.2083333327448, 26.2083333333333), class = c("XY", 
        "POINT", "sfg")), structure(c(-80.2083333327448, 26.0416666666667
        ), class = c("XY", "POINT", "sfg")), structure(c(-80.2083333327448, 
        26.0416666666667), class = c("XY", "POINT", "sfg")), 
        structure(c(-80.2083333327448, 26.0416666666667), class = c("XY", 
        "POINT", "sfg")), structure(c(-80.2083333327448, 26.0416666666667
        ), class = c("XY", "POINT", "sfg")), structure(c(-80.2083333327448, 
        26.0416666666667), class = c("XY", "POINT", "sfg"))), precision = 0, bbox = structure(c(xmin = -80.2083333327448, 
    ymin = 26.0416666666667, xmax = -80.2083333327448, ymax = 26.2083333333333
    ), class = "bbox"), crs = structure(list(input = "WGS 84", 
        wkt = "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"latitude\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"longitude\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L, class = c("sfc_POINT", 
    "sfc"))), row.names = c(NA, 9L), sf_column = "geometry", agr = structure(c(Info = NA_integer_, 
tmean = NA_integer_, CITYNAME = NA_integer_, Model = NA_integer_, 
Variable = NA_integer_, Datatype = NA_integer_, Resolution = NA_integer_, 
year = NA_integer_, month = NA_integer_, TMin = NA_integer_, 
TMax = NA_integer_), class = "factor", .Label = c("constant", 
"aggregate", "identity")), class = c("sf", "data.frame"))

Code:

library(tidyverse)
library(sf)

    sf_avg = sf %>% 
           summarise(AvgTMean = mean(tmean),
                     AvgTMin = mean(TMin),
                     AvgTMax = mean(TMax)) # Returns the same average value for both the points

# Test static map

qtm(df_df, fill ="AvgTMean", legend = TRUE)  
   tm_add_legend(labels = dff$AvgTMean)

CodePudding user response：

I had trouble with your example so I make a small reproducible example. You will probably need to adapt it for your needs.

# loading lib
library(sf)
library(dplyr)
# data is in tidy format so you have a replicate of 3 points 
pnt_attribute <- data.frame(year = rep(seq(from = 2022 - 40, to = 2021, by = 1),3) 
           , temp1 = c(rnorm(40, mean = 20, sd = 5) 
                    , rnorm(40, mean = 20, sd = 10)
                    , rnorm(40, mean = 15, sd = 5))
 )
# this create a table with 3 points replicate 40 times
pnts <- data.frame(
      id = c(rep(1, 40), 
             rep(2, 40), 
             rep(3,40))
      , x = c(rep(5, 40), 
              rep(1, 40), 
              rep(0, 40))
      , y = c(rep(2, 40), 
              rep(3, 40), 
              rep(1, 40)))
# combine the two df then transorm into sf 
rep_example <- st_as_sf(cbind(pnts,  pnt_attribute), coords = c("x", "y"))

Then you create your summary:

rep_example_agg <- rep_example |> 
                    group_by(id) |>
                    summarize(AvgTMean = mean(temp1, na.rm = TRUE),
                              AvgTMin  = min(temp1, na.rm = TRUE),
                              AvgTMax  = max(temp1, na.rm = TRUE))
# the result 
# Simple feature collection with 3 features and 4 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 0 ymin: 1 xmax: 5 ymax: 3
# CRS:           NA
# # A tibble: 3 × 5
# id AvgTMean AvgTMin AvgTMax geometry
# <dbl>    <dbl>   <dbl>   <dbl>  <POINT>
# 1     1     19.4    5.92    29.2    (5 2)
# 2     2     20.0   -1.34    40.8    (1 3)
# 3     3     14.6    3.37    26.1    (0 1)

I think you either needed the group_by of your data are in an other structure.

Good resource: https://geocompr.robinlovelace.net/attr.html#vector-attribute-aggregation