### Sidebar

spatial-ecology.org

Trainings:

Learn:

Data:

Community:
teachers
students
projects
Matera 2015
Vancouver 2015
Santa Barbara 2015
Site design
Hands on:
Installations

Donations
USD

GBP

wikistudholland:utwente_aart

# Modelling Wheat Harvest Dates in France

Aart van der Linden
Animal Production Systems group,
Plant Production Systems group,
Wageningen University

## INTRODUCTION

Temperature and daylength are the most important drivers for crop development. The rate of development of cereals can be described as a function of the growing degree days (GDD). Wheat requires 1600 GDD from sowing to maturity and harvest. The amount of GDD increases when temperature is higher than a base temperature (BT). The BT is 4.5 degrees Celsius for wheat.

If the amount of GDD is 300, and today's average temperature is 16.5 degrees Celsius, the GDD at the end of the day equals 300 + 16.5 - 4.5 = 312 GDD. The amount of GDD cannot decrease due to temperatures lower than the BT. Wheat is consired to be ready for harvest when the amount of GDD exceeds 1600.

Differences in harvest date of wheat are expected between regions due to differences in climate. The aim of this excercise is to assess when farmers can harvest their wheat in France, if the wheat variety has a 1600 GDD requirement and 4.5 degrees Celsius BT.

### General framework of the analysis

R is used as program language. The following steps have been taken to assess the harvest date of wheat in France:

1. Collection of daily average temperature data for 0.5 x 0.5 degrees grid cells for 1995 till 2013(http://www.ecad.eu/download/ensembles/download.php). The .gz data are unzipped to .nc data files. The package ncdf in R can read these data files. The datafile covers Europe, West-Asia and North-Africa.
2. Construction of the wheat development model based on the GDD and BT. The temperature data used are Julian day 1-300 of the year 1995. The wheat development model is run for all grids of France and the surroundig countries. The model covers longitudes of -5 to 10 degrees and latitudes of 41-52 degrees.
3. Construction of a raster with harvest date values (Julian day in 1995) and mapping based on this raster.

### Project objectives

Assessment wheat harvest dates in France (although the model is applicable to any European country, and more crops can be added)

### Vector data

`Two files were created that indicate the x and y coordinates of the map of France, using the R package 'mapdata'. The coordinates of France were form mapdata were split into polygons. Two polygons were constructed: the French mainland and the island Corsica. Some very small islands in the Atlantic ocean were removed.`

### Raster data

```Daily average temperature data, 0.5 x 0.5 grids, 1995-2013, for Europe, West-Asia and North-Africa.

### Other data sources

```http://en.wikipedia.org/wiki/Growing_degree-day

http://www.fao.org/docrep/006/y4011e/y4011e06.htm```

## Computational implementation

The following code is used to make a 0.5 degrees raster of France:

```# How to make a 0.5 x 0.5 raster image of France?
library(mapdata)
# Selection of the region of interest:
map('worldHires', region = 'France', xlim = c(-5,10), ylim = c(41,52)) # Creates a map of France
Francemapdata <- map('worldHires', region = 'France', xlim = c(-5,10), ylim = c(41,52))
plot(Francemapdata) # Visualises the coordinates

# Francemapdata cannot be put in the function rasterize()
# Therefore, we plot the longitude and latitude of the coordinates
print(Francemapdata\$x)
print(Francemapdata\$y)
Matrixfrance <- cbind(Francemapdata\$x,Francemapdata\$y)
# The coordinates are a bit messy (NA's, coordinates not ranked in subsequent order).
# Data have to be corrected manually for France. I saved the data for the French mainland
# and Corsica as .csv files.```

The .csv files from the previous code are now objects called 'mainland' and 'Corsica' in the code below.

```#####################################################################################################
# Course Spatio temporal modelling                                                                  #
# January 2014                                                                                      #
# Author: Aart van der Linden                                                                       #
#                                                                                                   #
# Aim of the model: Simulate maturation and harvest of wheat based on gridded daily temperature     #
# data for France                                                                                   #
#####################################################################################################

# Packages required
library(ncdf)
library(raster)
library(rasterVis)
library(sp)
library(colorspace)
library(maptools)
library(maps)
library(mapdata)

# Basic geographical information

longbase <- -5          # left border longitude of map (degrees)
latbase  <- 41          # bottom border latitude of map (degrees)

nrlonggrids <- 30       # amount of 0.5 grids to the east
nrlatgrids  <- 24       # amount of 0.5 grids to the north
nrmaxtime   <- 300      # time series (days)

# Basic crop parameters
TSUM  <- 1600           # degree days required for maturation of wheat (http://en.wikipedia.org/wiki/Growing_degree-day)
BASET <- 4.5            # base temperature wheat (http://www.fao.org/docrep/006/y4011e/y4011e06.htm)

#Opening the temperature data file

X <- open.ncdf("C:/Weather_data/tg_0.50deg_reg_1995-2013_v9.0.nc")   # Open a weather dataset for France with daily mean temperatures
# Dataset starts at 40 W longitude and at 25 N latitude

# Create a raster with grids that belong to France

Mainlandpolygon <- Polygons(list(Polygon(Mainland)),1)                   # Polygon for the French mainland
Corsicapolygon  <- Polygons(list(Polygon(Corsica)),2)                    # Polygon for the island Corsica

MapFrance <- SpatialPolygons(list(Mainlandpolygon, Corsicapolygon)  )    # Join the two polygons

plot(MapFrance, xlim=c(longbase,(longbase+nrlonggrids/2)), ylim=c(latbase,(latbase+nrlatgrids/2)))  # Plot the map of France with the mainland and Corsica

# rasterizing polygons: from a detailed polygon to a coarse grid

r1 <- raster(ncol=nrlonggrids, nrow=nrlatgrids)                          # Construct a raster with 30 x 24 grids
r2 <- raster()
res(r2)=c(0.5,0.5)
e <- extent(longbase, (longbase+nrlonggrids/2), latbase, (latbase+nrlatgrids/2))      # Geographical coordinates of France
rc <- crop(r2, e)                                                                     # Creates a cropped image rc

r1 <- rasterize(MapFrance, rc, fun='min', background = 0.0)               # Convert the map of France into a raster
plot(r1,xlim=c(longbase,(longbase+nrlonggrids/2)), ylim=c(latbase,(latbase+nrlatgrids/2)))  # Plot the raster of France with the mainland and Corsica, Figure 3.

#####################################################################
#                      Model maturity wheat                         #
#####################################################################

GDD        = array(dim=c(nrlonggrids,nrlatgrids,nrmaxtime+1))            # array for growing degree days
MAT        = array(dim=c(nrlonggrids,nrlatgrids,nrmaxtime))              # array for heading (0=no, 1 =yes)
MATDATE    = matrix(nrow=nrlonggrids, ncol=nrlatgrids)                   # time after 1 jan for heading in days

Tmean      = NULL                                                        # statements that time, longitude and latitude start at 1
time       = NULL                                                        # relative to the left-bottom corner
time[1]    = 1
long       = NULL
long[1]    = 1
lat        = NULL
lat[1]     = 1

for (x in 1:nrlonggrids) {                                               # model run for longitude (x), latitude(y) and time(i)
for (y in 1:nrlatgrids) {
for (i in 1:nrmaxtime) {

long[x+1] = long[x]+1                                                # longitude, latitude and time each increase 1 after a step
lat[y+1]  = lat[y]+1
time[i+1] = time[i]+1

GDD[x,y,1] = 0                                                       # inital value for GDD is zero
MAT[x,y,1] = 0                                                       # inital value for FLO is zero

Tmean[i] <- get.var.ncdf(X, start= c(2*longbase+81+long[x],2*latbase-50+lat[y],time[i]), count=c(long[x],lat[y],time[i]))  # reads Tmean from file

if(Tmean[i]>BASET) GDD[x,y,i+1] = GDD[x,y,i] + (Tmean[i]-BASET) else GDD[x,y,i+1] = GDD[x,y,i]   # calculation GDD over time, longitude and latitude
if(GDD[x,y,i]>TSUM) MAT[x,y,i] = 1 else MAT[x,y,i] = 0                                       # calculation FLO over time, longitude and latitude

MATDATE[x,y] = time[i]-sum(MAT[x,y,1:i])+1                                                # calculation FLOWERDATE over longitude and latitude
if(Tmean[1]<(-90)) MATDATE[x,y] = NaN                                                     # exclude grids without weather data

}
}
}

MATDATE <- rotate(MATDATE)                           # turn the flowerdate matrix three times 90 degrees to get the correct map
MATDATE <- rotate(MATDATE)
MATDATE <- rotate(MATDATE)

rharvest = NULL                                             # creates a raster with values for harvest dates
rharvest <- raster(ncol=nrlonggrids, nrow=nrlatgrids,
xmn = longbase, xmx = (longbase+nrlonggrids/2), ymn = latbase, ymx = (latbase+nrlatgrids/2))

values(rharvest) <- MATDATE

Legend <- seq(152, nrmaxtime, 30.5)                            # legend for map

boundaries <- map('worldHires',  fill=TRUE,                # get country borders
xlim=c(longbase,(longbase+nrlonggrids/2)), ylim=c(latbase,(latbase+nrlatgrids/2)),
plot=FALSE)

## read the map2SpatialPolygons help page for details
IDs <- sapply(strsplit(boundaries\$names, ":"), function(x) x[1])
Borders <- map2SpatialPolygons(boundaries, IDs=IDs, proj4string=CRS(projection(rharvest)))  # country borders in spatial polygon

# Figure 4

levelplot(rharvest,
main  ="Wheat harvest date (day of the year)",
margin=FALSE,
xlab  =list('Longitude', fontface='bold'),
ylab  =list('Latitude', rot=45, fontface='bold'),
at    =Legend,
par.settings= rasterTheme(region = c('lawngreen','green3','darkgreen','brown')),

) + layer(sp.polygons(Borders))

r1[r1 < 0.5] <- NA

boundaries <- map('worldHires',region='France',  fill=TRUE,
xlim=c(longbase,(longbase+nrlonggrids/2)), ylim=c(latbase,(latbase+nrlatgrids/2)),
plot=FALSE)

## read the map2SpatialPolygons help page for details
IDs <- sapply(strsplit(boundaries\$names, ":"), function(x) x[1])
Borders <- map2SpatialPolygons(boundaries, IDs=IDs, proj4string=CRS(projection(rharvest)))

# Figure 5.
main  ="Wheat harvest date (day of the year)",
margin=FALSE,
xlab  =list('Longitude', fontface='bold'),
ylab  =list('Latitude', rot=45, fontface='bold'),
at    =Legend,
par.settings= rasterTheme(region = c('lawngreen','green3','darkgreen','brown')),

) + layer(sp.polygons(Borders))

Colourhist <- rep(c('lawngreen','green3','darkgreen','brown'),c(30.5,30.5,30.5,30.5))

histogram(mask, main= "Histogram wheat harvest dates", xlab = "day of the year",
par.settings= rasterTheme(region = c('lawngreen','green3','darkgreen','brown')),
col="goldenrod1")

#gc() # empty RAM

proc.time()

# Plot a time series (example of the data), Figure 1.
BBB <- get.var.ncdf(X, start= c(72,33,1), count=c(1,1,6940))         # 72 = longitude -5, 33 = latitude 41 (calibrated on the island Corsica)
XXX <- c(seq(from = 1, to = 6940, by =1))                            # Make a time series
plot(BBB~XXX, ylim= c(-10,35), xlab="time", ylab="temperature (degrees Celsius)", pch=16, cex=0.3, col = "red")                                       # Plot the time series

AAA <- get.var.ncdf(X, start= c(72,33,1), count=c(30,24,1))          # Data for the map of France (turned version)

AAA[AAA<(-99.0)]<- NaN

rotate <- function(AAA) t(apply(AAA, 2, rev))                        # Define the rotation
AAA <-rotate(AAA)                                                    # Turn the map three times 90 degrees
AAA <-rotate(AAA)
AAA <-rotate(AAA)

r <- raster(ncol=30, nrow=24)                                        # Create raster, equal to the size of France
values(r) <- AAA                                                     # Put the data for France in the raster

# Figure 2.
plot(r)                                                              # Plot the raster```

## RESULTS and DISCUSSION

An example of a time series for one coordinate is given in figure 1, a map of France with average temperatures for one day is given in figure 2. Those two figures visualise the big dataset a bit.

Figure 1. Average daily temperature of the coordinate -5 W, 41 N, from 1995 to 2013.

Figure 2. Average temperature for France (longitude (-5 W - 10 E), latitude (41 N - 52N) at the 1st January 1995.

A raster of France, including the island Corsica was made based on the data of the mapdata package in r, using the command rasterize()

Figure 3. Grid cells (0.5 x 0.5 degrees) of France and Corsica. The colours indicate the number of the polygons (1= main land, 2 = Corsica)

The research question was when French farmers can harvest the wheat (Fig. 4).

Figure 4. Modelled harvest dates of wheat in France and neighbouring countries in 1995. Light green = June; Lawn green = July; Dark green = August and Brown red is September. White spots indicate places that are sea, ocean or land with harvest after September.

France is plotted in Fig. 4, but also the neighbouring countries. Therefore, we apply the mask of Fig. 3 to Fig. 4.

Figure 5. Modelled harvest dates of wheat in France in 1995.

Figure 6. Histogram of the modelled wheat harvest dates for France.

We can also compare the variability in harvest dates between subsequent years.

Figure 7. Wheat harvest in 1995 (upper figure) and 1996.

Another option is to plot a new region. Here the initial coordinate is set at -15 W longitude and 37 N latitude. Adapting those two numbers results in Fig. 8.

Figure 8. Wheat harvest date 1995 with (-15,37) as an initial coordinate