1

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: Spain map with interpolation layer of income 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)
1

Got it! The lightbulb sparked up this morning while reading this answer.

All I have to do is clip the raster using the mask() function with the raster and the polygon as arguments:

rr <- mask(residualgrid, FRA)

New_map

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.