The objective of the AGROMET project is to provide hourly 1km² gridded datasets of weather parameters with the best accuracy (i.e. spatialize hourly records from the stations on the whole area of Wallonia). Weather parameters will be predicted with explanatory variables. These explanatory variables are :

But other explanatory variables are discussed to be included in models. For example, CORINE land cover could be a great explanatory variables. Indeed, urban areas or agricultural areas could have an influence on weather parameters.

To include CORINE land cover data in models, we need to find them on the Internet and then, handle them in R.

Get CORINE land cover data on the Internet

CORINE land cover data are available for Europe, they are compile in large files with raster or vector format. They can be found on Copernicus (the European Union’s Earth Observation Programme). Their vector files are very large because they are databases.

However, we only need data from Wallonia. The research can be done on CORINE land cover data for Belgium. The Geoportal of the federal state of Belgium is a great source of informations and CORINE land cover data are available.

The CORINE land cover for 2012 is available as CLC_BE_shp_L08.zip with the shapefile format.

Then, you can use it in R.

Preparing the integration of CORINE land cover data in models

Reading shapefile in R

R can read shapefile thanks to maptools package which depends on sp package.

# Download shapefile
# Warning : go to https://ac.ngi.be/client-open/vndl5zdmH4wgScjyTVdS?language=en&openpath=ngi-standard-open%2FVectordata%2FCLC2012%2FCLC12_BE_shp_L08.zip&tab=dataaccess&auth=false&open=true&accesscode=vndl5zdmH4wgScjyTVdS
# and accept the conditions and then, download the file CLC12_BE.zip and unzip it in a folder
library(maptools)
corine.sp <- readShapePoly("../data/CLC12_BE.shp")

The shapefile will be loaded as a Large SpatialPolygonsDataFrame.

Clipping the shapefile to Wallonia limits

The data frame has data from all Belgium but we only want data from Wallonia. As a consequence, we need to clip the file.

To do that, we get a file with limits of Wallonia and we use it to delete data beyond the limits. Then, you have a lighter SpatialPolygonsDataFrame.

library(raster)
corine.wal.sp <- crop(corine.sp, wallonie.3812.sp)

Loading the legend and analyzing it

The shapefile has a legend with different codes. There are 47 labels corresponding to different type of land covering like urban areas or agricultural areas. Only 26 of them are in Wallonia.

This legend is available here.

# Download legend
download.file("http://www.eea.europa.eu/data-and-maps/data/corine-land-cover-2006-raster-1/corine-land-cover-classes-and/clc_legend.csv/at_download/file",
              destfile = "../data/clc_legend.csv")
legend <- read.csv(file = "../data/clc_legend.csv", header = TRUE, sep = ";")
# Legend codes for Wallonia
legend.code.wal <- data.frame(unique(corine.wal.sp$code_12))
# Legend for CORINE land cover in Wallonia
legend.wal <- subset(legend, CLC_CODE %in% legend.code.wal$unique.corine.wal.sp.code_12.)

Now we have the different land covers. We can group some of them if they are the same influence on weather parameters.

We decided to create 5 groups :

  • Artificial areas : it concerns fabrics, airports, urban areas, dump sites, extraction sites, roads, rails… where vegetation is insignificant and shade is prevalent.

  • Agricultural areas : it concerns crops, arable lands…

  • Herbaceous vegetation : pastures, natural grasslands, moors and heathlands

  • Forest : it concerns all the forests (coniferous, broad-leaved, mixed) and shrubs. They are a large source of shade that could influence weather parameters like temperature.

  • Water : it concerns water bodies and wetlands…

These groups corresponds to 5 influences that we estimated different. Our first tests will run with them but it is possible that we change our groupes later to adapt them to improve our models.

Identifying land covers near to each PAMESEB stations

Measures of PAMESEB stations could be influenced by land cover. That’s why we need to identify the environment of each station.

First, we get data from the position of every station in Lambert Belgium 2008 CRS.

stations.sp <- build_agromet_stations_points.sp.fun()
stations.sp <- spTransform(stations.sp, CRSobj = "+proj=lcc +lat_1=49.83333333333334 +lat_2=51.16666666666666 +lat_0=50.797815 +lon_0=4.359215833333333 +x_0=649328 +y_0=665262 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
stations.sf <- st_as_sf(stations.sp)

Then, we make a buffer around stations and we cross-reference data from CLC and stations. The distance of the buffer will be changed later to estimate the better distance for our models.

# Make a buffer around stations
# https://gis.stackexchange.com/questions/229453/create-a-circle-of-defined-radius-around-a-point-and-then-find-the-overlapping-a
# https://stackoverflow.com/questions/46704878/circle-around-a-geographic-point-with-st-buffer
stations.buff.sf <- st_buffer(x = stations.sf, dist = dist.num)
  
# Cross-reference data to find the different land covers in the buffer
st_crs(corine.wal.simple.sf) <- st_crs(stations.buff.sf)
class.buff.stations.sf <- st_intersection(corine.wal.simple.sf, stations.buff.sf)

Then, we calculate the prevalence of every land cover around each station.

This work gave me an exercise to fix a problem in my code, I refer you to this post.

##   sid                 CLASS common_area rate_cover
## 1   1  Artificials surfaces    1340.126   4.265755
## 2   1    Agricultural areas   20061.639  63.858179
## 3   1 Herbaceous vegetation    9999.809  31.830380
## 4   4  Artificials surfaces   11257.894  35.834990
## 5   4    Agricultural areas   20143.680  64.119324
## 6   7    Agricultural areas   23650.339  75.281367

Further experiments

Conversion of the data

The precedent table is raw and not clear. It can be transformed to be more readable.

##   sid Agricultural areas Herbaceous vegetation   Forest
## 1   1           63.85818              31.83038  0.00000
## 2  10           69.91902               0.00000 30.03530
## 3  13           89.50912               0.00000 10.44519
## 4  14           31.52348               0.00000  0.00000
## 5  15           99.95431               0.00000  0.00000
## 6  17           99.95431               0.00000  0.00000
##   Artificials surfaces
## 1             4.265755
## 2             0.000000
## 3             0.000000
## 4            68.430833
## 5             0.000000
## 6             0.000000

Extending the function to entire Wallonia

The objective of the project is to create a map with 16894 points, one for every km² in Wallonia. This work on CORINE land cover can be extended to these points taking a grid with all these points in the function previously created.

Obviously, the computation is longer because we have 16894 points whereas the first computation (on real stations) has only 29 points. It tooks approximately 7 minutes but it works perfectly.

Below, a result for 163 points (mapview cannot print 16894 points on the map). 

##       id Agricultural areas Herbaceous vegetation   Forest
## 1   vs_1           47.95830              34.02409 17.97192
## 2  vs_10            0.00000               0.00000 99.95431
## 3 vs_100           87.38790               0.00000  0.00000
## 4 vs_101           23.97087               0.00000  0.00000
## 5 vs_102           72.47716               0.00000  0.00000
## 6 vs_103           95.18502               0.00000  0.00000
##   Artificials surfaces
## 1             0.000000
## 2             0.000000
## 3            12.566414
## 4            75.983449
## 5            27.477157
## 6             4.769293