Palm oil and upland rice

Field experiment

The experimental set-up consisted of rice plots grown under full sun or artificial shade. A net placed above the rice plots were used as artificial shade. In order to get a realistic level of shade, the nets were set-up to intercept as much light as intercepted by palm trees.

Modelling

The ARCHIMED model was used to simulate the light transmitted by the palm trees to a plane at 1 meter above the ground. This modelling approach help to compute the light transmitted according to any planting design.

ARCHIMED initializations

path_archimed= "D:/OneDrive/Travail_AMAP/ARCHIMED/Archimed_feb2019"

timestep= 30*60 # 30 minutes
clearness= 0.75
sim_date= data.frame(date= as.POSIXct("2019-03-20"), clearness= 0.5)

Hour= seq(from= as.POSIXct("2019-01-01 05:00:00"), by= timestep, 
          to = as.POSIXct("2019-01-01 18:30:00"))

# Dummy temperatures in the day:
Tair= archimedR::temperature_dynamic(Tmax = 30, Tmin = 18, per = (24*60*60)/timestep,
                                     shift = 3*pi/2)[10:37]

Meteo= 
  expand.grid(date= sim_date$date, hour_start= Hour)%>%
  arrange(date,hour_start)%>%
  mutate(hour_end= format(hour_start+timestep,"%H:%M:%S"),
         hour_start= format(hour_start,"%H:%M:%S"),
         temperature= Tair,
         relativeHumidity= 90,
         clearness= clearness)%>%
  mutate(date= format(date, "%Y/%m/%d"), wind= 3)

data.table::fwrite(
  Meteo,sep=";",
  file = file.path(path_archimed,'app_parameters/meteo/meteo_palm_rice.csv'))

archimedR::set_config(file = file.path(path_archimed,"app_parameters/ArchimedConfiguration.properties"),
                      parameter = "meteoFile", value = "app_parameters/meteo/meteo_palm_rice.csv")

site_location= list(Latitude= 0.9261918, Altitude = 45)
set_config(file = file.path(path_archimed,"app_parameters/ArchimedConfiguration.properties"),
           parameter = "latitude", value = site_location$Latitude)
set_config(file = file.path(path_archimed,"app_parameters/ArchimedConfiguration.properties"),
           parameter = "altitude", value = site_location$Altitude)

archimedR::set_config(file = file.path(path_archimed,"app_parameters/ArchimedConfiguration.properties"),
                      parameter = "skySectors", value = 46)
archimedR::set_config(file = file.path(path_archimed,"app_parameters/ArchimedConfiguration.properties"),
                      parameter = "numberOfPixels", value = 200000)

sim_folder= "output_palm_rice"
set_config(file = file.path(path_archimed,"app_parameters/ArchimedConfiguration.properties"),
           parameter = "outputDirectory", value = sim_folder)

Making the 3D scene using VPalmr (or VPALM-IDE)

To simulate the different designs, we need first the 3D mock-ups of the palm trees (OPF files) and their planting design (OPS files).

We only need an average palm tree for this study, so we use the parameters from an average palm tree. The scene is made using VPalmr. A less advanced user may prefer to make the scenes using VPALM-IDE instead.

Importing the data

Several files are imported :

• The planting designs data.frames (*2).

• The parameters input to VPALM, for each physiological age.

Creating and writing the OPF and OPS:

planting_designs= list.files("Designs/planting_design", full.names = TRUE)
palm_param_files= list.files("Designs/Palm_parameters", full.names = TRUE)
param_df= 
  expand.grid(param_file= palm_param_files, 
              design= planting_designs,stringsAsFactors = FALSE)%>%
  mutate(name= basename(.data$design)%>%gsub(".csv","",.))

# Cleaning: 
if(dir.exists("Designs/Vpalm_outputs")){
  unlink("Designs/Vpalm_outputs", recursive = TRUE)
}
# Making the OPF and first OPS: 
designs= 
  lapply(1:length(palm_param_files), function(x){
    params= readRDS(param_df[x,1])
    design= data.table::fread(file = param_df[x,2], data.table = FALSE)
    # Making the scene:
    scene= 
      Vpalmr::make_scene(data = params,
                         path = file.path("Designs","Vpalm_outputs"), 
                         AMAPStudio = "Designs/vpalm.jar",
                         plot_design= design)
    # Remove whole average:
    list.files(file.path("Designs","Vpalm_outputs","scenes"),
               full.names = TRUE)%>%
      .[grep("Average_MAP",.)]%>%
      unlink(x = .)
    
    list.files(file.path("Designs","Vpalm_outputs","scenes","opf"),
               full.names = TRUE)%>%
      .[grep("Average_Average_MAP",.)]%>%
      unlink(x = .)
  })
# Rename scene folder:
file.rename("Designs/Vpalm_outputs/scenes","Designs/Vpalm_outputs/scenes_palm_rice")
# Rename OPS:
list.files(file.path("Designs","Vpalm_outputs","scenes_palm_rice"),
           full.names = TRUE)%>%
  .[grep(".ops",.)]%>%
  file.rename(., {file.path(dirname(.),paste0(basename(.)%>%gsub("\\.ops",'',.),"_",param_df[1,3],".ops"))})

Adding the other designs OPS

param_df_others= param_df[-c(1:length(palm_param_files)),]

lapply(1:nrow(param_df_others), function(x){
  design= data.table::fread(file = param_df_others[x,2], data.table = FALSE)
  params= readRDS(param_df_others[x,1])
  prog= unique(params$input$Parameter$Progeny)
  map= params$input$MAP_requested
  
  format_ops(design = design, Progeny = prog,
             map = map, average = TRUE)%>% 
    write_ops(file.path("Designs","Vpalm_outputs","scenes_palm_rice", 
                        paste0(prog, "_MAP_", map,"_",param_df_others[x,3], ".ops")),
              overwrite = TRUE)  
})

ARCHIMED simulations

First we have to copy the 3D scene onto the ARCHIMED installation folder

scene_path= file.path(path_archimed,"app_parameters/scene/scenes_palm_rice")
if(!dir.exists(scene_path)){dir.create(scene_path)}
file.copy("Designs/Vpalm_outputs/scenes_palm_rice",dirname(scene_path),recursive = TRUE)

Then, running the simulation for each scene

parameters= list("xMin","yMin","xMax","yMax")
design_files= list.files(scene_path)%>%.[grep('.ops',.)]
# Cleaning the Output folder:
unlink(file.path(path_archimed,sim_folder), recursive = T)

#  Reading the model manager: 
model_manager= fread(file.path(path_archimed,"app_parameters/models_manager/models_palmtree.csv"), data.table = F)

# NB: if ARCHIMED did not complete the simulation and returned an error, open the OPS file and 
# change the OPF files for ones that worked previously, run a simulation, change back the OPFs
# to the right ones and it may work.
tmp=
  lapply(design_files, function(x){
    # Set the scene boundaries:
    density= x%>%strsplit(.,"_")%>%unlist()%>%.[grep(".ops",.)]%>%gsub(".ops","",.)
    planting_design= data.table::fread(planting_designs[grep(density,basename(planting_designs))],
                                       data.table = FALSE)
    tmp=
      lapply(parameters, function(z){
        value= unique(planting_design[,grep(z,colnames(planting_design), ignore.case = TRUE)])
        set_config(file = file.path(path_archimed,"app_parameters/ArchimedConfiguration.properties"),
                   parameter = z, value = value)
      })
    # Set the ops name: 
    path= file.path('app_parameters', "scene","scenes_palm_rice", x)
    set_config(file = file.path(path_archimed,"app_parameters/ArchimedConfiguration.properties"),
               parameter = "file", value = path)
    
    # Updating the model manager:
    model_manager$Node_id[model_manager$Group=="pavement"]= 2 + nrow(planting_design)
    fwrite(model_manager, 
           file.path(path_archimed,"app_parameters/models_manager/models_palmtree.csv"),
           sep = ";", quote = FALSE)
    
    # Run ARCHIMED:
    run_archimed(exe = file.path(path_archimed,"archimed-lib_florian-1.0.0.jar"), memory= 16384, config= "Archimed")
  })%>%unlist()

# Test if the simulations ran successfully:
if(any(!tmp)){
  warning(paste("Error during ARCHIMED simulation for:\n *",paste(names(tmp)[!tmp], collapse = "\n * ")))
}

Light transmitted to the soil according to the distance to the tree

plane_df_step%>%
  ggplot(aes(x=dist_tree, y=irradiance))+
  facet_wrap(~paste0("Density: ",dens, ", Age: ",MAP)) +
  geom_point(aes(color= irradiance*10^6), show.legend = FALSE)+
  labs(title = 'Hour: {frame_time}', x = 'Distance to the closest palm tree (m)', 
       y = 'Irradiance (W m-2)')+
  gganimate::transition_time(Date)+
  gganimate::ease_aes('linear')

gganimate::anim_save(filename = 'points_light_rice_distance_to_tree.gif',path = "outputs")

Integrated over a day:

• Low density:

grid_df_day%>%
  filter(Treatment== "DA1_MAP_72_inter11xintra9_d101" |
           Treatment== "DA1_MAP_90_inter11xintra9_d101")%>%
  ggplot(aes(x=x, y=y))+
  facet_wrap(~paste0("Density: ",dens, ", Age: ",MAP)) +
  geom_tile(aes(fill=Intercepted))+
  coord_fixed()+
  labs(x = 'x (m)', y = 'y (m)',fill= "Intercepted (MJ m-2 d-1)") +
  theme(legend.position="bottom")+
  scale_fill_gradientn(
    colours = rev(RColorBrewer::brewer.pal(5,name = "OrRd")),
    breaks=seq(0,13,3), 
    limits = c(0,13))

ggplot2::ggsave(filename = 'Interception_daily_low_density.png',device = "png",
                path = "outputs",
                width = 17,height=12,units="cm",dpi=200,scale = 1.2)

• High density:

grid_df_day%>%
  filter(Treatment== "DA1_MAP_72_inter10xintra8_d125" |
           Treatment== "DA1_MAP_90_inter10xintra8_d125")%>%
  ggplot(aes(x=x, y=y))+
  facet_wrap(~paste0("Density: ",dens, ", Age: ",MAP)) +
  geom_tile(aes(fill=Intercepted))+
  coord_fixed()+
  labs(x = 'x (m)', y = 'y (m)',fill= "Intercepted (MJ m-2 d-1)") +
  theme(legend.position="bottom")+
  scale_fill_gradientn(
    colours = rev(RColorBrewer::brewer.pal(5,name = "OrRd")),
    breaks=seq(0,13,3), 
    limits = c(0,13))

ggplot2::ggsave(filename = 'Interception_daily_high_density.png',device = "png",
                path = "outputs",
                width = 17,height=12,units="cm",dpi=200,scale = 1.2)

• Distance to the tree

grid_df_day%>%
  ggplot(aes(x=dist_tree, y=Intercepted))+
  facet_wrap(~paste0("Density: ",dens, ", Age: ",MAP)) +
  geom_point(aes(color= Intercepted), show.legend = FALSE)+
  labs(title = 'Daily transmitted light', x = 'Distance to the closest palm tree (m)', 
       y = 'Irradiance (MJ m-2 d-1)')+
  scale_colour_gradientn(
    colours = rev(RColorBrewer::brewer.pal(5,name = "OrRd")),
    breaks=seq(0,13,3), 
    limits = c(0,13))

ggplot2::ggsave(filename = 'Interception_distance_to_tree_daily.png',device = "png",
                path = "outputs",
                width = 17,height=12,units="cm",dpi=200,scale = 1.2)

Light available for the rice:

We considered that the rice crop would be planted at a distance of 3m minimum to enable field workers to manage the palm, and computed hereafter the daily average light available on the strip:

grid_df_day_row= 
  grid_df_day%>%
  group_by(Treatment,x)%>%
  summarise(Intercepted= mean(Intercepted), 
            dens= unique(dens), MAP= unique(MAP),
            dist_tree_x= unique(dist_tree_x),
            x_tree_1= unique(x_tree_1),
            x_tree_2= unique(x_tree_2))

grid_df_day_row_poly=
  grid_df_day_row%>%
  filter(dist_tree_x<3)%>%
  mutate(tree= ifelse(x<(x_tree_2-3),"Palm 1","Palm 2"))

grid_df_day_row%>%
  ggplot()+
  facet_wrap(~paste0("Density: ",dens)) +
  geom_ribbon(data= grid_df_day_row_poly,
              aes(x= x, ymin=  -Inf, ymax= Inf,
                  y= Intercepted, group= tree), alpha= 0.3)+
  geom_line(aes(x=x, y= Intercepted, color= MAP), size= 1)+
  labs(x = 'x (m)', y= "Intercepted (MJ m-2 d-1)")

ggplot2::ggsave(filename = 'Interception_daily_transect.png',device = "png",
                path = "outputs",
                width = 17,height=12,units="cm",dpi=200,scale = 1.2)

The light transmitted to the interrow is in average :

interc= 
  grid_df_day_row%>%
  filter(dist_tree_x>3 & (x>x_tree_1&x<x_tree_2))%>%
  rename(Density= dens)%>%
  group_by(Density,MAP)%>%
  summarise(Transmitted= mean(Intercepted), sd= sd(Intercepted))
knitr::kable(interc, caption = "Radiation transmitted by the palm trees to the interrow (>3m from the plam)")

3D visualizations

Here is a visualization of the simulated light interception in the 3D scene from above:

Light interception in the 3D scene