Skip to contents

happign for foresters

Source: vignettes/web_only/happign_for_foresters.Rmd
happign_for_foresters.Rmd
library(happign)
library(sf)
library(tmap); tmap_mode("view") # Set map to interactive
library(dplyr)
library(ggplot2);sf_use_s2(FALSE) # Avoid problem with spherical geometry
library(purrr)
library(stars)
library(terra)

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 from the insee code. Fortunately, happign provides a table containing all insee codes (data("cog_2023")). Then, get_apicarto_commune is used to download shape from Apicarto commune.

data("cog_2023")

insee_code <- cog_2023[grepl("^Camors", cog_2023$LIBELLE),1]

borders <- get_apicarto_cadastre(insee_code, type = "commune")

tm_shape(borders)+
   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.

borders2 <- get_wfs(layer = "LIMITES_ADMINISTRATIVES_EXPRESS.LATEST:commune",
                    spatial_filter = "intersects",
                    ecql_filter = "nom_m LIKE 'CAMORS'")
#> Features downloaded : 1

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"

parcellaire <- get_wfs(x = borders,
                       spatial_filter = "intersects",
                       layer = parcellaire_layer)
#> Features downloaded : 1000...2000...3000...3490

tm_shape(borders)+
   tm_borders(col = "red", lwd = 2)+
tm_shape(parcellaire)+
   tm_polygons(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(borders, BDF1_layer, spatial_filter = "intersects")
#> Features downloaded : 102
BDF2 <- get_wfs(borders, BDF2_layer, spatial_filter = "intersects")
#> Features downloaded : 221

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(borders) +
   tm_borders(lwd = 2)

More calculations can be done as you can see below :

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]",
       fill = "Essence :")+
  theme(axis.text.x = element_blank())

Detect protected area

One information you really want when you work at forest management is if your zone of interest is inside protected area. The example code below is design to automatically test every layer starting with “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)

all_protected_area <- map(.x = protected_area_layers,
                          .f = ~ try(get_wfs(borders, .x, 
                                          spatial_filter = "intersects"))) |> 
   set_names(protected_area_layers) |> 
   discard(~ identical(length(.), 0L))

# Plot the result
tm_shape(all_protected_area[[1]])+
   tm_dots(group = "Point rencontre des secours en forêts", col = "red")+
tm_shape(all_protected_area[[2]])+
   tm_polygons(group = "Znieff 2", alpha = 0.8, col = "blue")+
tm_shape(borders,is.master = TRUE) +
   tm_borders(lwd = 2)

MNS, MNT and MNH…

It’s always good to know a more about terrain topologie. 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 DTM 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(borders, mnt_layer, res = 5, crs = 2154, rgb = FALSE)
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> Raster is saved at : /private/var/folders/0g/hj_q_pzx65bbjnslxz9n0src0000gn/T/Rtmp94yrpg/file368174a8314.tif
mns <- get_wms_raster(borders, mns_layer, res = 5, crs = 2154, rgb = FALSE)
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> Raster is saved at : /private/var/folders/0g/hj_q_pzx65bbjnslxz9n0src0000gn/T/Rtmp94yrpg/file36817708e4a2.tif

level_curve <- get_wfs(borders, "ELEVATION.CONTOUR.LINE:courbe",
                       spatial_filter = "intersects") |> 
   st_intersection(borders)
#> Features downloaded : 120
#> although coordinates are longitude/latitude, st_intersection assumes that they
#> are planar
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries

# 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

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(borders)+
   tm_borders(lwd = 2, col = "red")
#> stars object downsampled to 1027 by 974 cells. See tm_shape manual (argument raster.downsample)

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 20cm resolution possibility for IRC, let's take only the biggest parcels
biggest_parcels <- parcellaire |> 
   mutate(area = st_area(geometry)) |> 
   slice_max(area)

irc <- get_wms_raster(biggest_parcels,
                      res = 0.2,
                      layer = "ORTHOIMAGERY.ORTHOPHOTOS.IRC")
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> Raster is saved at : /private/var/folders/0g/hj_q_pzx65bbjnslxz9n0src0000gn/T/Rtmp94yrpg/file3681264a18ae.tif

# 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")

tm_shape(borders)+
   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")
#> stars object downsampled to 1082 by 924 cells. See tm_shape manual (argument raster.downsample)
#> Warning: Breaks contains positive and negative values. Better is to use
#> diverging scale instead, or set auto.palette.mapping to FALSE.

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)

tm_shape(borders)+
   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")
#> stars object downsampled to 1082 by 924 cells. See tm_shape manual (argument raster.downsample)

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(borders),
                              apikey = "ortho",
                              layer = "ORTHOIMAGERY.ORTHOPHOTOS.BDORTHO",
                              read_sf = F)
#> Warning: st_centroid assumes attributes are constant over geometries
#> Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
#> of_largest_polygon): st_centroid does not give correct centroids for
#> longitude/latitude data
info_sup$date_vol
#> [1] "2022-07-08Z"

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(borders, "BDTOPO_V3:cours_d_eau") |> 
   st_intersection(borders) 
#> Features downloaded : 45
#> although coordinates are longitude/latitude, st_intersection assumes that they
#> are planar
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
detail_hydro <- get_wfs(borders, "BDTOPO_V3:detail_hydrographique") |> 
   st_intersection(borders) 
#> Features downloaded : 25
#> although coordinates are longitude/latitude, st_intersection assumes that they
#> are planar
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries
# water detected by satellite
surf_hydro <- get_wfs(borders, "BDTOPO_V3:surface_hydrographique") |> 
   st_intersection(borders) 
#> Features downloaded : 89
#> although coordinates are longitude/latitude, st_intersection assumes that they
#> are planar
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries

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(borders)+
   tm_borders(lwd = 2)

What about history ?

The “Etat-major” map is a general map of France made, in its first version, in the 19th century. Here how to get it :

etat_major <- get_wms_raster(x = borders,
                             res = 1,
                             layer = "GEOGRAPHICALGRIDSYSTEMS.ETATMAJOR40")
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> Raster is saved at : /private/var/folders/0g/hj_q_pzx65bbjnslxz9n0src0000gn/T/Rtmp94yrpg/file36812d89a166.tif

# raster is very large and tmap cannot plot it. Convert to stars as below allow for quick plotting
# saving raster to disk allow full resolution visualisation : writeRaster(etat_major, "test.tif") 
etat_major <- st_as_stars(etat_major, ignore_file = TRUE)

tm_shape(etat_major)+
   tm_rgb()+
tm_shape(borders)+
   tm_borders(lwd = 3, col = "red")
#> stars object downsampled to 1028 by 974 cells. See tm_shape manual (argument raster.downsample)