Combine different projected raster and shapefile into Shiny leaflet with extract?

I have a small issue with Shiny. I hope you can help me out.

I had some issues with using leaflet in combination with Shiny. Plotting maps within a Shiny application using leaflet just works fine for me until now. However, I was faced with a problem using the extract function. I have a tiffile with temperature values and a shapefile with agricultural plots. For each plot I want to compute the mean temperature using extract. Below is a general description followed by the script used.

My input consists of a tiffile and a shapefile. This tiffile has is created in another script and was given a projected reference system called RD New:
+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs (EPSG 28992)

I imported this tiffile as a raster. The projection is NA. But based on the extent one can see that it is actually the RD New projection.

class : RasterLayer
dimensions : 4101, 3499, 14349399 (nrow, ncol, ncell)
resolution : 0.619901, 0.619901 (x, y)
extent : 121988.8, 124157.8, 464956, 467498.2 (xmin, xmax, ymin, ymax)
coord. ref. : NA
...

My second input is a shapefile containing the plots. I imported this one using readOGR.

class : SpatialPolygonsDataFrame
features : 148
extent : 4.904595, 4.936264, 52.17177, 52.19466 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
variables : 11
...

I wanted to use the extract function to compute the mean temperature per plot. I first added the RD New projection to the raster and the shapefile. Afterwards, I tried to run the extract function and add the result to leaflet within Shiny.

#Set working directory right
setwd("...")

#Libraries
library(shiny)
library(ggplot2)/
library(leaflet)
library(rgdal)
library(raster)
library(rgeos)
library(plyr)
library(RColorBrewer)

#Get tiffile soil temperature
soil_temp_raster<- raster('Predicted_temperature_raster/C00_Soiltemp_predicted.tif')

#Project raster to RD New
crs(soil_temp_raster) <- "+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.2369,50.0087,465.658,-0.406857330322398,0.350732676542563,-1.8703473836068,4.0812 +units=m +no_defs"

test_raster<- projectRaster(soil_temp_raster, crs = CRS("+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.2369,50.0087,465.658,-0.406857330322398,0.350732676542563,-1.8703473836068,4.0812 +units=m +no_defs"))
test_raster_lowres<- aggregate(test_raster, fac=5)

#Import brp2017 and reproject it to RD New
brp2017polygons<- readOGR(dsn = "BRP2017", layer = "A01_BRP2017_wgs84")
test_shapefile<- spTransform(brp2017polygons, CRSobj = CRS("+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.2369,50.0087,465.658,-0.406857330322398,0.350732676542563,-1.8703473836068,4.0812 +units=m +no_defs"))

#Combine raster and shapefile into extract function
extract_test<- extract(x=test_raster_lowres, y=test_shapefile, fun=mean, na.rm=T, df = T, sp = T)

##Define UI ----
ui <- fluidPage(

#Outputs
mainPanel(
leafletOutput(outputId = "mymap", height= 600, width = '100%'))
)

server <- function(input, output, session) {

output$mymap <- renderLeaflet({
leaflet(extract_test) %>% addTiles() %>%
addPolygons(weight = 1, smoothFactor = 0.5,
opacity = 1.0, fillOpacity = 0.5, popup=~values)
})
}

Create a Shiny app object
shinyApp(ui = ui, server = server)


The final result looks like this:
Shiny application

Something has gone wrong with projection I guess. Could someone help me out?

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.