<div dir="ltr">Hello Sir/madam, <div><br></div><div>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, <a href="https://github.com/devanmcg/IntroRangeR/blob/master/12_RasGIS/Interpolation.R">https://github.com/devanmcg/IntroRangeR/blob/master/12_RasGIS/Interpolation.R</a> . </div><div>Is there a way in which these code could be used for Mapping in Pacific Region. Below is the edited version with errors</div><div><br></div><div>#packages <br>if (!require("pacman")) install.packages("pacman")<br>pacman::p_load(tidyverse, sf, kknn)<br>#<br># Create a shorthand reference to Fiji Islands CRS<br>crs_etrs = "+proj=merc +lon_0=150 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"<br>#<br># Fiji geo data<br># download<br>#Draft map of Fiji<br>world <- rnaturalearth::ne_countries(scale = "Large", returnclass = "sf")<br><br>Fiji <- ggplot(data=world) +<br>  geom_sf() + <br>  coord_sf(<br>    crs = 3832, # <a href="https://epsg.io/3832">https://epsg.io/3832</a><br>   xlim = c(2984265.06, 3539584.72), # limits are taken from projected bounds<br>   ylim = c(-2224162.41, -1811688.88)  # of EPSG:3832<br>  )+ theme_bw()<br><br>#Extract Fiji Boundary data<br>Fiji <- FJI %>% <br>  filter(iso_a3 == "FJI")<br>plot(Fiji)<br><br><br><br>Fiji <- Fiji %>%<br>  group_by(iso_a3) %>% <br>  summarize()<br># simplify<br>Fiji %>% <br>  st_simplify(preserveTopology = T, dTolerance = 0.01)<br><br>#load point data from the data file stored <br>data_53YearTemp <- read.csv("C://Users/Sownal/Documents/53YearsTemp.csv")<br>View(data_53YearTemp)<br><br>data$long <- as.numeric(data$Longitude)<br>data$lat <- as.numeric(data$Latitude)<br><br>#remove NA valuues in the spatial Data Frame<br>data_53YearTemp <- na.omit(data_53YearTemp) <br><br>#convert 53YearTemp to a Dataframe<br>#spatial points dataframe<br>data.sp1 <- SpatialPointsDataFrame(data_53YearTemp[,c(3,4)], data_53YearTemp[,-c(3,4)])<br>View(data.sp1)<br>#point data to sf<br>data.sp1 %>%<br>  st_as_sf(coords = c("lng", "lat"),<br>           crs = "+proj=merc +lon_0=150 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs ")<br><br><br>#reproject data to the specificed crs or epsg country code<br># CRS: ETRS<br>etrs <- "+proj=merc +lon_0=150 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs "<br>Fiji %>%<br>  st_transform(etrs)<br>data.sp_train %>%<br>  st_transform(etrs)<br><br>#remove NA valuues in the spatial Data Frame<br>data_53YearTemp <- na.omit(data_53YearTemp) <br><br>#cliping dataframe to buffer Fiji <br>Fiji_buffered <- Fiji %>%<br>  st_buffer(dist = 10000) # in meters, 10km buffer<br>data.sp1 <- data.sp1[Fiji_buffered, ]<br><br># specify grid width in pixels<br>width_in_pixels <- 300<br># dx is the width of a grid cell in meters<br>dx <- ceiling( (st_bbox(Fiji_buffered)["xmax"] - <br>                  st_bbox(Fiji_buffered)["xmin"]) / width_in_pixels)<br># dy is the height of a grid cell in meters<br># because we use quadratic grid cells, dx == dy<br>dy <- dx<br># calculate the height in pixels of the resulting grid<br>height_in_pixels <- floor( (st_bbox(Fiji_buffered)["ymax"] - <br>                              st_bbox(Fiji_buffered)["ymin"]) / dy)<br><br>grid <- st_make_grid(Fiji_buffered, <br>                     cellsize = dx,<br>                     n = c(width_in_pixels, height_in_pixels),<br>                     what = "centers"<br>)<br><br>rm(dx, dy, height_in_pixels, Fiji_buffered)<br><br>plot(grid)<br># View grid<br>ggplot() +<br>  geom_sf(data=grid) +<br>  geom_sf(data=Fiji, fill=NA, color="blue", lwd=2)<br># Prepare data for interpolation<br>#<br># Create tibble of the TempData sf object<br>TEmp_input <- data_53YearTemp %>%<br>  tibble::as.tibble(data_53YearTemp = .$pronunciation_id, <br>         lon = st_coordinates(.)[, 4], <br>         lat = st_coordinates(.)[, 3]) %>%<br>  select(Year, lon, lat)<br>#training data set for the interpolations<br>data.sp_train <- SpatialPointsDataFrame(data_53YearTemp[,c(3,4)], data_53YearTemp[,-c(3,4)])<br>View(data.sp_train)<br><br><br><br># config<br>k <- 1000 # "k" for k-nearest-neighbour-interpolation<br><br># specify function which is executed for each tile of the grid<br>computeGrid <- function(grid, data.sp1, knn) {<br>  # create empty result data frame<br>  Temp_result <- data.frame(Temp = as.factor(NA), <br>                                lon = st_coordinates(grid)[, 1], <br>                                lat = st_coordinates(grid)[, 2])<br>  # run KKNN<br>  Temp_kknn <- kknn::kknn(data.sp1_train ~ ., <br>                              train = data.sp1_train, <br>                              test = Temp_result, <br>                              kernel = "gaussian", <br>                              k = knn)<br>  # bring back to result data frame<br>  # only retain the probability of the dominant dialect at that grid cell<br>  Temp_result %>%<br>    # extract the interpolated dialect at each grid cell with the <br>    # kknn::fitted function<br>    mutate(Temp = fitted(Temp_kknn),<br>           # only retain the probability of the interpolated dialect,<br>           # discard the other 7<br>           prob = apply(Temp_kknn$prob, <br>                        1, <br>                        function(x) max(x)))<br>  return(Temp_result)<br>}<br><br># specify the number of cores below (adapt if you have fewer cores or<br># want to reserve some computation power to other stuff)<br>#registerDoParallel(cores = 5)<br><br># specify number of batches and resulting size of each batch (in grid cells)<br>no_batches <- 60 # increase this number if you run into memory problems<br><br>batch_size <- ceiling(length(grid) / no_batches)<br><br># PARALLEL COMPUTATION<br>start.time <- Sys.time()<br>Temp_result <- foreach(<br>  batch_no = 1:no_batches, <br>  # after each grid section is computed, rbind the resulting df into<br>  # one big dialects_result df<br>  .combine = rbind, <br>  # the order of grid computation doesn't matter: this speeds it up even more<br>  .inorder = FALSE) %dopar% {<br>    # compute indices for each grid section, depending on batch_size and current<br>    # batch<br>    start_idx <- (batch_no - 1) * batch_size + 1<br>    end_idx <- batch_no * batch_size<br>    # specify which section of the grid to interpolate, using regular<br>    # subsetting<br>    grid_batch <- grid[start_idx:ifelse(end_idx > length(grid), <br>                                        length(grid),<br>                                        end_idx)]<br>    # apply the actual computation to each grid section<br>    df <- computeGrid(grid_batch, data.sp1_train, k)<br>  }<br>end.time <- Sys.time()<br>time.taken <- end.time - start.time<br>time.taken<br><br><br># convert resulting df back to sf object, but do not remove raw geometry cols<br>Temp_raster <- st_as_sf(data.sp_train, <br>                            coords = c("lon", "lat"),<br>                            crs = etrs,<br>                            remove = F)<br># clip raster to Fiji again<br>Temp_raster <- Temp_raster[Fiji, ]<br>rm(Temp_result, k)<br></div><div><br></div><div>Sample data is also provided. Any assistance would be greatly appreciated.</div><div><br></div><div>Yours sincerely </div><div>sownalc</div></div>