[R] Raster map in R

Tue Mar 8 10:11:02 CET 2022

On 08/03/2022 08:21, Eliza Botto wrote:
> deaR expeRts,
>
> I have the following data
>
>> dput(Tuto)
> structure(list(X = c(-114.028, -114.011, -114.442, -113.937,
> -114.187, -114.083, -113.949, -114.15, -114.068, -114.203, -113.958,
> -114.248, -114.18, -114.14, -114.071, -114.042, -114.187, -114.03,
> -113.97, -113.824, -114.084, -114.152, -114.468, -113.935, -113.994,
> -114.048, -114.188, -114.129, -114.071, -113.934, -114.001, -114.075,
> -114.121, -113.958, -114.039, -113.963, -114.062, -114.183, -114.118,
> -114.119, -113.954, -114.051, -113.988, -114.194, -114.025),
>      Y = c(51.431, 51.279, 51.167, 51.165, 51.155, 51.14, 51.126,
>      51.126, 51.116, 51.115, 51.092, 51.091, 51.084, 51.082, 51.076,
>      51.069, 51.062, 51.052, 51.051, 51.048, 51.044, 51.039, 51.037,
>      51.034, 51.03, 51.029, 51.021, 51.006, 51.005, 50.99, 50.983,
>      50.966, 50.958, 50.948, 50.929, 50.916, 50.911, 50.908, 50.908,
>      50.907, 50.877, 50.877, 50.86, 50.854, 50.826), a = c(0.838,
>      0.901, 0.953, 0.902, 0.782, 0.938, 0.884, 0.879, 0.932, 0.947,
>      0.965, 0.828, 0.923, 0.892, 0.884, 0.897, 0.912, 0.988, 0.901,
>      0.855, 0.999, 0.846, 0.845, 0.798, 0.749, 0.753, 0.762, 0.646,
>      0.729, 0.544, 0.265, 0.449, 0.334, 0.36, 0.325, 0.337, 0.249,
>      0.114, 0.149, 0.173, 0.175, 0.184, 0.219, 0.106, 0.148),
>      n = c(0.653, 0.625, 0.641, 0.63, 0.656, 0.619, 0.628, 0.634,
>      0.63, 0.604, 0.598, 0.617, 0.632, 0.635, 0.637, 0.646, 0.619,
>      0.613, 0.63, 0.615, 0.604, 0.639, 0.598, 0.593, 0.583, 0.606,
>      0.594, 0.653, 0.63, 0.577, 0.624, 0.626, 0.676, 0.641, 0.629,
>      0.579, 0.603, 0.607, 0.66, 0.614, 0.618, 0.574, 0.552, 0.62,
>      0.599)), row.names = c(19L, 22L, 18L, 3L, 45L, 20L, 14L,
> 43L, 40L, 44L, 37L, 12L, 42L, 41L, 39L, 38L, 33L, 8L, 36L, 16L,
> 32L, 31L, 5L, 34L, 35L, 9L, 13L, 30L, 29L, 4L, 24L, 27L, 28L,
> 23L, 25L, 21L, 26L, 6L, 15L, 2L, 1L, 7L, 11L, 10L, 17L), class = "data.frame")

Here's my attempt:

I used the terra library.

First I visually checked the distribution of the points to show that you 
should use interpolation to convert the points to a raster.

Then, as an example I used the terra::interpolate function to get a 
raster interpolation using the simple IDW from gstat. This **might not** 
be the correct interpolation method. That is something you will have to 
determine.

Note: I changed your X and Y columns to lowercase to stay in line with 
the examples in gstat

library(terra)

# Points data.frame
Tuto <- structure(list(x = c(-114.028, -114.011, -114.442, -113.937,
-114.187, -114.083, -113.949, -114.15, -114.068, -114.203, -113.958,
-114.248, -114.18, -114.14, -114.071, -114.042, -114.187, -114.03,
-113.97, -113.824, -114.084, -114.152, -114.468, -113.935, -113.994,
-114.048, -114.188, -114.129, -114.071, -113.934, -114.001, -114.075,
-114.121, -113.958, -114.039, -113.963, -114.062, -114.183, -114.118,
-114.119, -113.954, -114.051, -113.988, -114.194, -114.025),
     y = c(51.431, 51.279, 51.167, 51.165, 51.155, 51.14, 51.126,
     51.126, 51.116, 51.115, 51.092, 51.091, 51.084, 51.082, 51.076,
     51.069, 51.062, 51.052, 51.051, 51.048, 51.044, 51.039, 51.037,
     51.034, 51.03, 51.029, 51.021, 51.006, 51.005, 50.99, 50.983,
     50.966, 50.958, 50.948, 50.929, 50.916, 50.911, 50.908, 50.908,
     50.907, 50.877, 50.877, 50.86, 50.854, 50.826), a = c(0.838,
     0.901, 0.953, 0.902, 0.782, 0.938, 0.884, 0.879, 0.932, 0.947,
     0.965, 0.828, 0.923, 0.892, 0.884, 0.897, 0.912, 0.988, 0.901,
     0.855, 0.999, 0.846, 0.845, 0.798, 0.749, 0.753, 0.762, 0.646,
     0.729, 0.544, 0.265, 0.449, 0.334, 0.36, 0.325, 0.337, 0.249,
     0.114, 0.149, 0.173, 0.175, 0.184, 0.219, 0.106, 0.148),
     n = c(0.653, 0.625, 0.641, 0.63, 0.656, 0.619, 0.628, 0.634,
     0.63, 0.604, 0.598, 0.617, 0.632, 0.635, 0.637, 0.646, 0.619,
     0.613, 0.63, 0.615, 0.604, 0.639, 0.598, 0.593, 0.583, 0.606,
     0.594, 0.653, 0.63, 0.577, 0.624, 0.626, 0.676, 0.641, 0.629,
     0.579, 0.603, 0.607, 0.66, 0.614, 0.618, 0.574, 0.552, 0.62,
     0.599)), row.names = c(19L, 22L, 18L, 3L, 45L, 20L, 14L,
43L, 40L, 44L, 37L, 12L, 42L, 41L, 39L, 38L, 33L, 8L, 36L, 16L,
32L, 31L, 5L, 34L, 35L, 9L, 13L, 30L, 29L, 4L, 24L, 27L, 28L,
23L, 25L, 21L, 26L, 6L, 15L, 2L, 1L, 7L, 11L, 10L, 17L), class = 
"data.frame")

Tuto_pts <- vect(Tuto, geom=c("x", "y"), crs="EPSG:4326")
plot(Tuto_pts)
# Clearly not evenly distributed, interpolation needed

# Prepare target raster, with given extent and resolution of 0.01 degree
Tuto_ext <- ext(-114.5, -113.7, 50.800, 51.600)
res = 0.01
Tuto_rast <- rast(crs = "EPSG:4326", extent=Tuto_ext, resolution=res, 
vals=0)

# IDW model from gstat package

library(gstat)
model_idw_a <- gstat(id = "a", formula = a~1, locations=~x+y, data=Tuto,
                    nmax=7, set=list(idp = .5))
Tuto_a <- interpolate(Tuto_rast, model_idw_a, debug.level=0, index=1)

model_idw_n <- gstat(id = "n", formula = n~1, locations=~x+y, data=Tuto,
                    nmax=7, set=list(idp = .5))
Tuto_n <- interpolate(Tuto_rast, model_idw_n, debug.level=0, index=1)

# Merge two interpolations into one multiband raster
Tuto_multiband = c(Tuto_a, Tuto_n)
plot(Tuto_multiband)

HTH,

Micha

> In the given data,
>
> X-> Latitude
>
> Y-> Longitude
>
> a-> Factor 1
>
> n-> Factor 2
>
> I want to draw a raster map of these values within the following limits  (xmn=-114.5, xmx=-113.7, ymn=50.800, ymx=51.600). I have tried through "raster" library but failed as I wasn't able to properly generate the data.
> I will be thankful if I can get any kind of help. Thank-you very much in advance.
>
> Eliza
>
>
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918