Raster map in R
If Micha's reply doesn't satisfy, the r-sig-geo list would be a better place to post this. https://stat.ethz.ch/mailman/listinfo/r-sig-geo Bert Gunter "The trouble with having an open mind is that people keep coming along and sticking things into it." -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
On Tue, Mar 8, 2022 at 1:13 AM Micha Silver <tsvibar at gmail.com> wrote:
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 at 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
______________________________________________ R-help at 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.