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
DataTemp.csv
Description: MS-Excel spreadsheet
_______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo