User Tools

Site Tools


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. 
http://www.ecad.eu/download/ensembles/download.php

Other data sources

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

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

METHOD

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
Mainland <- read.csv(file="M:/R/Francemainland.csv",head=FALSE,sep=",")  # Data file mainland from package "mapdata"
Corsica  <- read.csv(file="M:/R/Francecorsica.csv",head=FALSE,sep=",")   # Data file mainland from package "mapdata"
 
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
mask <- mask(rharvest, r1)
 
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.
levelplot(mask, 
          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

wikistudholland/utwente_aart.txt · Last modified: 2021/01/20 20:36 (external edit)