[R] Raster map in R
Bert Gunter
bgunter@4567 @end|ng |rom gm@||@com
Tue Mar 8 15:54:21 CET 2022
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 using 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 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
>
> ______________________________________________
> 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.
More information about the R-help
mailing list