Plotting Spatial Points Data Frame

[Jouline]

I would like to plot the prediction I obtained from my heavy metal study.
I got this object Ni.krig that contains the prediction from krigage,
and this clean.donnees.Ni represent my data points.

I would like to plot a map with both prediction grid and measurement points.

As an example I made this with IDW prediction before.

But in order to make krigage possible, I had to switch from my initial dataframe to a SpatialPointsDataFrame (also named S4) :

> typeof(clean.donnees.Ni)
[1] "S4"
> typeof(Ni.krig)
[1] "S4"

> str(Ni.krig)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame':	546 obs. of  2 variables:
  .. ..$ var1.pred: num [1:546] 16.6 16.6 16.6 16.6 16.6 ...
  .. ..$ var1.var : num [1:546] 28.8 28.8 28.8 28.8 28.8 ...
  ..@ coords.nrs : int [1:2] 1 2
  ..@ coords     : num [1:546, 1:2] 45553 50553 55553 60553 65553 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:546] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "x" "y"
  ..@ bbox       : num [1:2, 1:2] 45553 68096 170553 168096
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr NA

> str(clean.donnees.Ni) 
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame':	2593 obs. of  12 variables:
  .. ..$ FID_1     : int [1:2593]  587 537 222 552 477 209 13 554 72 221 ...
  .. ..$ ID        : chr [1:2593] "02/LA06306" "02/LA05627" "00941J" "02/LA05824" ...
  .. ..$ Ni        : num [1:2593] 8 9 10.4 16 14 11.3 9.05 11 15.1 9.63 ...
  .. ..$ Zn        : num [1:2593] 27.3 28.6 52.1 54.9 60.2 ...
  .. ..$ Cr        : num [1:2593] 16.9 25 20.5 25 27 ...
  .. ..$ SIGLE_PEDO: chr [1:2593] "(u)-Ldc2_3" "(x)Ldc" "(x)Lda1" "(x)Ldc" ...
  .. ..$ region_etm: int [1:2593] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ GS        : int [1:2593] 3 3 3 3 3 3 4 3 3 3 ...
  .. ..$ Cdbx.Ni   : num [1:2593] 2.99 3.23 3.53 4.54 4.21 ...
  .. ..$ Cdbx.Zn   : num [1:2593] 4.2 4.28 5.28 5.37 5.53 ...
  .. ..$ Cdbx.Cr   : num [1:2593] 5.62 7.1 6.31 7.1 7.42 ...
  .. ..$ Ni.res    : num [1:2593] -4.16 -3.51 -3.15 3.26 1.25 ...
  ..@ coords.nrs : int [1:2] 2 3
  ..@ coords     : num [1:2593, 1:2] 87194 90548 75270 90628 89151 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2593] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "x" "y"
  ..@ bbox       : num [1:2, 1:2] 70855 73121 166591 161300
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr NA

# Im only making krigage with the $ Ni in the data

I realised I could use ggplot2 anymore because S4 is a different object type.

So I investigate a little bit and found this question :

Plotting spatial points data in R?

I tried this :

library(sp)
library(sf)
library(ggmap)
proj4string(Ni.krig) <- CRS("+init=epsg:28992")


Ni.krig <- st_as_sf(Ni.krig,coords = 1:2)

Ni.krig <- Ni.krig %>% 
  st_transform(crs = 4326)

Krig_map <- get_stamenmap(
  bbox = unname(st_bbox(Ni.krig)),
  zoom = 3, maptype = 'toner-lite', source = 'stamen') %>% ggmap()

Krig_map + 
  geom_sf(
    data = Ni.krig, 
    aes(size = var1.pred), 
    color = 'red', alpha = 0.5,
    show.legend = 'point', inherit.aes = F
  )

but it gaves me this as result :

Of course, my CRS is 31370 so the location is not good but I wonder what I should put in this line : Ni.krig <- Ni.krig %>% st_transform(crs = ???)

I am more familiar with ggplot2 so maybe is you know a way to use this package instead of those manipulations ?

Something like :

ggplot() + 
  geom_tile(data=Ni.krig@data, 
            aes(x=Ni.krig@coords[1,] , y=Ni.krig@coords[2,], fill = Ni.predkrig)) 

Thank you !

[technocrat]

Consider using the sf package to create a simple features data frame and then plotting with ggplot, using the geom_sf() feature.

suppressPackageStartupMessages({
  library(ggplot2)
  library(sf)
})

demo(nc, ask = FALSE, echo = FALSE)
#> Reading layer `nc.gpkg' from data source 
#>   `/home/roc/R/x86_64-pc-linux-gnu-library/4.1/sf/gpkg/nc.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 100 features and 14 fields
#> Attribute-geometry relationship: 0 constant, 8 aggregate, 6 identity
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
ggplot() + 
  geom_sf(data = nc, aes(fill = BIR74)) + 
  scale_y_continuous(breaks = 34:36)