# [R-sig-Geo] Point Interpolation and Analysis

sownal chand @own@|ch@nd|m@ @end|ng |rom gm@||@com
Tue Mar 15 06:37:52 CET 2022

```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>
```

More information about the R-sig-Geo mailing list