In paul-carteron/happign: R Interface to 'IGN' Web Services
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.align = "center",
fig.height = 4,
fig.width = 6
)
library(happign)
library(sf)
library(dplyr)
library(ggplot2)
library(terra)
library(tmap)
tmap_mode("view")
First choose a zone of interest
For the example we will work with the Camors forest. First of all we need the commune border that can be obtained with get_apicarto_commune from the insee code. For convenience, happign provides a table containing all insee codes (data("com_2025")).
data("com_2025")
#search communes name Camors and extract insee code
insee_code <- com_2025[grepl("^Camors", com_2025$LIBELLE_COM),"COM"]
camors_border <- get_apicarto_cadastre(insee_code, type = "commune")
tm_shape(camors_border)+
tm_borders()+
tm_text("nom_com")
An other way of getting borders without Apicarto is to use ECQL language to directly query IGN WFS geoservers.
camors_border2 <- get_wfs(
layer = "ADMINEXPRESS-COG-CARTO.LATEST:commune",
ecql_filter = "nom_m LIKE 'CAMORS'"
)
# "nom_m" attribut can be found using get_wfs_attributes(layer = "ADMINEXPRESS-COG-CARTO.LATEST:commune")
Download cadastral parcel
Cadastral parcels are essential for any forest manager. Here how to download it with get_wfs.
layers <- get_layers_metadata("wfs", "parcellaire")
parcellaire_layer <- layers[15,1] # "CADASTRALPARCELS.PARCELLAIRE_EXPRESS:parcelle"
camors_parcels <- get_wfs(
x = camors_border,
spatial_filter = "intersects",
layer = parcellaire_layer)
tm_shape(camors_border)+
tm_borders(col = "red", lwd = 2)+
tm_shape(camors_parcels)+
tm_polygons(fill_alpha = 0)
Rq :
- IGN WFS can return a maximum of 1000 features. The
get_wfs() function overrides this limit by performing several consecutive requests as the console show
Downloading BD Forêt
The first interesting layer for forester is the "BD Forêt" which is all vegetation type assigned to each area greater than or equal to 0.5 ha (5,000 m²). There is two layer for forest : the old one BD Forêt V1 and the new one BD Forêt V2 that can be accessed with "burger menu" on the top left of interactive map below.
layers <- get_layers_metadata("wfs", "environnement")
BDF1_layer <- layers[4,1]
BDF2_layer <- layers[5,1]
BDF1 <- get_wfs(camors_border, BDF1_layer, spatial_filter = "intersects")
BDF2 <- get_wfs(camors_border, BDF2_layer, spatial_filter = "intersects")
# tmap V4 code
tm_shape(BDF1) +
tm_polygons(fill = "libelle",
popup.vars = names(BDF1)[1:(ncol(BDF1)-2)],
fill.legend = tm_legend_hide())+
tm_shape(BDF2) +
tm_polygons(fill = "tfv",
fill_alpha = 0.5,
popup.vars = names(BDF2)[1:(ncol(BDF2)-2)],
fill.legend = tm_legend_hide()) +
tm_shape(camors_border) +
tm_borders(lwd = 2)
# # tmap v3 code
# tm_shape(BDF1) +
# tm_polygons(col = "libelle",
# popup.vars = names(BDF1)[1:(ncol(BDF1)-2)],
# legend.show = FALSE)+
# tm_shape(BDF2) +
# tm_polygons(col = "tfv",
# alpha = 0.5,
# popup.vars = names(BDF2)[1:(ncol(BDF2)-2)],
# legend.show = FALSE) +
# tm_shape(camors_border) +
# tm_borders(lwd = 2)
More calculations can be done to get a better understanding of forest type. Below, areas per species are calculated :
forest_type_BDF2 <- BDF2 |>
mutate(area = as.numeric(st_area(geometry))) |>
st_drop_geometry() |>
group_by(essence) |>
summarise(sum_area = sum(area)/10000) |>
arrange(desc(sum_area)) |>
mutate(essence = as.factor(essence))
ggplot()+
geom_col(data = forest_type_BDF2,
aes(x = rev(reorder(essence, sum_area)),
y = sum_area,
fill = as.factor(essence)))+
theme_bw()+
labs(title = "Surface couverte par essences [ha]",
y = "Surface [ha]",
x = "",
fill = "Essence :")+
theme(axis.text.x = element_blank())
Detect protected area
One information you really want when you're a forest manager is if your zone of interest is inside protected area. The example code below is design to automatically test every layer containing "PROTECTED" so you can be sure that you have all of them.
Again, you can click on map, point and shape for more informations.
protected_area_layers <- get_layers_metadata("wfs", "environnement") |>
filter(grepl(".*PROTECTED.*", Name)) |>
pull(Name)
get_wfs_on_camors <- function(protected_area_layer){
message(protected_area_layer,"...")
layer <- try(
get_wfs(
camors_border,
protected_area_layer,
spatial_filter = "intersects"
)
)
return(layer)
}
try_all_protected_area <- lapply(protected_area_layers, get_wfs_on_camors) |>
setNames(protected_area_layers)
# remove empty layers
all_protected_area <- Filter(function(x) !is.null(x), try_all_protected_area)
# Plot the result
tm_shape(camors_border) +
tm_borders(lwd = 2)+
tm_shape(all_protected_area[[1]])+
tm_polygons(fill_alpha = 0.8, col = "blue")
MNS, MNT and MNH...
Terrain topologie is quite useful for forest management. IGN offers MNT and MNS for download. As a reminder, the MNT corresponds to the surface of the ground and the MNS to the real surface (in our case, the trees). It is thus easy to find the height of the trees by subtracting the MNT from the MNS.
layers <- get_layers_metadata("wms-r", "altimetrie")
mnt_layer <- layers[3,1] # "ELEVATION.ELEVATIONGRIDCOVERAGE.HIGHRES"
mns_layer <- layers[4,1] # "ELEVATION.ELEVATIONGRIDCOVERAGE.HIGHRES.MNS"
mnt <- get_wms_raster(camors_border, mnt_layer, res = 15, crs = 2154, rgb = FALSE)
mns <- get_wms_raster(camors_border, mns_layer, res = 15, crs = 2154, rgb = FALSE)
# Add level curve for flexing
level_curve <- get_wfs(camors_border, "ELEVATION.CONTOUR.LINE:courbe",
spatial_filter = "intersects") |>
st_intersection(camors_border)
# Calculate digital height model i.e. tree height
mnh <- mns - mnt
mnh[mnh < 0] <- NA # Remove negative value
mnh[mnh > 50] <- 40 # Remove height more than 50m
# tmap v4 code
tm_shape(mnh) +
tm_raster(col.scale = tm_scale_continuous(
values = "carto.earth", value.na = "grey"
),
col.legend = tm_legend(title = "Height", showNA = F)) +
tm_shape(level_curve)+
tm_lines(col="grey25", col_alpha = 0.8)+
tm_shape(camors_border)+
tm_borders(lwd = 2, col = "red")
# tmap v3
# tm_shape(mnh) +
# tm_raster(style = "cont",
# title = "Height",
# palette = "-Spectral",
# colorNA = "grey",
# showNA = F) +
# tm_shape(level_curve)+
# tm_lines(col = "black")+
# tm_shape(camors_border)+
# tm_borders(lwd = 2, col = "red")
NDVI
The code below present the calculation of the NDVI. All informations and palette come from this website
The value range of an NDVI is -1 to 1. It is (Near Infrared - Red) / (Near Infrared + Red) :
- Water has a low reflectance in red, but almost no NIR (near infrared) reflectance. So the difference will be small and negative, and the sum will be small, and NDVI large and negative.
- Plants have a low reflectance in red, and a strong NIR reflectance. So the difference will be large and positive, and the sum will be just about the same as the difference, so NDVI will be large and positive.
Categories are somewhat arbitrary, and you can find various rules of thumb, such as:
- Negative values of NDVI (values approaching -1) correspond to water. Values close to zero (-0.1 to 0.1) generally correspond to barren areas of rock, sand, or snow. Low, positive values represent shrub and grassland (approximately 0.2 to 0.4), while high values indicate temperate and tropical rainforests (values approaching 1).
- Very low values of NDVI (0.1 and below) correspond to water, barren areas of rock, sand, or snow. Moderate values represent shrub and grassland (0.2 to 0.3), while high values indicate temperate and tropical rainforests (0.6 to 0.8).
# To show the 80cm resolution possibility for IRC, let's take only the biggest parcels
biggest_parcels <- camors_parcels |>
mutate(area = st_area(geometry)) |>
slice_max(area)
irc <- get_wms_raster(biggest_parcels,
res = 0.8,
layer = "ORTHOIMAGERY.ORTHOPHOTOS.IRC")
# calculate ndvi from near_infrared and infrared
ndvi_fun <- function(nir, red){
(nir - red) / (nir + red)
}
ndvi <- lapp(irc[[c(1, 2)]], fun = ndvi_fun)
# palette for plotting
breaks_ndvi <- c(-1,-0.2,-0.1,0,0.025 ,0.05,0.075,0.1,0.125,0.15,0.175,0.2 ,0.25 ,0.3 ,0.35,0.4,0.45,0.5,0.55,0.6,1)
palette_ndvi <- c("#BFBFBF","#DBDBDB","#FFFFE0","#FFFACC","#EDE8B5","#DED99C","#CCC782","#BDB86B","#B0C261","#A3CC59","#91BF52","#80B347","#70A340","#619636","#4F8A2E","#407D24","#306E1C","#216112","#0F540A","#004500")
#tmap v4 code
tm_shape(biggest_parcels, is.master = TRUE)+
tm_borders(lwd = 2, col = "blue")+
tm_shape(ndvi)+
tm_raster(
col.scale = tm_scale_continuous(
limits = c(-1, 1),
values = palette_ndvi,
),
col.legend = tm_legend(title = "NDVI"))
# # tmap v3 code
# tm_shape(camors_border)+
# tm_borders(lwd = 2, col = "red")+
# tm_shape(ndvi)+
# tm_raster(stretch.palette = F,
# style = "cont",
# title = "NDVI",
# breaks = breaks_ndvi,
# palette = palette_ndvi,
# colorNA = NULL)+
# tm_shape(biggest_parcels, is.master = TRUE)+
# tm_borders(lwd = 2, col = "blue")
The gloss index
The gloss index represents the average of the image glosses. This index is therefore sensitive to the brightness of the soil, related to its moisture and the presence of salts on the surface. It characterizes especially the albedo (solar radiation that is reflected back to the atmosphere). The gloss index allows us to estimate whether the observed surface feature is light or dark.
# calculate gloss_index from near_infrared and infrared
gloss_fun <- function(nir, red){
sqrt(red^2 + nir^2)
}
gloss_index <- lapp(irc[[c(1, 2)]],
fun = gloss_fun)
# tmap v4 code
tm_shape(biggest_parcels)+
tm_borders(lwd = 2, col = "blue")+
tm_shape(gloss_index)+
tm_raster(col.legend = tm_legend(title = "GLOSS INDEX"))
# # tmap v3
# tm_shape(camors_border)+
# tm_borders(lwd = 2, col = "red")+
# tm_shape(gloss_index)+
# tm_raster(style = "cont",
# title = "GLOSS INDEX")+
# tm_shape(biggest_parcels, is.master = T)+
# tm_borders(lwd = 2, col = "blue")
Info on raster
For some resources, it is essential to have additional information in order to use them properly. For example, in the forest field, the date the picture was taken is important because our object of study is dynamic. The function get_location_info() allows to get, if they exist, additional information on the raster layers. Warning : only point are supported.
Rq :
- Note that the function returns a shape by default. If only the metadata is required, you can speed up the response time by setting read_sf to FALSE.
- To find all queryable layers from an apikey you can use
are_queryable(apikey)
# info_sup <- get_location_info(x = st_centroid(camors_border),
# apikey = "ortho",
# layer = "ORTHOIMAGERY.ORTHOPHOTOS.BDORTHO",
# read_sf = F)
# info_sup$date_vol
The call to the function get_location_info indicates that this orthophoto was taken the 2019-05-14.
In some cases, the function can be used to retrieve resources that are not explicitly available. For example, there are no vector layers of public forests, only a raster layer. However, when additional information is requested, the outline of the forest is returned.
# x <- st_sfc(st_point(c(-3.549957, 47.83396)), crs = 4326) # Carnoet forest
#
# forest <- get_location_info(x,
# apikey = "environnement",
# layer = "FORETS.PUBLIQUES",
# read_sf = TRUE)
# tm_shape(forest)+
# tm_borders()+
# tm_shape(x)+
# tm_dots(size = 0.05, col = "red")
Last but not least... BD Topo
BD topo from IGN covers in a coherent way all the geographical and administrative entities of the national territory. So you can find in it :
- Administrative (boundaries and administrative units);
- Addresses (mailing addresses) ;
- Building (constructions) ;
- Hydrography (water-related features) ;
- Named places (place or locality with a toponym describing a natural space or inhabited place);
- Land use (vegetation, foreshore, hedge);
- Services and activities (utilities, energy storage and transportation, industrial sites);
- Transportation (road, rail or air infrastructure, routes);
- Regulated areas (most of the areas are subject to specific regulations).
For the example below I choose to download all water-related data :
cour_eau <- get_wfs(camors_border, "BDTOPO_V3:cours_d_eau") |>
st_intersection(camors_border)
detail_hydro <- get_wfs(camors_border, "BDTOPO_V3:detail_hydrographique") |>
st_intersection(camors_border)
# water detected by satellite
surf_hydro <- get_wfs(camors_border, "BDTOPO_V3:surface_hydrographique") |>
st_intersection(camors_border)
tm_shape(cour_eau)+
tm_lines(col = "blue")+
tm_shape(detail_hydro)+
tm_dots(col = "red")+
tm_shape(surf_hydro)+
tm_polygons("steelblue")+
tm_shape(camors_border)+
tm_borders(lwd = 2)
paul-carteron/happign documentation built on June 15, 2025, 8:24 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.align = "center", fig.height = 4, fig.width = 6 )
library(happign) library(sf) library(dplyr) library(ggplot2) library(terra) library(tmap) tmap_mode("view")
First choose a zone of interest
For the example we will work with the Camors forest. First of all we need the commune border that can be obtained with get_apicarto_commune from the insee code. For convenience, happign provides a table containing all insee codes (data("com_2025")).
data("com_2025") #search communes name Camors and extract insee code insee_code <- com_2025[grepl("^Camors", com_2025$LIBELLE_COM),"COM"] camors_border <- get_apicarto_cadastre(insee_code, type = "commune") tm_shape(camors_border)+ tm_borders()+ tm_text("nom_com")
An other way of getting borders without Apicarto is to use ECQL language to directly query IGN WFS geoservers.
camors_border2 <- get_wfs( layer = "ADMINEXPRESS-COG-CARTO.LATEST:commune", ecql_filter = "nom_m LIKE 'CAMORS'" ) # "nom_m" attribut can be found using get_wfs_attributes(layer = "ADMINEXPRESS-COG-CARTO.LATEST:commune")
Download cadastral parcel
Cadastral parcels are essential for any forest manager. Here how to download it with get_wfs.
layers <- get_layers_metadata("wfs", "parcellaire") parcellaire_layer <- layers[15,1] # "CADASTRALPARCELS.PARCELLAIRE_EXPRESS:parcelle" camors_parcels <- get_wfs( x = camors_border, spatial_filter = "intersects", layer = parcellaire_layer) tm_shape(camors_border)+ tm_borders(col = "red", lwd = 2)+ tm_shape(camors_parcels)+ tm_polygons(fill_alpha = 0)
Rq :
- IGN WFS can return a maximum of 1000 features. The
get_wfs()function overrides this limit by performing several consecutive requests as the console show
Downloading BD Forêt
The first interesting layer for forester is the "BD Forêt" which is all vegetation type assigned to each area greater than or equal to 0.5 ha (5,000 m²). There is two layer for forest : the old one BD Forêt V1 and the new one BD Forêt V2 that can be accessed with "burger menu" on the top left of interactive map below.
layers <- get_layers_metadata("wfs", "environnement") BDF1_layer <- layers[4,1] BDF2_layer <- layers[5,1] BDF1 <- get_wfs(camors_border, BDF1_layer, spatial_filter = "intersects") BDF2 <- get_wfs(camors_border, BDF2_layer, spatial_filter = "intersects") # tmap V4 code tm_shape(BDF1) + tm_polygons(fill = "libelle", popup.vars = names(BDF1)[1:(ncol(BDF1)-2)], fill.legend = tm_legend_hide())+ tm_shape(BDF2) + tm_polygons(fill = "tfv", fill_alpha = 0.5, popup.vars = names(BDF2)[1:(ncol(BDF2)-2)], fill.legend = tm_legend_hide()) + tm_shape(camors_border) + tm_borders(lwd = 2)
# # tmap v3 code # tm_shape(BDF1) + # tm_polygons(col = "libelle", # popup.vars = names(BDF1)[1:(ncol(BDF1)-2)], # legend.show = FALSE)+ # tm_shape(BDF2) + # tm_polygons(col = "tfv", # alpha = 0.5, # popup.vars = names(BDF2)[1:(ncol(BDF2)-2)], # legend.show = FALSE) + # tm_shape(camors_border) + # tm_borders(lwd = 2)
More calculations can be done to get a better understanding of forest type. Below, areas per species are calculated :
forest_type_BDF2 <- BDF2 |> mutate(area = as.numeric(st_area(geometry))) |> st_drop_geometry() |> group_by(essence) |> summarise(sum_area = sum(area)/10000) |> arrange(desc(sum_area)) |> mutate(essence = as.factor(essence)) ggplot()+ geom_col(data = forest_type_BDF2, aes(x = rev(reorder(essence, sum_area)), y = sum_area, fill = as.factor(essence)))+ theme_bw()+ labs(title = "Surface couverte par essences [ha]", y = "Surface [ha]", x = "", fill = "Essence :")+ theme(axis.text.x = element_blank())
Detect protected area
One information you really want when you're a forest manager is if your zone of interest is inside protected area. The example code below is design to automatically test every layer containing "PROTECTED" so you can be sure that you have all of them.
Again, you can click on map, point and shape for more informations.
protected_area_layers <- get_layers_metadata("wfs", "environnement") |> filter(grepl(".*PROTECTED.*", Name)) |> pull(Name) get_wfs_on_camors <- function(protected_area_layer){ message(protected_area_layer,"...") layer <- try( get_wfs( camors_border, protected_area_layer, spatial_filter = "intersects" ) ) return(layer) } try_all_protected_area <- lapply(protected_area_layers, get_wfs_on_camors) |> setNames(protected_area_layers) # remove empty layers all_protected_area <- Filter(function(x) !is.null(x), try_all_protected_area) # Plot the result tm_shape(camors_border) + tm_borders(lwd = 2)+ tm_shape(all_protected_area[[1]])+ tm_polygons(fill_alpha = 0.8, col = "blue")
MNS, MNT and MNH...
Terrain topologie is quite useful for forest management. IGN offers MNT and MNS for download. As a reminder, the MNT corresponds to the surface of the ground and the MNS to the real surface (in our case, the trees). It is thus easy to find the height of the trees by subtracting the MNT from the MNS.
layers <- get_layers_metadata("wms-r", "altimetrie") mnt_layer <- layers[3,1] # "ELEVATION.ELEVATIONGRIDCOVERAGE.HIGHRES" mns_layer <- layers[4,1] # "ELEVATION.ELEVATIONGRIDCOVERAGE.HIGHRES.MNS" mnt <- get_wms_raster(camors_border, mnt_layer, res = 15, crs = 2154, rgb = FALSE) mns <- get_wms_raster(camors_border, mns_layer, res = 15, crs = 2154, rgb = FALSE) # Add level curve for flexing level_curve <- get_wfs(camors_border, "ELEVATION.CONTOUR.LINE:courbe", spatial_filter = "intersects") |> st_intersection(camors_border) # Calculate digital height model i.e. tree height mnh <- mns - mnt mnh[mnh < 0] <- NA # Remove negative value mnh[mnh > 50] <- 40 # Remove height more than 50m # tmap v4 code tm_shape(mnh) + tm_raster(col.scale = tm_scale_continuous( values = "carto.earth", value.na = "grey" ), col.legend = tm_legend(title = "Height", showNA = F)) + tm_shape(level_curve)+ tm_lines(col="grey25", col_alpha = 0.8)+ tm_shape(camors_border)+ tm_borders(lwd = 2, col = "red")
# tmap v3 # tm_shape(mnh) + # tm_raster(style = "cont", # title = "Height", # palette = "-Spectral", # colorNA = "grey", # showNA = F) + # tm_shape(level_curve)+ # tm_lines(col = "black")+ # tm_shape(camors_border)+ # tm_borders(lwd = 2, col = "red")
NDVI
The code below present the calculation of the NDVI. All informations and palette come from this website The value range of an NDVI is -1 to 1. It is (Near Infrared - Red) / (Near Infrared + Red) :
- Water has a low reflectance in red, but almost no NIR (near infrared) reflectance. So the difference will be small and negative, and the sum will be small, and NDVI large and negative.
- Plants have a low reflectance in red, and a strong NIR reflectance. So the difference will be large and positive, and the sum will be just about the same as the difference, so NDVI will be large and positive.
Categories are somewhat arbitrary, and you can find various rules of thumb, such as:
- Negative values of NDVI (values approaching -1) correspond to water. Values close to zero (-0.1 to 0.1) generally correspond to barren areas of rock, sand, or snow. Low, positive values represent shrub and grassland (approximately 0.2 to 0.4), while high values indicate temperate and tropical rainforests (values approaching 1).
- Very low values of NDVI (0.1 and below) correspond to water, barren areas of rock, sand, or snow. Moderate values represent shrub and grassland (0.2 to 0.3), while high values indicate temperate and tropical rainforests (0.6 to 0.8).
# To show the 80cm resolution possibility for IRC, let's take only the biggest parcels biggest_parcels <- camors_parcels |> mutate(area = st_area(geometry)) |> slice_max(area) irc <- get_wms_raster(biggest_parcels, res = 0.8, layer = "ORTHOIMAGERY.ORTHOPHOTOS.IRC") # calculate ndvi from near_infrared and infrared ndvi_fun <- function(nir, red){ (nir - red) / (nir + red) } ndvi <- lapp(irc[[c(1, 2)]], fun = ndvi_fun) # palette for plotting breaks_ndvi <- c(-1,-0.2,-0.1,0,0.025 ,0.05,0.075,0.1,0.125,0.15,0.175,0.2 ,0.25 ,0.3 ,0.35,0.4,0.45,0.5,0.55,0.6,1) palette_ndvi <- c("#BFBFBF","#DBDBDB","#FFFFE0","#FFFACC","#EDE8B5","#DED99C","#CCC782","#BDB86B","#B0C261","#A3CC59","#91BF52","#80B347","#70A340","#619636","#4F8A2E","#407D24","#306E1C","#216112","#0F540A","#004500") #tmap v4 code tm_shape(biggest_parcels, is.master = TRUE)+ tm_borders(lwd = 2, col = "blue")+ tm_shape(ndvi)+ tm_raster( col.scale = tm_scale_continuous( limits = c(-1, 1), values = palette_ndvi, ), col.legend = tm_legend(title = "NDVI"))
# # tmap v3 code # tm_shape(camors_border)+ # tm_borders(lwd = 2, col = "red")+ # tm_shape(ndvi)+ # tm_raster(stretch.palette = F, # style = "cont", # title = "NDVI", # breaks = breaks_ndvi, # palette = palette_ndvi, # colorNA = NULL)+ # tm_shape(biggest_parcels, is.master = TRUE)+ # tm_borders(lwd = 2, col = "blue")
The gloss index
The gloss index represents the average of the image glosses. This index is therefore sensitive to the brightness of the soil, related to its moisture and the presence of salts on the surface. It characterizes especially the albedo (solar radiation that is reflected back to the atmosphere). The gloss index allows us to estimate whether the observed surface feature is light or dark.
# calculate gloss_index from near_infrared and infrared gloss_fun <- function(nir, red){ sqrt(red^2 + nir^2) } gloss_index <- lapp(irc[[c(1, 2)]], fun = gloss_fun) # tmap v4 code tm_shape(biggest_parcels)+ tm_borders(lwd = 2, col = "blue")+ tm_shape(gloss_index)+ tm_raster(col.legend = tm_legend(title = "GLOSS INDEX"))
# # tmap v3 # tm_shape(camors_border)+ # tm_borders(lwd = 2, col = "red")+ # tm_shape(gloss_index)+ # tm_raster(style = "cont", # title = "GLOSS INDEX")+ # tm_shape(biggest_parcels, is.master = T)+ # tm_borders(lwd = 2, col = "blue")
Info on raster
For some resources, it is essential to have additional information in order to use them properly. For example, in the forest field, the date the picture was taken is important because our object of study is dynamic. The function get_location_info() allows to get, if they exist, additional information on the raster layers. Warning : only point are supported.
Rq :
- Note that the function returns a shape by default. If only the metadata is required, you can speed up the response time by setting read_sf to FALSE.
- To find all queryable layers from an apikey you can use
are_queryable(apikey)
# info_sup <- get_location_info(x = st_centroid(camors_border), # apikey = "ortho", # layer = "ORTHOIMAGERY.ORTHOPHOTOS.BDORTHO", # read_sf = F) # info_sup$date_vol
The call to the function get_location_info indicates that this orthophoto was taken the 2019-05-14.
In some cases, the function can be used to retrieve resources that are not explicitly available. For example, there are no vector layers of public forests, only a raster layer. However, when additional information is requested, the outline of the forest is returned.
# x <- st_sfc(st_point(c(-3.549957, 47.83396)), crs = 4326) # Carnoet forest # # forest <- get_location_info(x, # apikey = "environnement", # layer = "FORETS.PUBLIQUES", # read_sf = TRUE) # tm_shape(forest)+ # tm_borders()+ # tm_shape(x)+ # tm_dots(size = 0.05, col = "red")
Last but not least... BD Topo
BD topo from IGN covers in a coherent way all the geographical and administrative entities of the national territory. So you can find in it :
- Administrative (boundaries and administrative units);
- Addresses (mailing addresses) ;
- Building (constructions) ;
- Hydrography (water-related features) ;
- Named places (place or locality with a toponym describing a natural space or inhabited place);
- Land use (vegetation, foreshore, hedge);
- Services and activities (utilities, energy storage and transportation, industrial sites);
- Transportation (road, rail or air infrastructure, routes);
- Regulated areas (most of the areas are subject to specific regulations).
For the example below I choose to download all water-related data :
cour_eau <- get_wfs(camors_border, "BDTOPO_V3:cours_d_eau") |> st_intersection(camors_border) detail_hydro <- get_wfs(camors_border, "BDTOPO_V3:detail_hydrographique") |> st_intersection(camors_border) # water detected by satellite surf_hydro <- get_wfs(camors_border, "BDTOPO_V3:surface_hydrographique") |> st_intersection(camors_border) tm_shape(cour_eau)+ tm_lines(col = "blue")+ tm_shape(detail_hydro)+ tm_dots(col = "red")+ tm_shape(surf_hydro)+ tm_polygons("steelblue")+ tm_shape(camors_border)+ tm_borders(lwd = 2)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.