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 <tsvi...@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@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@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. ______________________________________________ R-help@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.