[R-sig-Geo] Point Interpolation and Analysis

[sownal chand]

Hello Sir/madam,

I have been trying to interpolate this data and make an animation by
following the codes given by Devan using Dr Timo's Dialect analysis in
Germany,
https://github.com/devanmcg/IntroRangeR/blob/master/12_RasGIS/Interpolation.R
.
Is there a way in which these code could be used for Mapping in
Pacific Region. Below is the edited version with errors

#packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, sf, kknn)
#
# Create a shorthand reference to Fiji Islands CRS
crs_etrs = "+proj=merc +lon_0=150 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m
+no_defs"
#
# Fiji geo data
# download
#Draft map of Fiji
world <- rnaturalearth::ne_countries(scale = "Large", returnclass = "sf")

Fiji <- ggplot(data=world) +
  geom_sf() +
  coord_sf(
    crs = 3832, # https://epsg.io/3832
   xlim = c(2984265.06, 3539584.72), # limits are taken from projected
bounds
   ylim = c(-2224162.41, -1811688.88)  # of EPSG:3832
  )+ theme_bw()

#Extract Fiji Boundary data
Fiji <- FJI %>%
  filter(iso_a3 == "FJI")
plot(Fiji)



Fiji <- Fiji %>%
  group_by(iso_a3) %>%
  summarize()
# simplify
Fiji %>%
  st_simplify(preserveTopology = T, dTolerance = 0.01)

#load point data from the data file stored
data_53YearTemp <- read.csv("C://Users/Sownal/Documents/53YearsTemp.csv")
View(data_53YearTemp)

data$long <- as.numeric(data$Longitude)
data$lat <- as.numeric(data$Latitude)

#remove NA valuues in the spatial Data Frame
data_53YearTemp <- na.omit(data_53YearTemp)

#convert 53YearTemp to a Dataframe
#spatial points dataframe
data.sp1 <- SpatialPointsDataFrame(data_53YearTemp[,c(3,4)],
data_53YearTemp[,-c(3,4)])
View(data.sp1)
#point data to sf
data.sp1 %>%
  st_as_sf(coords = c("lng", "lat"),
           crs = "+proj=merc +lon_0=150 +k=1 +x_0=0 +y_0=0 +datum=WGS84
+units=m +no_defs ")


#reproject data to the specificed crs or epsg country code
# CRS: ETRS
etrs <- "+proj=merc +lon_0=150 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m
+no_defs "
Fiji %>%
  st_transform(etrs)
data.sp_train %>%
  st_transform(etrs)

#remove NA valuues in the spatial Data Frame
data_53YearTemp <- na.omit(data_53YearTemp)

#cliping dataframe to buffer Fiji
Fiji_buffered <- Fiji %>%
  st_buffer(dist = 10000) # in meters, 10km buffer
data.sp1 <- data.sp1[Fiji_buffered, ]

# specify grid width in pixels
width_in_pixels <- 300
# dx is the width of a grid cell in meters
dx <- ceiling( (st_bbox(Fiji_buffered)["xmax"] -
                  st_bbox(Fiji_buffered)["xmin"]) / width_in_pixels)
# dy is the height of a grid cell in meters
# because we use quadratic grid cells, dx == dy
dy <- dx
# calculate the height in pixels of the resulting grid
height_in_pixels <- floor( (st_bbox(Fiji_buffered)["ymax"] -
                              st_bbox(Fiji_buffered)["ymin"]) / dy)

grid <- st_make_grid(Fiji_buffered,
                     cellsize = dx,
                     n = c(width_in_pixels, height_in_pixels),
                     what = "centers"
)

rm(dx, dy, height_in_pixels, Fiji_buffered)

plot(grid)
# View grid
ggplot() +
  geom_sf(data=grid) +
  geom_sf(data=Fiji, fill=NA, color="blue", lwd=2)
# Prepare data for interpolation
#
# Create tibble of the TempData sf object
TEmp_input <- data_53YearTemp %>%
  tibble::as.tibble(data_53YearTemp = .$pronunciation_id,
         lon = st_coordinates(.)[, 4],
         lat = st_coordinates(.)[, 3]) %>%
  select(Year, lon, lat)
#training data set for the interpolations
data.sp_train <- SpatialPointsDataFrame(data_53YearTemp[,c(3,4)],
data_53YearTemp[,-c(3,4)])
View(data.sp_train)



# config
k <- 1000 # "k" for k-nearest-neighbour-interpolation

# specify function which is executed for each tile of the grid
computeGrid <- function(grid, data.sp1, knn) {
  # create empty result data frame
  Temp_result <- data.frame(Temp = as.factor(NA),
                                lon = st_coordinates(grid)[, 1],
                                lat = st_coordinates(grid)[, 2])
  # run KKNN
  Temp_kknn <- kknn::kknn(data.sp1_train ~ .,
                              train = data.sp1_train,
                              test = Temp_result,
                              kernel = "gaussian",
                              k = knn)
  # bring back to result data frame
  # only retain the probability of the dominant dialect at that grid cell
  Temp_result %>%
    # extract the interpolated dialect at each grid cell with the
    # kknn::fitted function
    mutate(Temp = fitted(Temp_kknn),
           # only retain the probability of the interpolated dialect,
           # discard the other 7
           prob = apply(Temp_kknn$prob,
                        1,
                        function(x) max(x)))
  return(Temp_result)
}

# specify the number of cores below (adapt if you have fewer cores or
# want to reserve some computation power to other stuff)
#registerDoParallel(cores = 5)

# specify number of batches and resulting size of each batch (in grid cells)
no_batches <- 60 # increase this number if you run into memory problems

batch_size <- ceiling(length(grid) / no_batches)

# PARALLEL COMPUTATION
start.time <- Sys.time()
Temp_result <- foreach(
  batch_no = 1:no_batches,
  # after each grid section is computed, rbind the resulting df into
  # one big dialects_result df
  .combine = rbind,
  # the order of grid computation doesn't matter: this speeds it up even
more
  .inorder = FALSE) %dopar% {
    # compute indices for each grid section, depending on batch_size and
current
    # batch
    start_idx <- (batch_no - 1) * batch_size + 1
    end_idx <- batch_no * batch_size
    # specify which section of the grid to interpolate, using regular
    # subsetting
    grid_batch <- grid[start_idx:ifelse(end_idx > length(grid),
                                        length(grid),
                                        end_idx)]
    # apply the actual computation to each grid section
    df <- computeGrid(grid_batch, data.sp1_train, k)
  }
end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken


# convert resulting df back to sf object, but do not remove raw geometry
cols
Temp_raster <- st_as_sf(data.sp_train,
                            coords = c("lon", "lat"),
                            crs = etrs,
                            remove = F)
# clip raster to Fiji again
Temp_raster <- Temp_raster[Fiji, ]
rm(Temp_result, k)

Sample data is also provided. Any assistance would be greatly appreciated.

Yours sincerely
sownalc

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20220315/856f8f3b/attachment.html>

-------------- next part --------------
A non-text attachment was scrubbed...
Name: DataTemp.csv
Type: application/vnd.ms-excel
Size: 5203 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20220315/856f8f3b/attachment.xlb>