library(happign)
library(sf)
library(tmap); tmap_mode("view") # Set map to interactive
library(dplyr)
library(ggplot2)
library(purrr)

First choose a zone of interest

For the example we will work with the Gâvre forest. The starting point come from google map : you just have to select a random point inside the forest by using right click on map.

point <- st_sfc(st_point(c(-1.7977020162748798, 47.52634370369097)), crs = 4326)

tm_shape(point)+
tm_dots(size = 0.1)+
tm_basemap("OpenStreetMap")

Extract commune borders

Now that we have point, getting shape of commune is just a matter of finding the good layer name from IGN. The one we are after is called “ADMINEXPRESS-COG-CARTO.LATEST:commune” from the “Administrative” category. Then we use the get_wfs() function to download the shape.

As you can see below, there are two communes close to our point (tip : click on the interactive map below to get more info). For following examples we will only use “Le Gâvre” shape.

apikey <- get_apikeys()[1]  # "administratif"

# The layer_name I use come from all_layers = get_layers_metadata(apikey, "wfs") that return all of them layer_name

borders <- get_wfs(shape = point,
apikey = apikey,
layer_name = layer_name)

tm_shape(borders)+
tm_polygons(lwd = 2,
col = "nom",
id = "nom",
popup.vars = "population")+
tm_basemap("OpenStreetMap")

# Select only commune to reduce downloading time for hight resolution raster
borders <- borders %>%
slice(1)

apikey <- get_apikeys()[14] # "parcellaire"

parcellaire <- get_wfs(shape = borders,
apikey = apikey,
layer_name = name_parcellaire_layer) %>%
st_transform(st_crs(borders)) %>%
st_intersection(borders)
tm_basemap("OpenStreetMap")`