I was going through this very interesting code in order to make a fancy interactive map of Spain and a couple of questions popped into my mind:
How could I add the surface layer so that instead of a rectangle defined by grd it is mapped to the shape of Spain? Or if this is impossible (maybe I don't understand how the Kriging modelization is done), is there any way to "cut" out everything that exceeds the borders? I thought about changing the grd so that the newdata parameter from the idw function remained inside the Spanish borders, but no success so far, it get errors when trying to convert to pixels using gridded(grd)
Is there any way to add a button to the leaflet map to activate/deactivate the layers?
I have been doing research to find out, but I am a newbie regarding GIS and right now I'm a bit lost... any help please?
Photo of my actual results so far: Sample code from the related question:
library(ggplot2)
library(gstat)
library(sp)
library(maptools)
require(rgeos)
library(rgdal)
library(foreign)
library(maptools)
library(knitr)
library(jsonlite)
library(raster)
library(leaflet)
#récupérer les données a partir de json
dataJson <- fromJSON(readLines('jsonallcoor.js'))
# Téléchargement du fond de carte
FRA=readShapePoly("FRA_adm0.shp")
#concour (interpolation idw)
x=data.frame(0)
y=data.frame(0)
t=data.frame(0)
for ( i in 1:103) {
x[i]<- as.data.frame.factor(dataJson[[i]]$longitude)# define x & y as longitude and latitude
y[i]<- as.data.frame.factor(dataJson[[i]]$latitude)
t[i]<- as.data.frame.factor(dataJson[[i]]$temperature)
}
#changer le nom de dataframe
x<- as.numeric(x)
y <- as.numeric(y)
temperature<- as.numeric(t)
frame=as.data.frame(cbind(x,y,temperature))
frame.xy = frame[c("x", "y")]
coordinates(frame.xy ) <- ~x+y
plot(frame.xy )
#Define the grid extent
x.range <- as.numeric(c(-4.445833,9.484722)) # min/max longitude of the interpolation area
y.range <- as.numeric(c(40.50306,51.2)) # min/max latitude of the interpolation area
#Create a data frame from all combinations of the supplied vectors or factors. See the description of the return value for precise details of the way this is done. Set spatial coordinates to create a Spatial object. Assign gridded structure:
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.1),
y = seq(from = y.range[1], to = y.range[2], by = 0.1))
# expand points to grid
coordinates(grd) <- ~x + y
#to pixel
gridded(grd) <- TRUE
#Plot the weather station locations and interpolation grid:
plot(grd, cex = 1.5, col = "grey")
#les lon et lat
points(frame.xy, pch =1, col = "black", cex = 0.1)
#Interpolate surface and fix the output
idw <- idw(formula = temperature ~ 1, locations = frame.xy,
newdata = grd) # apply idw model for the data
residual_grid = raster(idw, "var1.pred")
contourr <- rasterToContour(residual_grid)
library(leaflet)
## Initialisation
m <- leaflet(padding = 0)
m <- addTiles(m)## Ajout des pays-
m <- addPolygons(map = m, data = FRA, opacity = 100,
color = "#FAFCFA",
weight = 0.25,popup = NULL,
options = list(clickable = FALSE),
fill = T, fillColor = "#B3C4B3",
fillOpacity = 100)
## Ajout des cercles
for ( i in 1:103) {
m <-addCircleMarkers(map = m,
lng = dataJson[[i]]$longitude,
lat = dataJson[[i]]$latitude,
radius=5, weight = 0.25,
stroke = T, opacity = 100,
fill = T, fillColor = "#920000",
fillOpacity = 100,
popup = dataJson[[i]]$station,
color = "white")
}
m <- fitBounds(map = m,
lng1 = -4.445833,
lat1 = 40.50306,
lng2 = 9.484722,
lat2 = 50.56417)
m<-addPolylines(map=m,data =contourr,fillOpacity=5,fillColor = "transparent",opacity=10,weight=1)
m<-addRectangles(map=m,
lng1=9.484722, lat1=40.50306,
lng2=-4.445833, lat2=51.2,
fillColor = "transparent"
)
## Dimensions de la carte
m$width <- 1000
m$height <-700
library(htmlwidgets)
saveWidget(m, 'contour.html', selfcontained = FALSE)