risk_yield_loss.R 27.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
################################################################################
## This script calculates the risk of forecast yields falling below           ##
## historical reference yields by a user-specified threshold.                 ##
##                                                                            ##
## It returns:                                                                ##
## - Gridded data per crop and irrigation status of                           ##
##   * Risk of yield loss below threshold (%)                                 ##
##   * Harvested area affected of such loss (ha) grouped by risk category     ##
##     (0-20%, 20-40%, 40-60%, 60-80%, 80-100%)                               ##
##                                                                            ##
## Input variables (tunable knobs):                                           ##
## - start year of historical reference period                                ##
## - end year of historical reference period                                  ##
## - yield loss threshold (%)                                                 ##
## - forecast year to compare to historical reference yield                   ##
## - irrigation intervention scenario                                         ##
## - fertilizer intervention scenario                                         ##
## - sowing date intervention scenario                                        ##
##                                                                            ##
## Further settings:                                                          ##
## - experiment name                                                          ##
## - loss probability threshold (%), currently not used, for aggregation to   ##
##   spatial units                                                            ##
################################################################################
25

26
27
28
29
30
rm(list = ls())

################################################################################
## Default values for command line parameters                                 ##
# Experiment name (must be defined in rsync_setup.R)
31
experiment <- "August2021_experiment"
32
# Reference period
33
34
hist_ref_startyear <- 2016
hist_ref_endyear <- 2020
35
# Threshold for yield loss compared to average yield during reference period
36
yield_loss <- 10 # in percent
37
38
39
40
# Fraction (in %) of weather options that must show yield loss of at least
# yield_loss %, not used at the moment
loss_probability <- (2/3) * 100
# Forecast year (2021 or 2022) which is compared to historical yield
41
forecast_year <- 2021
42
# Intervention scenario options
43
44
45
irrigation_scenario <- "reference"
fertilizer_scenario <- "reference"
sowing_scenario <- "reference"
46
################################################################################
47

48
49
################################################################################
## Load R packages                                                            ##
50
library(ncdf4)
51
52
library(raster)
################################################################################
53

54
55
56
57
58
59
60
61
62
63
64
65
66
67
################################################################################
## Setup of risk categories                                                   ##
risk_category_bnds <- seq(0, 100, by = 20)
risk_category_names <- character(length(risk_category_bnds) - 1)
for (r in seq_along(risk_category_bnds)[- 1]) {
  risk_category_names[r - 1] <- paste0(
    "Risk: ",
    risk_category_bnds[r - 1],
    "-",
    risk_category_bnds[r],
    "%"
  )
}
################################################################################
68

69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
################################################################################
## Read command line parameters (if any were provided)                        ##
# Will always take the first value if parameters are provided more than once
if (length(commandArgs(trailingOnly = TRUE)) > 0) {
  if (any(grepl("hist_ref_startyear", commandArgs(trailingOnly = TRUE)))) {
    param <- strsplit(
      grep(
        "hist_ref_startyear",
        commandArgs(trailingOnly = TRUE),
        value = TRUE
      )[1],
      split = "="
    )
    startyear <- as.integer(
      grep("hist_ref_startyear", unlist(param), invert = TRUE, value = TRUE)
    )
    if (is.finite(startyear)) {
      print(
        paste(
          "Parameter hist_ref_startyear provided as command line argument:",
          startyear
        )
      )
92
93
      hist_ref_startyear <- startyear
    } else {
94
      stop(paste("Invalid parameter hist_ref_startyear:", sQuote(startyear)))
95
96
    }
  }
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
  if (any(grepl("hist_ref_endyear", commandArgs(trailingOnly = TRUE)))) {
    param <- strsplit(
      grep(
        "hist_ref_endyear",
        commandArgs(trailingOnly = TRUE),
        value = TRUE
      )[1],
      split = "="
    )
    endyear <- as.integer(
      grep("hist_ref_endyear", unlist(param), invert = TRUE, value = TRUE)
    )
    if (is.finite(startyear)) {
      print(
        paste(
          "Parameter hist_ref_endyear provided as command line argument:",
          endyear
        )
      )
116
117
      hist_ref_endyear <- endyear
    } else {
118
      stop(paste("Invalid parameter hist_ref_endyear:", sQuote(endyear)))
119
120
    }
  }
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
  if (any(grepl("yield_loss", commandArgs(trailingOnly = TRUE)))) {
    param <- strsplit(
      grep("yield_loss", commandArgs(trailingOnly = TRUE), value = TRUE)[1],
      split = "="
    )
    loss <- as.double(
      grep("yield_loss", unlist(param), invert = T, value = TRUE)
    )
    if (is.finite(loss)) {
      print(
        paste(
          "Parameter yield_loss provided as command line argument:",
          loss
        )
      )
136
137
      yield_loss <- loss
    } else {
138
      stop(paste("Invalid parameter yield_loss:", sQuote(loss)))
139
140
    }
  }
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
  if (any(grepl("loss_probability", commandArgs(trailingOnly = TRUE)))) {
    param <- strsplit(
      grep(
        "loss_probability",
        commandArgs(trailingOnly = TRUE),
        value = TRUE
      )[1],
      split = "="
    )
    loss <- as.double(
      grep("loss_probability", unlist(param), invert = TRUE, value = TRUE)
    )
    print(
      paste(
        "Parameter loss_probability provided as command line argument:",
        loss
      )
    )
159
160
    loss_probability <- loss
  }
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
  if (any(grepl("forecast_year", commandArgs(trailingOnly = TRUE)))) {
    param <- strsplit(
      grep("forecast_year", commandArgs(trailingOnly = TRUE), value = TRUE)[1],
      split = "="
    )
    year <- as.integer(
      grep("forecast_year", unlist(param), invert = TRUE, value = TRUE)
    )
    if (is.finite(year)) {
      print(
        paste(
          "Parameter forecast_year provided as command line argument:",
          year
        )
      )
176
177
      forecast_year <- year
    } else {
178
      stop(paste("Invalid parameter forecast_year:", sQuote(year)))
179
180
    }
  }
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
  if (any(grepl("irrigation_scenario", commandArgs(trailingOnly = TRUE)))) {
    param <- strsplit(
      grep(
        "irrigation_scenario",
        commandArgs(trailingOnly = TRUE),
        value = TRUE
      )[1],
      split = "="
    )
    scen <- grep(
      "irrigation_scenario",
      unlist(param),
      invert = TRUE, value = TRUE
    )
    if (nchar(scen) > 0) {
      print(
        paste(
          "Parameter irrigation_scenario provided as command line argument:",
          scen
        )
      )
202
203
      irrigation_scenario <- scen
    } else {
204
      stop(paste("Invalid parameter irrigation_scenario:", sQuote(scen)))
205
206
    }
  }
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
  if (any(grepl("fertilizer_scenario", commandArgs(trailingOnly = TRUE)))) {
    param <- strsplit(
      grep(
        "fertilizer_scenario",
        commandArgs(trailingOnly = TRUE),
        value = TRUE
      )[1],
      split = "="
    )
    scen <- grep(
      "fertilizer_scenario",
      unlist(param),
      invert = TRUE, value = TRUE
    )
    if (nchar(scen) > 0) {
      print(
        paste(
          "Parameter fertilizer_scenario provided as command line argument:",
          scen
        )
      )
228
229
      fertilizer_scenario <- scen
    } else {
230
      stop(paste("Invalid parameter fertilizer_scenario:", sQuote(scen)))
231
232
    }
  }
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
  if (any(grepl("sowing_scenario", commandArgs(trailingOnly = TRUE)))) {
    param <- strsplit(
      grep(
        "sowing_scenario",
        commandArgs(trailingOnly = TRUE),
        value = TRUE
      )[1],
      split = "="
    )
    scen <- grep("sowing_scenario", unlist(param), invert = TRUE, value = TRUE)
    if (nchar(scen) > 0) {
      print(
        paste(
          "Parameter sowing_scenario provided as command line argument:",
          scen
        )
      )
250
251
      sowing_scenario <- scen
    } else {
252
      stop(paste("Invalid parameter sowing_scenario:", sQuote(scen)))
253
254
    }
  }
255
  if (any(grepl("^experiment=", commandArgs(trailingOnly = TRUE)))) {
256
    param <- strsplit(
257
      grep("^experiment=", commandArgs(trailingOnly = TRUE), value = TRUE)[1],
258
259
      split = "="
    )
260
    exp <- grep("^experiment$", unlist(param), invert = TRUE, value = TRUE)
261
262
263
264
    if (nchar(exp) > 0) {
      print(
        paste("Parameter experiment provided as command line argument:", exp)
      )
265
266
      experiment <- exp
    } else {
267
      stop(paste("Invalid parameter experiment:", sQuote(exp)))
268
269
270
    }
  }
}
271
################################################################################
272

273
274
################################################################################
## Rsync functionality                                                        ##
275
source("rsync_setup.R")
276
## Utility functions to read/process data                                     ##
277
278
279
280
source("get_forecast_options.R")
source("get_weather_options.R")
source("download_yield_data.R")
source("average_yield.R")
281
################################################################################
282

283
284
285
286
287
################################################################################
## Valid experiment?                                                          ##
if (!experiment %in% names(historical_remote_dir) ||
  !experiment %in% names(forecast_remote_dir)
) {
288
289
  stop(paste("Invalid experiment", sQuote(experiment)))
}
290
################################################################################
291

292
293
294
295
################################################################################
## Calculations of historical reference                                       ##
# Download historical time series data
if (!file.exists("output")) {
296
297
  dir.create("output")
}
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
hist_timeseries_filename <- download_yield_data(
  runtype = "historical",
  localname = "output/hist_timeseries.nc",
  experiment = experiment,
  irrigation_scenario = NULL,
  fertilizer_scenario = NULL,
  sowing_scenario = NULL,
  weather_year = NULL
)
if (!is.character(hist_timeseries_filename) ||
  !file.exists(hist_timeseries_filename)
) {
  stop(
    paste(
      "Error downloading historical time series data to calculate",
      "reference yield."
    )
  )
316
} else {
317
318
319
320
321
322
  print(
    paste(
      "Historical yield time series downloaded to temporary file",
      hist_timeseries_filename
    )
  )
323
}
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
# Calculate average yield for reference period
print(
  paste0(
    "Calculate average yield for reference period ",
    hist_ref_startyear, "-", hist_ref_endyear
  )
)
hist_ref_yield <- average_yield(
  filename = hist_timeseries_filename,
  avg_variable = "crop_yield",
  avg_startyear = hist_ref_startyear,
  avg_endyear = hist_ref_endyear,
  multiple_per_year = "mean",
  multiple_weight_variable = "harvested_area"
)
## Calculate threshold for yield loss
if (yield_loss < 0 || yield_loss > 100) {
  stop(
    paste(
      "yield_loss threshold must be between 0 and 100%.\n",
      "Provided value:", yield_loss
    )
  )
347
}
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
yield_loss_threshold <- hist_ref_yield * (100 - yield_loss) / 100
# Determine spatial extent and resolution of historical yield data
# Note: this will give a warning due to the structure of the variable but still
# works as expected
print("Expect a NetCDF warning about missing dimensions here...")
tmpraster <- raster(hist_timeseries_filename, band = 1, var = "crop_yield")
yield_raster <- raster(
  extent(tmpraster),
  resolution = res(tmpraster),
  crs = "+proj=longlat +datum=WGS84"
)
rm(tmpraster)
# Compare latitude axis direction of NetCDF
nc <- nc_open(hist_timeseries_filename)
hist_timeseries_flip <- (
  (nc$dim$lat$vals[1] < nc$dim$lat$vals[2]) !=
  (yFromRow(yield_raster, 1) < yFromRow(yield_raster, 2))
)
# Flip NetCDF data if axis has opposite direction
if (hist_timeseries_flip) {
  bakdim <- dim(yield_loss_threshold)
  bakdimnames <- dimnames(yield_loss_threshold)
  yield_loss_threshold <- yield_loss_threshold[, nc$dim$lat$len:1, , ]
  dim(yield_loss_threshold) <- bakdim
  dimnames(yield_loss_threshold) <- bakdimnames
  if (!is.null(bakdimnames)[[2]]) {
    dimnames(yield_loss_threshold)[[2]] <- rev(bakdimnames[[2]])
  }
}
nc_close(nc)
# Delete dowloaded file to conserve space
379
file.remove(hist_timeseries_filename)
380
################################################################################
381

382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
################################################################################
## Calculations of forecast yields and yield losses                           ##
# Get all available weather options for selected forecast scenario options
forecast_weather_options <- get_weather_options(
  experiment = experiment,
  remote_dir = forecast_remote_dir[experiment],
  irrigation_scenario = irrigation_scenario,
  fertilizer_scenario = fertilizer_scenario,
  sowing_scenario = sowing_scenario,
  include_filenames = FALSE
)
# Variables to collect information across forecast ensemble
forecast_below_threshold_count <- array(
  dim = dim(yield_loss_threshold),
  dimnames = dimnames(yield_loss_threshold)
)
ha_forecast_below_threshold <- array(
  dim = dim(yield_loss_threshold),
  dimnames = dimnames(yield_loss_threshold)
)
harvested_area_all <- array(
  dim = dim(yield_loss_threshold),
  dimnames = dimnames(yield_loss_threshold)
)
harvested_area_all_count <- array(0,
  dim = dim(yield_loss_threshold),
  dimnames = dimnames(yield_loss_threshold)
)
# Check how many weather options fall below yield_loss_threshold
for (wy in forecast_weather_options) {
  forecast_filename <- download_yield_data(
    runtype = "forecast",
    localname = "output/forecast.nc",
    experiment = experiment,
    irrigation_scenario = irrigation_scenario,
    fertilizer_scenario = fertilizer_scenario,
    sowing_scenario = sowing_scenario,
    weather_year = wy
  )
  if (!is.character(forecast_filename) || !file.exists(forecast_filename)) {
    stop(
      paste(
        "Error downloading forecast data for weather option:",
        sQuote(wy)
      )
    )
428
  } else {
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
    print(
      paste(
        "Forecast data for weather option", sQuote(wy),
        "downloaded to temporary file", forecast_filename
      )
    )
  }
  print(
    paste(
      "Calculate yield change in",
      forecast_year,
      "compared to historical reference yield"
    )
  )
  forecast_yield <- average_yield(
    filename = forecast_filename,
    avg_variable = "crop_yield",
    avg_startyear = forecast_year,
    avg_endyear = forecast_year,
    multiple_per_year = "sum",
    multiple_weight_variable = NULL
  )
  harvested_area <- average_yield(
    filename = forecast_filename,
    avg_variable = "harvested_area",
    avg_startyear = forecast_year,
    avg_endyear = forecast_year,
    multiple_per_year = "sum",
    multiple_weight_variable = NULL
  )
  # Compare latitude axis direction of NetCDF
  nc <- nc_open(forecast_filename)
  forecast_flip <- (
    (nc$dim$lat$vals[1] < nc$dim$lat$vals[2]) !=
    (yFromRow(yield_raster, 1) < yFromRow(yield_raster, 2))
  )
  # Flip NetCDF data if axis has opposite direction
  if (forecast_flip) {
    bakdim <- dim(forecast_yield)
    bakdimnames <- dimnames(forecast_yield)
    forecast_yield <- forecast_yield[, nc$dim$lat$len:1, , ]
    dim(forecast_yield) <- bakdim
    dimnames(forecast_yield) <- bakdimnames
    if (!is.null(bakdimnames)[[2]]) {
      dimnames(forecast_yield)[[2]] <- rev(bakdimnames[[2]])
    }
    bakdim <- dim(harvested_area)
    bakdimnames <- dimnames(harvested_area)
    harvested_area <- harvested_area[, nc$dim$lat$len:1, , ]
    dim(harvested_area) <- bakdim
    dimnames(harvested_area) <- bakdimnames
    if (!is.null(bakdimnames)[[2]]) {
      dimnames(harvested_area)[[2]] <- rev(bakdimnames[[2]])
    }
483
  }
484
485
486
  nc_close(nc)

  # Delete downloaded file to ensure next weather option is downloaded correctly
487
  file.remove(forecast_filename)
488
489

  if (any(dim(forecast_yield) != dim(yield_loss_threshold))) {
490
491
    stop("Forecast and reference yield data have different dimensions")
  }
492
493
494
495
496
497
498
  if (any(dim(harvested_area) != dim(yield_loss_threshold))) {
    stop(
      paste(
        "Forecast harvested area and reference yield data have different",
        "dimensions"
      )
    )
499
  }
500
501
502
503
504
505
506
507
  # Harmonize sequence of crops and cropping systems (rainfed/irrigated)
  if (!is.null(dimnames(forecast_yield)[[3]]) &&
    !is.null(dimnames(yield_loss_threshold)[[3]])
  ) {
    if (length(intersect(
      dimnames(forecast_yield)[[3]],
      dimnames(yield_loss_threshold)[[3]]
    )) != dim(yield_loss_threshold)[3]) {
508
509
510
511
512
      stop("Crop mismatch between forecast and reference yield data")
    }
    bakdim <- dim(forecast_yield)
    bakdimnames <- dimnames(forecast_yield)
    bakdimnames[[3]] <- dimnames(yield_loss_threshold)[[3]]
513
514
515
516
517
    crop_match <- match(
      dimnames(yield_loss_threshold)[[3]],
      dimnames(forecast_yield)[[3]]
    )
    forecast_yield <- forecast_yield[, , crop_match, ]
518
519
520
    dim(forecast_yield) <- bakdim
    dimnames(forecast_yield) <- bakdimnames
  }
521
522
523
524
525
526
527
  if (!is.null(dimnames(forecast_yield)[[4]]) &&
    !is.null(dimnames(yield_loss_threshold)[[4]])
  ) {
    if (length(intersect(
      dimnames(forecast_yield)[[4]],
      dimnames(yield_loss_threshold)[[4]]
    )) != dim(yield_loss_threshold)[4]) {
528
529
530
531
532
      stop("Cropping system mismatch between forecast and reference yield data")
    }
    bakdim <- dim(forecast_yield)
    bakdimnames <- dimnames(forecast_yield)
    bakdimnames[[4]] <- dimnames(yield_loss_threshold)[[4]]
533
534
535
536
537
    irr_match <- match(
      dimnames(yield_loss_threshold)[[4]],
      dimnames(forecast_yield)[[4]]
    )
    forecast_yield <- forecast_yield[, , , irr_match]
538
539
540
    dim(forecast_yield) <- bakdim
    dimnames(forecast_yield) <- bakdimnames
  }
541
542
543
544
545
546
547
548
549
550
551
552
553
  if (!is.null(dimnames(harvested_area)[[3]]) &&
    !is.null(dimnames(yield_loss_threshold)[[3]])
  ) {
    if (length(intersect(
      dimnames(harvested_area)[[3]],
      dimnames(yield_loss_threshold)[[3]]
    )) != dim(yield_loss_threshold)[3]) {
      stop(
        paste(
          "Crop mismatch between forecast harvested area and reference",
          "yield data"
        )
      )
554
555
556
557
    }
    bakdim <- dim(harvested_area)
    bakdimnames <- dimnames(harvested_area)
    bakdimnames[[3]] <- dimnames(yield_loss_threshold)[[3]]
558
559
560
561
562
    crop_match <- match(
      dimnames(yield_loss_threshold)[[3]],
      dimnames(harvested_area)[[3]]
    )
    harvested_area <- harvested_area[, , crop_match, ]
563
564
565
    dim(harvested_area) <- bakdim
    dimnames(harvested_area) <- bakdimnames
  }
566
567
568
569
570
571
572
573
574
575
576
577
578
  if (!is.null(dimnames(harvested_area)[[4]]) &&
    !is.null(dimnames(yield_loss_threshold)[[4]])
  ) {
    if (length(intersect(
      dimnames(harvested_area)[[4]],
      dimnames(yield_loss_threshold)[[4]]
    )) != dim(yield_loss_threshold)[4]) {
      stop(
        paste(
          "Cropping system mismatch between forecast harvested area and",
          "reference yield data"
        )
      )
579
580
581
582
    }
    bakdim <- dim(harvested_area)
    bakdimnames <- dimnames(harvested_area)
    bakdimnames[[4]] <- dimnames(yield_loss_threshold)[[4]]
583
584
585
586
587
    irr_match <- match(
      dimnames(yield_loss_threshold)[[4]],
      dimnames(harvested_area)[[4]]
    )
    harvested_area <- harvested_area[, , , irr_match]
588
589
590
    dim(harvested_area) <- bakdim
    dimnames(harvested_area) <- bakdimnames
  }
591
  if (!identical(is.na(forecast_yield), is.na(harvested_area))) {
592
593
    stop("Mismatch between forecast crop yield and harvested area")
  }
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
  # Crops/cells with significant yield loss
  forecast_below_threshold <- which(forecast_yield <= yield_loss_threshold)
  # Determine whether to set aggregation variables (because they are still NA)
  # or update them if set already
  set_value <- intersect(
    which(is.na(forecast_below_threshold_count)),
    forecast_below_threshold
  )
  update_value <- intersect(
    which(!is.na(forecast_below_threshold_count)),
    forecast_below_threshold
  )
  set_value_ha <- intersect(
    which(is.na(ha_forecast_below_threshold)),
    forecast_below_threshold
  )
  update_value_ha <- intersect(
    which(!is.na(ha_forecast_below_threshold)),
    forecast_below_threshold
  )
614
  crop_present <- which(!is.na(harvested_area))
615
616
617
618
619
620
621
622
623
  set_value_ha_sum <- intersect(which(is.na(harvested_area_all)), crop_present)
  update_value_ha_sum <- intersect(
    which(!is.na(harvested_area_all)),
    crop_present
  )
  # Update/set aggregation variables
  # Number of forecasts with yield_loss above threshold
  forecast_below_threshold_count[update_value] <-
    forecast_below_threshold_count[update_value] + 1
624
  forecast_below_threshold_count[set_value] <- 1
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
  # Harvested area of crops with yield loss above threshold
  ha_forecast_below_threshold[set_value_ha] <-
    harvested_area[set_value_ha]
  ha_forecast_below_threshold[update_value_ha] <-
    ha_forecast_below_threshold[update_value_ha] +
    harvested_area[update_value_ha]
  # Harvested area for all crops/cells regardless of yield loss
  harvested_area_all[set_value_ha_sum] <- harvested_area[set_value_ha_sum]
  harvested_area_all[update_value_ha_sum] <-
    harvested_area_all[update_value_ha_sum] +
    harvested_area[update_value_ha_sum]
  # Counter for harvested area for all crops
  harvested_area_all_count[crop_present] <-
    harvested_area_all_count[crop_present] + 1
  # Clean up
  rm(
    forecast_below_threshold,
    set_value,
    update_value,
    set_value_ha,
    update_value_ha,
    crop_present,
    set_value_ha_sum,
    update_value_ha_sum,
    forecast_yield,
    harvested_area
  )
}
# Calculate risk in each cell
yield_loss_risk <- forecast_below_threshold_count /
  length(forecast_weather_options) * 100
# Calculate harvested area at risk (average across weather options where yield
# loss is above threshold)
ha_forecast_below_threshold <-
  ha_forecast_below_threshold / forecast_below_threshold_count
zero <- which(forecast_below_threshold_count == 0)
ha_forecast_below_threshold[zero] <- NA
# Calculate harvested area regardless of yield loss risk (average across all
# weather options)
harvested_area_all_count[which(is.na(harvested_area_all))] <- NA
harvested_area_all <- harvested_area_all / harvested_area_all_count
harvested_area_all[which(harvested_area_all_count == 0)] <- NA
# Calculate risk category
yield_loss_risk_category <- array(
  dim = dim(yield_loss_risk),
  dimnames = dimnames(yield_loss_risk)
)
for (r in seq_along(risk_category_names)) {
  set_value <- which(yield_loss_risk >= risk_category_bnds[r])
  yield_loss_risk_category[set_value] <- r
}
if (min(risk_category_bnds) == 0) {
  # Add areas for which no risk was found
  set_value <- which(is.na(yield_loss_risk) & !is.na(harvested_area_all))
  ha_forecast_below_threshold[set_value] <- harvested_area_all[set_value]
  yield_loss_risk_category[set_value] <- 1
681
682
}

683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
# Split ha_forecast_below_threshold into bands by risk category
# This way, harvested areas can be related to risk category easily in CauseMos
ha_by_risk_category <- array(
  dim = c(
    dim(ha_forecast_below_threshold),
    length(risk_category_names)
  ),
  dimnames = c(
    dimnames(ha_forecast_below_threshold),
    list(risk_category_names)
  )
)
for (r in seq_along(risk_category_names)) {
  set_value <- which(yield_loss_risk_category == r)
  # Change index from single-band array to corresponding index in multi-band
  # array
  set_value_tgt <- set_value +
    (r - 1) * length(ha_forecast_below_threshold)
  # Assign only harvested areas with the corresponding risk category to the
  # respective band in ha_by_risk_category
  ha_by_risk_category[set_value_tgt] <-
    ha_forecast_below_threshold[set_value]
}
706

707
################################################################################
708

709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
################################################################################
## Prepare gridded output                                                     ##
# Gridded risk of yield loss of at least yield_loss percent
lon_dim <- ncdim_def(
  name = "longitude",
  units = "degrees_east",
  vals = xFromCol(yield_raster),
  longname = "Longitude"
)
lat_dim <- ncdim_def(
  name = "latitude",
  units = "degrees_north",
  vals = yFromRow(yield_raster),
  longname = "Latitude"
)
crop_dim <- ncdim_def(
  name = "crop",
  units = "",
  vals = 1:dim(yield_loss_risk)[3],
  longname = "Crop index",
  create_dimvar = FALSE
)
irrigation_dim <- ncdim_def(
  name = "irrigation",
  units = "",
  vals = 1:dim(yield_loss_risk)[4],
  longname = "Irrigation status index",
  create_dimvar = FALSE
)
# Time dimension only one year set to forecast_year
time_dim <- ncdim_def(
  name = "time",
  units = "year",
  vals = forecast_year,
  longname = "Forecast year"
)
yield_loss_risk_var <- ncvar_def(
  name = "yield_loss_risk",
  units = "%",
  dim = list(lon_dim, lat_dim, crop_dim, irrigation_dim, time_dim),
  missval = 1e20,
  longname = paste0(
    "Risk of yield loss of at least ",
    yield_loss, "%",
    " below historical reference yield"
  ),
  prec = "float",
)
# Gridded harvested area at risk of yield loss below threshold
# This has harvested areas split by risk category to be able to use risk
# category as a qualifier in CauseMos
risk_category_dim <- ncdim_def(
  name = "risk_category",
  units = "",
  vals = seq_along(risk_category_names),
  longname = "Risk category index",
  create_dimvar = FALSE
)
harvested_area_var <- ncvar_def(
  name = "harvested_area_at_risk",
  units = "ha",
  dim = list(
    lon_dim,
    lat_dim,
    crop_dim,
    irrigation_dim,
    risk_category_dim,
    time_dim
  ),
  missval = 1e20,
  longname = paste0(
    "Harvested area at risk of yield loss of at least ",
    yield_loss, "%",
    " below historical reference yield"
  ),
  prec = "float"
)
# Name variables
nchar_dim <- ncdim_def(
  name = "nchar",
  units = "",
  vals = seq_len(
    max(
      nchar(
        c(dimnames(yield_loss_risk)[[3]], dimnames(yield_loss_risk)[[4]],
          risk_category_names
        )
      )
    ) + 1
  ),
  create_dimvar = FALSE
)
crop_var <- ncvar_def(
  name = "NameCrop",
  units = "",
  dim = list(nchar_dim, crop_dim),
  missval = NULL,
  longname = "Crop name",
  prec = "char"
)
irrigation_var <- ncvar_def(
  name = "NameIrrigation",
  units = "",
  dim = list(nchar_dim, irrigation_dim),
  missval = NULL,
  longname = "Irrigation status",
  prec = "char"
)
risk_category_var <- ncvar_def(
  name = "NameCategory",
  units = "",
  dim = list(nchar_dim, risk_category_dim),
  missval = NULL,
  longname = "Risk category",
  prec = "char"
)
                          
# Generate output NetCDF files
Sebastian Ostberg's avatar
Sebastian Ostberg committed
827
828
if (!dir.exists("../output"))
  dir.create("../output")
829
830
831
832
# First file
print(
  paste(
    "Save gridded yield loss risk indicator variable to",
Sebastian Ostberg's avatar
Sebastian Ostberg committed
833
    "../output/yield_loss_risk.nc"
834
835
836
  )
)
nc <- nc_create(
Sebastian Ostberg's avatar
Sebastian Ostberg committed
837
  filename = "../output/yield_loss_risk.nc",
838
839
840
841
842
843
844
845
846
847
  vars = list(
    yield_loss_risk_var,
    crop_var,
    irrigation_var
  )
)
# Second file
print(
  paste(
    "Save gridded crop area at risk variable to",
Sebastian Ostberg's avatar
Sebastian Ostberg committed
848
    "../output/harvested_area_at_risk.nc"
849
850
851
  )
)
nc2 <- nc_create(
Sebastian Ostberg's avatar
Sebastian Ostberg committed
852
  filename = "../output/harvested_area_at_risk.nc",
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
  vars = list(
    harvested_area_var,
    crop_var,
    irrigation_var,
    risk_category_var
  )
)

# Write name attributes
# First file
ncvar_put(nc, crop_var, dimnames(yield_loss_risk)[[3]])
ncvar_put(nc, irrigation_var, dimnames(yield_loss_risk)[[4]])
# Second file
ncvar_put(nc2, crop_var, dimnames(yield_loss_risk)[[3]])
ncvar_put(nc2, irrigation_var, dimnames(yield_loss_risk)[[4]])
ncvar_put(nc2, risk_category_var, risk_category_names)
# Write yield_loss_risk (first file)
ncvar_put(nc, yield_loss_risk_var, c(yield_loss_risk))
# Write harvested_area_var (second file)
ncvar_put(nc2, harvested_area_var, c(ha_by_risk_category))
# Finish NetCDFs
nc_close(nc)
nc_close(nc2)
################################################################################