-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathwrite_caland_inputs.r
2354 lines (2104 loc) · 138 KB
/
write_caland_inputs.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
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
379
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
428
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
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
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
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
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
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# write_caland_inputs.r
####
# California Natural and Working Lands Carbon and Greenhouse Gas
# Model (CALAND) Copyright (c) 2020, The Regents of the University of
# California, through Lawrence Berkeley National Laboratory (subject to
# receipt of any required approvals from the U.S. Dept. of Energy). All
# rights reserved.
# If you have questions about your rights to use or distribute this software,
# please contact Berkeley Lab's Intellectual Property Office at
# [email protected].
#
# NOTICE. This Software was developed under funding from the U.S. Department
# of Energy and the U.S. Government consequently retains certain rights. As
# such, the U.S. Government has been granted for itself and others acting on
# its behalf a paid-up, nonexclusive, irrevocable, worldwide license in the
# Software to reproduce, distribute copies to the public, prepare derivative
# works, and perform publicly and display publicly, and to permit others to do so.
####
# This software and its associated input data are licensed under a modified BSD open source license
# Please see license.txt for details
# Generate the `caland()` input files (.xls) based on the gis stats and parameters from the literature
######################################### Overview of `write_caland_inputs()` ##################################################
# The `write_caland_inputs()` function is defined in the write\_caland\_inputs.r file. Its purpose is to read
# the raw data files found in the caland/raw_data directory, and organize them into the detailed input files
# needed to run the `CALAND()` model (i.e., one carbon input file (.xls) and an individual scenario input file (.xls)
# for each defined scenario). Each of the input files that `write_caland_inputs()` creates are Excel workbooks, which are
# comprised of individual worksheets, one for each data table. The raw data files define the management scenarios;
# initial 2010 areas and carbon densities; annual area changes; ecosystem carbon fluxes; management parameters; parameters
# for land conversion to cultivated or developed lands; wildfire parameters; wildfire areas; mortality fractions; and climate
# change scalars. There is a suite of settings (i.e., arguments) with various options that you choose when running
# `write_caland_inputs()`, such as climate (historical, RCP 4.5, or RCP 8.5) and whether you are controlling for wildire and
# land use land cover change (LULCC) for sensitivity testing of individual practices.
# The `CALAND()` model operates on 940 land categories, but some of the raw data are not land category-specific.
# Thus, one of the main purposes of `write_caland_inputs()` is to disaggregate all non-land-category-specific raw data into
# individual land categories, and to save these expanded data tables in the carbon and scenario `CALAND()` input files. On
# the other hand, some of the raw data are land category-specific, including the initial 2010 areas and carbon densities,
# forest vegetation carbon fluxes, forest management carbon parameters, and some management areas.
# CAUTION: Creating new raw data files and using this function to create new input files for `CALAND()` is a complicated task.
# All of the raw data have to be carefully assembled to ensure it is processed as intended. Therefore, creating new input files
# is not recommended unless you have a good understanding of the required data and how the model works.
####################################### Raw data input files to `write_caland_inputs()`#######################################
# The inputs to `write_caland_inputs()`are referred to as raw data, as they are generally not region- and ownership-specific
# and must be disaggregated to all the 940 individual land categories simulated by `CALAND()`. The raw data include both .xls
# files and .csv files, which are in caland/raw\_data/. The raw data used to generate the CALAND results reported in the
# Draft California 2030 Natural and Working Lands Implementation Plan (2019) are in parentheses next to each input file below.
# Carbon parameter raw data file (e.g., lc_params.xls)
# The carbon parameter raw data file is an Excel workbook comprised of 8 individual worksheets with the following names:
# vegc\_uptake, soilc\_accum, conversion2ag\_urban, forest\_manage, dev\_manage, grass\_manage, ag\_manage, and wildfire.
# They represent vegetation carbon fluxes (historic climate); soil carbon fluxes (historic climate); land conversion C
# parameters; management C parameters for forest, developed lands, rangelands, and cultivated lands; and carbon parameters
# for low, med, and high wildfire severity, respectively.
# Column headers are on row 12
# Only non-zero values are included (missing values are zero)
# More rows can be added that define new management practices as appropriate
# Scenarios raw data file (e.g., nwl\_scenarios\_v6\_ac.xls)
# The scenarios raw data file is an Excel file comprised of at least one scenario of management practices and corresponding
# areas or fractions (one scenario per worksheet). Note the practice determines whether an area or fraction is used; it
# is not arbitrary.
# 10 columns each: Region, Land_Type, Ownership, Management, start\_year, end\_year, start\_area, end\_area, start\_frac,
# end\_frac.
# Either area (start\_area and end\_area) or fractional (start\_frac, end\_frac) values are used depending on the practice;
# the other two columns are NA. If area is required, the units can be in ha or ac, which must be specified by the
# `units_scenario` argument.
# Region and Ownership can be "All" or a specific name
# Each record defines management for a specific time period; outside of the defined time periods the management values are
# zero
# Each scenario worksheet will result in the creation of a single input scenario file (.xls), which will be comprised of
# additional worksheets of data tables derived from the other raw data files described below.
# The worksheet name for each scenario in the scenarios raw data file will determine the scenario input file name that is
# created.
# Notes on restoration practices of the following land types:
# Fresh marsh: comes out of only private and state land in the Delta, so make sure that these ownerships are designated.
# Coastal marsh comes out of only private and state land in the coastal regions, so make sure these are available, from
# cultivated
# Meadow restoration is only in the sierra cascades, and in private, state, and usfs nonwilderness, from shrub, grass,
# savanna, and woodland
# Woodland restoration comes from cultivated and grassland
# Afforestation is forest area expansion and comes from shrubland and grassland
# Reforestation is forest area expansion that comes from shrubland only to match non-regeneration of forest
# Climate raw data file (e.g., climate\_c\_scalars\_iesm\_rcp85.csv)
# The climate raw data file is a .csv file of annual climate scalars (fraction) for vegetation and soil carbon fluxes.
# The .csv file consits of four ID columns: Region, Land_Type, Ownership, and Component (Soil or Vegetation); and one
# column for each year, starting in year 2010.
# These fractions will directly scale the vegetation and soil C flux values under a historic climate in `CALAND()`
# There must be a valid climate file, even if historic climate is designated `CLIMATE = HIST`, in which case all
# climate scalars will be changed to 1.
# Wildfire raw data file (e.g., fire\_area\_canESM2\_85\_bau\_2001\_2100.csv)
# The wildfire raw data file is a .csv file of annual wildfire areas (units = ha).
# The .csv file consists of two ID columns (Region and Ownership), and a fire area column for each year, starting
# in year 2001. All possible Region-Ownership combinations are included (81 total), even if they don't exist.
# The wildfire areas from 2001 to 2015 are averaged by `write_caland_inputs()` to compute the initial 2010 wildfire areas
# written to the scenario input file. If historic climate is designated, `CLIMATE = HIST`, the average 2001 to 2015
# burn areas in the RCP 8.5 wildfire raw data are used for initial 2010 burn areas and throughout the entire simulation
# (same value, no trend).
# Wildfire severity fractions of the annual wildfire areas are hard-wired in `write_caland_inputs()`; and are not
# dependent on the selected climate. The high, medium, and low severity fractions are uniform across the State of
# California, starting with high = 0.26, medium = 0.29, and low = 0.45. The high severity fraction increases annually
# at a historical rate (i.e., 0.0027) and the low and medium fractions increase proportionally in order to account for
# the total wildfire area each year.
# Non-regeneration areas are parameterized as a function of the high severity fraction and a threshold distance from
# burn edge; however this is handled in [`CALAND()`](#arguments-in-CALAND) with the `NR_Dist` argument.
# Mortality raw data file (e.g., mortality\_annual\_july\_2018.csv)
# The mortality raw data file is a .csv file of mortality fractions (units = fraction), representing the annual
# mortality rate as a fraction of live biomass carbon (above-ground main canopy and roots).
# The .csv file consits of three ID columns (Region, Land_Type, and Ownership) and one data column (Mortality_frac).
# Land type must be specific but region and/or ownership can be "All".
# Valid land types include any woody land type that have C accumulation in vegetation (i.e., Forest, Shrubland,
# Savanna, Woodland, and Developed). However, if rows do not exist for all valid land types,
# write_caland_inputs() will assign it a default annual mortality fraction of 0.01.
# Unique treatment of Forest mortality: A doubled forest mortality trend is computed by `write_caland_inputs()`
# for 10 years (2015-2024) to emulate the recent and expected forest mortality due to insects and drought. This
# is the default setting, but can be changed in the arguments.
# Unique treatment of Developed mortality: Mortality in Developed lands is processed differently from others in
# `CALAND()` because the urban system is highly managed, and allows for more control of what happens to the dead
# biomass; the mortality is transferred to the above-ground harvest for dead removal management.
# New area GIS raw data file (e.g., CALAND_Area_Changes_2010_to_2101.csv)
# The new area GIS raw data file is a .csv file of annual area changes (sq m) for each land category, starting in
# year 2010.
# The .csv file consists of one ID column (Landcat) and one annual area change column for each year.
# With the exception of Seagrass, Fresh marsh is the only land type that does not initially exist in this GIS raw
# data files. The Fresh marsh land categories are added by `write_caland_inputs()` with initial areas of 0,
# which only increases by restoration.
# Original area GIS raw data files
# CALAND v1 and v2 used two remote-sensing based area GIS files (i.e., area_lab_sp9_own9_2001lt15_sqm_stats.csv
# and area_lab_sp9_own9_2010lt15_sqm_stats.csv), but these files were not used in the Draft California 2030
# Natural and Working Lands Climate Change Implementation Plan.
# The original area GIS files include two .csv files of areas (sq m) for each land category in 2001 and 2010,
# respectively.
# The .csv files consist of six ID columns (no headers, but rows 1 through 6 are Region integer code, Region,
# Land_type integer code, Land_type, Ownership integer code, and Ownership, respectively), and one area change
# column.
# The land category (Land_cat) integer codes will be calculated and then matched with the integer codes in the
# carbon GIS raw data files.
# The year has to be in the name of the file; currently `write_caland_inputs()` assumes that the years are 2001 and 2010.
# Carbon GIS raw data files (e.g., gss_soc_tpha_sp9_own9_2010lt15_stats.csv, lfc_agc_se_tpha_sp9_own9_2010lt15_stats.csv,
# lfc_agc_tpha_sp9_own9_2010lt15_stats.csv, lfc_bgc_se_tpha_sp9_own9_2010lt15_stats.csv,
# lfc_bgc_tpha_sp9_own9_2010lt15_stats.csv, lfc_ddc_se_tpha_sp9_own9_2010lt15_stats.csv,
# lfc_ddc_tpha_sp9_own9_2010lt15_stats.csv, lfc_dsc_se_tpha_sp9_own9_2010lt15_stats.csv,
# lfc_dsc_tpha_sp9_own9_2010lt15_stats.csv, lfc_ltc_se_tpha_sp9_own9_2010lt15_stats.csv,
# lfc_ltc_tpha_sp9_own9_2010lt15_stats.csv, lfc_usc_se_tpha_sp9_own9_2010lt15_stats.csv,
# lfc_usc_tpha_sp9_own9_2010lt15_stats.csv)
# The carbon GIS raw data files are .csv files (13 total) of univariate statistics of the initial C densities
# (units = Mg C per ha) of the seven carbon pools (i.e., soil, above-ground biomass, below-ground biomass (roots),
# down dead biomass, standing dead biomass, litter biomass, and understory biomass). These files were created by
# processing GIS raster data with a 30 m2 resolution. Thus, the statisitics are based on carbon density data
# associated with 30m2 cells within each 2010 land category boundary.
# Each .csv file consists of 14 columns in the following order: zone (i.e., land category code), label (no data),
# valid cell count (i.e., count of non-missing cells), null cell count (count of missing cells), min, max,
# range, mean, mean of absolute values (mean_of_abs), stddev (standard deviation), variance, coefficient of
# variation (coeff_var), sum, absolute sum (sum_abs).
# Only the min, max, mean, and standard deviation are used.
# Vegetation C densities are from the California Air Resources Board inventory, one value per carbon pool for
# each land category.
# Soil organic C densities are from gSSURGO; any zero values in the original data set were assumed not to be
# valid, and were therefore omitted in calculating the land category averages.
# In the soil C density raster file (i.e., the file from which the soil C raw data file is summarizing) there were
# 17601 ha (mostly rivers) with zeros, and 12236544 ha (mostly desert) with missing data. The missing values
# will be converted to zero by `write_caland_inputs()`. However, all of these zero values are excluded from
# calculating the averages.
############################################ Arguments in `write_caland_inputs()`###############################################
# 1. `c_file`: Assign a name to the carbon input file that `write_caland_inputs()` will create for `CALAND()` (needs to be .xls);
# default is `c_file = "carbon_input_nwl.xls"`.
# 2. `inputs_dir`: The directory within "./inputs" to put the generated input files; default is `c_file = ""` so they go into
# "./inputs".
# 3. `parameter_file`: Carbon parameter raw data file (needs to be .xls file); default is `parameter_file = "lc_params.xls"`.
# 4. `scenarios_file`: Scenarios raw data file with region- and/or ownership-generic management scenarios that will be expanded
# upon in the scenario files created for `CALAND()` (needs to be .xls file with one scenario per worksheet); default is
# `scenarios_file = "nwl_scenarios_v6_ac.xls"`.
# 5. `units_scenario`: Specify units for input areas in scenarios_file (must be `"ac"` or `"ha"`); default is
# `units_scenario = "ac"`.
# 6. `climate_c_file`: Climate raw data file containing either RCP 4.5 or RCP 8.5 climate scalars for vegetation and soil carbon
# fluxes (needs to be .csv); the RCP must match the RCP associated with the wildfire raw data file. If historic climate is
# desired (`CLIMATE = "HIST"`), assign either RCP 4.5 or RCP 8.5 climate raw data file, as there needs to be a valid file from
# which `write_caland_inputs()` can convert all values to 1 (i.e., no future climate effect). Default is `climate_c_file =
# "climate_c_scalars_iesm_rcp85.csv"`.
# 7. `fire_area_file`: The wildfire raw data file containing annual burned area by region-ownership (needs to be .csv); the RCP
# must match the climate raw data file; needs to be RCP 8.5 for historic climate (`CLIMATE = "HIST"`). Default is
# `fire_area_file = "fire_area_canESM2_85_bau_2001_2100.csv"`.
# 8. `mortality_file`: Mortality raw data file containing mortality rates as annual fraction of above-ground biomass; includes
# values for woody land types that have C accumulation in vegetation; any missing valid land types will be assigned a default
# annual mortality fraction of 0.01. Default is `mortality_file = "mortality_annual_july_2018.csv"`.
# 9. `area_gis_files_new`: New area GIS raw data file of GIS statistics of annual area change by land category (sq m) (needs to
# be a .csv); default is `area_gis_files_new = "CALAND_Area_Changes_2010_to_2101.csv"`.
# 10. `area_gis_files_orig`: Original area GIS raw data files of GIS statistics of areas by land category (sq m) that are compared
# to compute annual area changes (concatenate two .csv file names, one for each year); default is `area_gis_files_orig =
# c("area_lab_sp9_own9_2001lt15_sqm_stats.csv", "area_lab_sp9_own9_2010lt15_sqm_stats.csv")`.
# 11. `land_change_method`: Assigning `"Landcover"` will use the original method of remote sensing landcover change from 2001 to
# 2010 (i.e., files assigned to `area_gis_files_orig`), or assigning `"Landuse_Avg_Annual"` will use the average LUCAS-modeled
# annual area changes from 2010 to 2050 based on projected change in cultivated and developed areas (i.e., file assigned to
# `area_gis_files_new`); default is `land_change_method = "Landuse_Avg_Annual"`.
# 12. `carbon_gis_files`: Carbon GIS raw data files of GIS statistics of carbon density by carbon pool and land category
# (Mg C per ha) (concatentate 13 .csv file names); default is:
# `carbon_gis_files = c("gss_soc_tpha_sp9_own9_2010lt15_stats.csv",
# "lfc_agc_se_tpha_sp9_own9_2010lt15_stats.csv", "lfc_agc_tpha_sp9_own9_2010lt15_stats.csv",
# "lfc_bgc_se_tpha_sp9_own9_2010lt15_stats.csv", "lfc_bgc_tpha_sp9_own9_2010lt15_stats.csv",
# "lfc_ddc_se_tpha_sp9_own9_2010lt15_stats.csv", "lfc_ddc_tpha_sp9_own9_2010lt15_stats.csv",
# "lfc_dsc_se_tpha_sp9_own9_2010lt15_stats.csv", "lfc_dsc_tpha_sp9_own9_2010lt15_stats.csv",
# "lfc_ltc_se_tpha_sp9_own9_2010lt15_stats.csv", "lfc_ltc_tpha_sp9_own9_2010lt15_stats.csv",
# "lfc_usc_se_tpha_sp9_own9_2010lt15_stats.csv", "lfc_usc_tpha_sp9_own9_2010lt15_stats.csv")`.
# 13. `scen_tag`: (optional) A scenario file name tag that will be appended to the scenario names in the worksheets in the scenarios
# raw data file (`scenarios_file`); default is `scen_tag = "default"`.
# 14. `start_year`: The initial year of the simulation, which must match the year of the carbon GIS raw data files. Additionally,
# the new area file should include this year, and one of the original area GIS files needs to be for this year. Default is
# `start_year = 2010`.
# 15. `end_year`: this is the final year output desired from the `CALAND()` simulation (matches the `CALAND()` end_year argument);
# default is `end_year = 2101`.
# 16. `CLIMATE`: projected climate (RCP 4.5 or RCP 8.5) or historical climate (`"PROJ"` or `"HIST"`, respectively); affects
# wildfire areas and vegetation and soil carbon flux values; the projected climate scenario is determined by the fire and climate
# input files; `"HIST"` uses RCP 8.5 wildfire average areas from 2001 to 2015 for historical average burn area, and sets all
# climate scalars to 1. Deafult is `CLIMATE = "HIST"`.
# 17. `forest_mort_fact`: Assign a number that will be used to adjust the forest mortality during the period specified by
# `forest_mort_adj_first` and `forest_mort_adj_last`; default is `forest_mort_fact = 2`.
# 18. `forest_mort_adj_first`: the first year to adjust forest mortality; default is `forest_mort_adj_first = 2015`.
# 19. `forest_mort_adj_last`: the last year to adjust forest mortality; default is `forest_mort_adj_last = 2024`.
# 20. `control_wildfire_lulcc_file`: .csv file with flags (1 or 0) for each scenario in `scenarios_file` to control wildfire and
# LULCC. This is for sensitivity testing of individual practices or processes in `CALAND()`. Default is
# `control_wildfire_lulcc_file = "individual_proposed_sims_control_lulcc_wildfire_aug2018.csv"`. Note this will not be used
# unless `control_wildfire_lulcc = TRUE`.
# 21. `control_wildfire_lulcc`: if TRUE, read `control_wildfire_lulcc_file` to control wildfire and LULCC; default is
# `control_wildfire_lulcc = FALSE`.
############################################ Outputs from `write_caland_inputs()`############################################
# Two output files are written to caland/inputs/ (unless a sub-directory is specified differently from the default
# `inputs_dir = ""`). There is one carbon output file and one scenario output file for each scenario in the `scenarios_file`:
# (1) Carbon .xls file: `c_file`
# Carbon densities are in Mg C per ha (t C per ha)
# Factors (or scalars) are fractions
# (2) Scenario .xls file(s): <scenario\_name>\_`scen_tag`\_<climate>.xls
# One for each scenario in the `scenarios_file`
# <scenario\_name> will be named based on the worksheet names in the `scenarios_file`
# <climate> will be named based on the climate associated with the `climate_c_file` unless historical climate is chosen
# (`CLIMATE == "HIST"`), in which case "HIST" will be used.
# areas are in ha
################################################# Start script ############################################################
# this enables java to use up to 4GB of memory for reading and writing excel files
options(java.parameters = "-Xmx4g" )
# Load all the required packages
libs <- c( "XLConnect" )
for( i in libs ) {
if( !require( i, character.only=T ) ) {
cat( "Couldn't load", i, "\n" )
stop( "Use install.packages() to download this library\nOr use the GUI Package Installer\nInclude dependencies,
and install it for local user if you do not have root access\n" )
}
library( i, character.only=T )
}
write_caland_inputs <- function(scen_tag = "default", c_file = "carbon_input_nwl.xls", start_year = 2010, end_year = 2101,
CLIMATE = "PROJ", parameter_file = "lc_params.xls", scenarios_file = "nwl_scenarios_v6_ac.xls",
units_scenario = "ac",
inputs_dir = "",
climate_c_file = "climate_c_scalars_iesm_rcp85.csv",
fire_area_file = "fire_area_canESM2_85_bau_2001_2100.csv",
mortality_file = "mortality_annual_july_2018.csv",
area_gis_files_new = "CALAND_Area_Changes_2010_to_2101.csv", land_change_method = "Landuse_Avg_Annual",
area_gis_files_orig = c("area_lab_sp9_own9_2001lt15_sqm_stats.csv", "area_lab_sp9_own9_2010lt15_sqm_stats.csv"),
carbon_gis_files = c("gss_soc_tpha_sp9_own9_2010lt15_stats.csv", "lfc_agc_se_tpha_sp9_own9_2010lt15_stats.csv",
"lfc_agc_tpha_sp9_own9_2010lt15_stats.csv", "lfc_bgc_se_tpha_sp9_own9_2010lt15_stats.csv",
"lfc_bgc_tpha_sp9_own9_2010lt15_stats.csv", "lfc_ddc_se_tpha_sp9_own9_2010lt15_stats.csv",
"lfc_ddc_tpha_sp9_own9_2010lt15_stats.csv", "lfc_dsc_se_tpha_sp9_own9_2010lt15_stats.csv",
"lfc_dsc_tpha_sp9_own9_2010lt15_stats.csv", "lfc_ltc_se_tpha_sp9_own9_2010lt15_stats.csv",
"lfc_ltc_tpha_sp9_own9_2010lt15_stats.csv", "lfc_usc_se_tpha_sp9_own9_2010lt15_stats.csv",
"lfc_usc_tpha_sp9_own9_2010lt15_stats.csv"),
forest_mort_fact = 2,
forest_mort_adj_first = 2015,
forest_mort_adj_last = 2024,
control_wildfire_lulcc = FALSE,
control_wildfire_lulcc_file = "individual_proposed_sims_control_lulcc_wildfire_aug2018.csv") {
cat("Start write_caland_inputs at", date(), "\n")
num_c_in_files = length(carbon_gis_files)
# give warning if the fire_area_file does not correspond correctly with the historic CLIMATE setting
if (CLIMATE == "HIST" & fire_area_file != "fire_area_canESM2_85_bau_2001_2100.csv") {
stop("Historic climate (CLIMATE = 'HIST') must be paired with fire_area_file = 'fire_area_canESM2_85_bau_2001_2100.csv'\n")
}
# give warning if the fire_area_file does not correspond correctly with the projected CLIMATE setting
if (CLIMATE == "PROJ") {
if ( (climate_c_file == "climate_c_scalars_iesm_rcp85.csv" & fire_area_file != "fire_area_canESM2_85_bau_2001_2100.csv") |
(climate_c_file == "climate_c_scalars_iesm_rcp45.csv" & fire_area_file != "fire_area_canESM2_45_bau_2001_2100.csv") ){
stop("Projected climate and fire files must match\n")
}
}
# reference year for calculating area changes from start year
ref_year = 2001
diff_years = start_year - ref_year
scen_end_year = end_year - 1
in_dir = "raw_data/"
out_dir = paste0("inputs/", inputs_dir)
if(substr(out_dir,nchar(out_dir), nchar(out_dir)) != "/") { out_dir = paste0(out_dir, "/") }
dir.create(out_dir, recursive=TRUE)
xltag = ".xls"
# set the climate tag for the input file name
if (CLIMATE == "HIST") {
climtag = "_hist"
} else {
if (climate_c_file == "climate_c_scalars_iesm_rcp85.csv") {
climtag = "_RCP85"
} else if (climate_c_file == "climate_c_scalars_iesm_rcp45.csv") {
climtag = "_RCP45"
} else {climtag = "_climNA"}
}
c_file_out = paste0(out_dir, c_file)
# c_map_file_out = "local_files/carbon_density_map_source.xlsx"
if (land_change_method == "Landcover") {
scen_head_file = paste0(in_dir, "scenario_headers.xlsx")
} else {
scen_head_file = paste0(in_dir, "scenario_headers_usgs_area_change.xlsx")
}
param_head_file = paste0(in_dir, "parameter_headers.xlsx")
# convert acres to hectares
ac2ha = 0.404685642
# the column headers are on line 12, for both input and output files
start_row = 12
# the header goes from row 1 to row 10, and is only the first column
last_head_row = 10
# convert sq m to ha
sqm2ha = 1.0/10000
# output dataframe lists
out_scen_sheets = c("area_2010", "annual_net_area_change", "annual_managed_area", "annual_wildfire_area", "annual_mortality", "veg_climate_scalars", "soil_climate_scalars")
out_c_sheets = c("sum_allorgc_2010", "sum_biomassc_2010", "agcmain_2010", "bgcmain_2010", "usc_2010", "dsc_2010", "ddc_2010", "ltc_2010",
"soc_2010", "vegc_uptake", "soilc_accum", "conversion2ag_urban", "forest_manage", "dev_manage", "grass_manage", "ag_manage", "wildfire")
out_c_tags = c("allorgc", "biomassc", "agc", "bgc", "usc", "dsc", "ddc", "ltc", "soc", "vegc_uptake", "soilc_accum", "conversion2ag_urban", "forest_manage",
"dev_manage", "grass_manage", "ag_manage", "wildfire")
num_scen_sheets = length(out_scen_sheets)
num_c_sheets = length(out_c_sheets)
out_scen_df_list <- list()
out_c_df_list <- list()
in_c_df_list <- list()
accum_df_list <- list()
man_df_list <- list()
# regions
reg_names = c("Central_Coast", "Central_Valley", "Delta", "Deserts", "Eastside", "Klamath", "North_Coast", "Sierra_Cascades", "South_Coast")
num_reg = length(reg_names)
# land types
lt_names = c("Water", "Ice", "Barren", "Sparse", "Desert", "Shrubland", "Grassland", "Savanna", "Woodland", "Forest", "Meadow",
"Coastal_marsh", "Fresh_marsh", "Cultivated", "Developed_all", "Seagrass")
num_lt = length(lt_names)
# ownerships
own_names = c("BLM", "DoD", "Easement", "Local_gov", "NPS", "Other_fed", "Private", "State_gov", "USFS_nonwild")
num_own = length(own_names)
# useful out_c_df_list indices
allc_ind = 1
biomassc_ind = 2
cpool_start = 3
cpool_end = 9
params_start = 10
params_end = 17
vegcuptake_ind = 10
soilcaccum_ind = 11
conversion_ind = 12
forest_man_ind = 13
ag_manage_ind = 16
wildfire_ind = 17
# useful indices of the input parameter files
param_start_col = c(4, 4, 2, 5, 3, 3, 4, 2)
# some default parameters
# fresh marsh land type code
# this comes out of Delta Cultivated land only
fresh_marsh_name = "Fresh_marsh"
fresh_marsh_code = 120
delta_code = 3
cult_code = 130
# seagrass
# Seagrass estimated extent range is 4451-6070 ha (NOAA)
seagrass_reg_name = "Ocean"
seagrass_lt_name = "Seagrass"
seagrass_own_name = "Other_fed"
seagrass_reg_code = 10
seagrass_lt_code = 200
seagrass_own_code = 6
seagrass_start_area_ha = 5261.00
###### mortality
# default mortality for woody land types with vegc_uptake
# valid land types: Shrubland, Savanna, Woodland, Forest, Developed_all
# all others need to be 0, although CALAND does check this
# this default value is just used to initialize the table
mortality_default = 0.01
###### wildfire (ha)
# assume that the intensities are: High, Medium, Low
# these intensities must match those in the parameter file
# input fire area is per region-ownership, and distributed annually in caland to land types proportionally
# the initial year fire area is the annual average of 2001-2015
# BAU fire will be determined by the 2001-2015 trend of the input file, applied to the initial area
# climate driven fire area will be determined by the annual values in the provided data file, except for initial year burn area
# do not use these any more! but useful numbers for comparison
# average burned area of 2000-2015 from CALFIRE fire perimeters dataset
#wildfire_mean = 243931.10
#wildfire_stddev = 151439.00
#wildfire_ann_val = wildfire_mean
# the years to average for the initial year burn area
initial_fire_start = 2001
initial_fire_end = 2015
num_initial_fire_years = initial_fire_end - initial_fire_start + 1
# initial fractions of burn area assigned to each severity (these must sum to 1)
high_fire_frac = 0.26
low_fire_frac = 0.45
med_fire_frac = 0.29
# BAU trend in high severity fraction, change in fraction of hs burn area per year
# other two classes will decrease proportionally
hs_fire_trend = 0.0027
######## Developed_all above and below ground C density
# regional averages of cities, based on urban forest
# mcpherson et al 2017
# the reported values have been mapped to caland in the file urban_forest_data.xlsx
# the split between above and below is based on the paper assumption that 28% of total live biomass in in roots
# only the total live biomass is reported in the paper, and we were originally told the root values were not included
# we recently learned that the root values are included in the reported values
# use the above and below
# use the std err as the uncertainty (in the stddev column)
# add and subtract the stderr to get max and min
# the values are in the order of reg_names above
dev_all_agc_mean = c(14.27, 11.23, 11.23, 1.60, 7.49, 7.49, 17.32, 7.49, 6.23)
dev_all_bgc_mean = c(5.55, 4.37, 4.37, 0.62, 2.91, 2.91, 6.74, 2.91, 2.42)
dev_all_agc_stddev = c(0.33, 0.07, 0.07, 0.08, 0.11, 0.11, 0.59, 0.11, 0.08)
dev_all_bgc_stddev = c(0.01, 0.001, 0.001, 0.003, 0.003, 0.003, 0.02, 0.003, 0.002)
dev_all_agc_min = dev_all_agc_mean - dev_all_agc_stddev
dev_all_agc_max = dev_all_agc_mean + dev_all_agc_stddev
dev_all_bgc_min = dev_all_bgc_mean - dev_all_bgc_stddev
dev_all_bgc_max = dev_all_bgc_mean + dev_all_bgc_stddev
if (FALSE) { # don't need this any more
# now sum them into above only and propagate the error
dev_all_agc_mean = dev_all_agc_mean + dev_all_bgc_mean
dev_all_agc_stddev = sqrt(dev_all_agc_stddev^2 + dev_all_bgc_stddev^2)
dev_all_agc_min = dev_all_agc_mean - dev_all_agc_stddev
dev_all_agc_max = dev_all_agc_mean + dev_all_agc_stddev
dev_all_bgc_mean[] = 0
dev_all_bgc_stddev[] = 0
dev_all_bgc_min[] = 0
dev_all_bgc_max[] = 0
}
# these are the original statewide values from bjorkman et al 2015
#dev_all_agc_min = 0.13
#dev_all_agc_max = 47.01
#dev_all_agc_mean = 10.7
#dev_all_agc_stddev = 10.19
# Grassland, Savanna, Woodland soil c values
# statewide average of estimated total col soil c density from reviewed lit
# silver et al 2010
# keep the original stddev
range_soc_min = 47.0
range_soc_max = 246.0
range_soc_mean = 116.0
### forest regional npp values, all biomass (hudiburg et al 2009, appendix Table A2-A)
# units converted to MgC/ha/yr
# these are used to get relative regional breakdown to apply to the statewide input values
# north coast is from coast range
# eastside is from east cascades
# klamath is from klamath mountains
# sierra cascades is from sierra nevada
# deserts is from central basin
# south coast is from oak-chapparal
# central coast is from oak-chapparal
# Delta and central valley use the oak-chapparal because it is very close to the unweighted average of 4.76
reg_names = c("Central_Coast", "Central_Valley", "Delta", "Deserts", "Eastside", "Klamath", "North_Coast", "Sierra_Cascades", "South_Coast")
reg_vals = c(4.7, 4.7, 4.7, 2.2, 3.3, 6.0, 7.8, 4.6, 4.7)
forest_npp = data.frame(Region=reg_names, npp=reg_vals)
forest_npp$Land_Type = "Forest"
### Delta adjustments
# commenting out this section - all cultivated input values are being updated in lc_params
# Delta Cultivated peatland soil c accum values MgC/ha/yr (negative value is soil c loss) (knox 2015)
# average these? Or just use corn?
# these are consistent with Hatala 2012 and Teh 2011, but are more complete in that the harvest values are included
# corn (NEE = 2.78 MgC/ha/yr soil c loss and harvest = 2.93 MgC/ha/yr)
#soil_c_accum_peat_corn = -5.71
#soil_c_accum_peat_corn_ci = 0.34
# rice (NEE = -0.5 MgC/ha/yr (which is a soil c gain) and harvest = 1.62 MgC/ha/yr)
# methane C lost is 0.053 MgC/ha/yr
#soil_c_accum_peat_rice = -1.17
#soil_c_accum_peat_rice_ci = 1.03
# average values because not all Delta land is peatland or corn
#soil_c_accum_peat_mean = -3.44
#soil_c_accum_peat_min = -5.71
#soil_c_accum_peat_max = -1.17
#soil_c_accum_peat_stddev = 3.21
#soil_c_accum_peat_avg_ci = 0.69
# If the cultivated peatland values are used for the Delta, then the SoilCaccum_frac in ag_manage needs to be
# updated also assume the same mean benefit of 0.5 MgC/ha/yr
# for rice, this would give: 0.67/1.17=0.57
# for corn, this would give: 5.21/5.71=0.91
# for average, this would give: 2.94/3.44=0.85
#soil_c_accum_frac_peat = 0.85
#####################################################################################################################
################################# process the original lancover area files ##########################################
#####################################################################################################################
# first determine which file is the start year file
# the columns are: region code, region name, land type code, land type name, ownership code, ownership name, area
styr_ind = grep(start_year, area_gis_files_orig)
refyr_ind = grep(ref_year, area_gis_files_orig)
start_area_in = read.csv(paste0(in_dir, area_gis_files_orig[styr_ind]), header=FALSE, stringsAsFactors=FALSE)
start_area_in[,c(1,3,5)] <- as.numeric(unlist(start_area_in[,c(1,3,5)]))
ref_area_in = read.csv(paste0(in_dir, area_gis_files_orig[refyr_ind]), header=FALSE, stringsAsFactors=FALSE)
ref_area_in[,c(1,3,5)] <- as.numeric(unlist(ref_area_in[,c(1,3,5)]))
colnames(start_area_in) <- c("reg_code", "reg_name", "lt_code", "lt_name", "own_code", "own_name", "area_sqm")
colnames(ref_area_in) <- c("reg_code", "reg_name", "lt_code", "lt_name", "own_code", "own_name", "area_sqm")
# generate the land category codes and add them to the area in tables
# region*100000 + landtype*100 + ownership
# spatial and ownership can range from 0-99 classes (but zero is not used here)
# land type can range from 0 to 999
# aggregate land type (evt) is currently enumerated as multiples of 10, with water==0
start_area_in$lcat_code = start_area_in$reg_code * 100000 + start_area_in$lt_code * 100 + start_area_in$own_code
ref_area_in$lcat_code = ref_area_in $reg_code * 100000 + ref_area_in $lt_code * 100 + ref_area_in $own_code
# strip off the mismatched categories and diagnose these area differences
start_na = start_area_in[is.na(start_area_in$reg_code) & is.na(start_area_in$lt_code) & is.na(start_area_in$own_code),]
ref_na = ref_area_in[is.na(ref_area_in$reg_code) & is.na(ref_area_in$lt_code) & is.na(ref_area_in$own_code),]
cat("\nstart NA area (ha) =", start_na$area_sqm * sqm2ha, "\n")
cat("ref NA area (ha) =", ref_na$area_sqm * sqm2ha, "\n")
cat("start-ref NA area (ha) =", start_na$area_sqm * sqm2ha - ref_na$area_sqm * sqm2ha, "\n")
start_lcat_mismatch = start_area_in[is.na(start_area_in$lcat_code) & !(is.na(start_area_in$reg_code) & is.na(start_area_in$lt_code) & is.na(start_area_in$own_code)),]
ref_lcat_mismatch = ref_area_in[is.na(ref_area_in$lcat_code) & !(is.na(ref_area_in$reg_code) & is.na(ref_area_in$lt_code) & is.na(ref_area_in$own_code)),]
start_lcat_mismatch_ha = sum(start_lcat_mismatch$area_sqm * sqm2ha)
ref_lcat_mismatch_ha = sum(ref_lcat_mismatch$area_sqm * sqm2ha)
cat("\nstart lcat mismatch area (ha) =", start_lcat_mismatch_ha, "\n")
cat("ref lcat mismatch area (ha) =", ref_lcat_mismatch_ha, "\n")
cat("start-ref lcat mismatch area (ha) =", start_lcat_mismatch_ha - ref_lcat_mismatch_ha, "\n")
start_area_in = start_area_in[!is.na(start_area_in$lcat_code),]
ref_area_in = ref_area_in[!is.na(ref_area_in$lcat_code),]
# add a column for the area in ha (the values are only precise to 100 sqm)
start_area_in$start_area_ha = start_area_in$area_sqm * sqm2ha
ref_area_in$ref_area_ha = ref_area_in$area_sqm * sqm2ha
# get some necessary sums
start_area_sum_ha = sum(start_area_in$start_area_ha)
ref_area_sum_ha = sum(ref_area_in$ref_area_ha)
cat("\nstart total land area (ha) =", start_area_sum_ha, "\n")
cat("ref total land area (ha) =", ref_area_sum_ha, "\n")
cat("start-ref total land area (ha) =", start_area_sum_ha - ref_area_sum_ha, "\n")
# these are the region-ownership areas, which do not change over time
# the reference and start are nearly identical, so use the start year values
# because the normalization of area needs to be consistent with the start area
ref_area_reg_own = aggregate(ref_area_ha ~ reg_name + own_name, ref_area_in, FUN=sum)
names(ref_area_reg_own)[ncol(ref_area_reg_own)] <- "ref_reg_own_area_ha"
start_area_reg_own = aggregate(start_area_ha ~ reg_name + own_name, start_area_in, FUN=sum)
names(start_area_reg_own)[ncol(start_area_reg_own)] <- "start_reg_own_area_ha"
# calculate the annual net area change
# need to merge the data to deal with mismatches in existing categories between the years
# first subtract 2001 from 2010
# if the area is zero in 2010 and the annual loss is < 0, set the annual loss to zero, and proportionally redistribute the loss to the other land types
# this means that losses in other land types will be increased and gains in other land types will be decreased
# adjust the differences by adjusting (reducing in this case) the 2001 areas proportionally by land category within ownership class and region (so that the total areas match and total differences are zero)
# this done by adjusting the non-zero differences directly
# divide the differences by start-ref year-difference to annualize them
diff_area_in = merge(start_area_in, ref_area_in, by=c("lcat_code"), all=TRUE)
diff_area_in$start_area_ha = replace(diff_area_in$start_area_ha, is.na(diff_area_in$start_area_ha), 0)
diff_area_in$ref_area_ha = replace(diff_area_in$ref_area_ha, is.na(diff_area_in$ref_area_ha), 0)
diff_area_in$diff_area_ha = diff_area_in$start_area_ha - diff_area_in$ref_area_ha
diff_area_sum_ha = sum(diff_area_in$diff_area_ha)
# fill and drop duplicate and unnecessary columns
diff_area_in$area_sqm.x = NULL
diff_area_in$area_sqm.y = NULL
diff_area_in$reg_name.x[is.na(diff_area_in$reg_name.x)] = diff_area_in$reg_name.y[is.na(diff_area_in$reg_name.x)]
diff_area_in$lt_name.x[is.na(diff_area_in$lt_name.x)] = diff_area_in$lt_name.y[is.na(diff_area_in$lt_name.x)]
diff_area_in$own_name.x[is.na(diff_area_in$own_name.x)] = diff_area_in$own_name.y[is.na(diff_area_in$own_name.x)]
diff_area_in$reg_code.x[is.na(diff_area_in$reg_code.x)] = diff_area_in$reg_code.y[is.na(diff_area_in$reg_code.x)]
diff_area_in$lt_code.x[is.na(diff_area_in$lt_code.x)] = diff_area_in$lt_code.y[is.na(diff_area_in$lt_code.x)]
diff_area_in$own_code.x[is.na(diff_area_in$own_code.x)] = diff_area_in$own_code.y[is.na(diff_area_in$own_code.x)]
diff_area_in$reg_code.y = NULL
diff_area_in$reg_name.y = NULL
diff_area_in$lt_code.y = NULL
diff_area_in$lt_name.y = NULL
diff_area_in$own_code.y = NULL
diff_area_in$own_name.y = NULL
colnames(diff_area_in) <- c("lcat_code", "reg_code", "reg_name", "lt_code", "lt_name", "own_code", "own_name", "start_area_ha", "ref_area_ha", "diff_area_ha")
# get the region-ownership areas for the start year
diff_area_calc = merge(diff_area_in, start_area_reg_own, by=c("reg_name", "own_name"), all.x=TRUE)
diff_area_calc = diff_area_calc[order(diff_area_calc$lcat_code),]
# redistribute the unavailable losses, based on the start area propoertions
unavail_loss = diff_area_calc[diff_area_calc$diff_area_ha<0 & diff_area_calc$start_area_ha==0,]
unavail_loss_agg = aggregate(diff_area_ha ~ reg_name + own_name, unavail_loss, FUN=sum)
names(unavail_loss_agg)[ncol(unavail_loss_agg)] <- "start_reg_own_diff_area_ha"
diff_area_calc = merge(diff_area_calc, unavail_loss_agg, by=c("reg_name", "own_name"), all.x=TRUE)
diff_area_calc = diff_area_calc[order(diff_area_calc$lcat_code),]
diff_area_calc$start_reg_own_diff_area_ha = replace(diff_area_calc$start_reg_own_diff_area_ha, is.na(diff_area_calc$start_reg_own_diff_area_ha), 0)
diff_area_calc$diff_area_clean_ha = diff_area_calc$diff_area_ha
diff_area_calc$diff_area_clean_ha[diff_area_calc$diff_area_ha<0 & diff_area_calc$start_area_ha==0] = 0
diff_area_calc$diff_area_clean_ha = diff_area_calc$diff_area_clean_ha + diff_area_calc$start_reg_own_diff_area_ha * diff_area_calc$start_area_ha / diff_area_calc$start_reg_own_area_ha
diff_area_clean_sum_ha = sum(diff_area_calc$diff_area_clean_ha)
# adjust for differences due to total area mismatch
# get the region-ownership change area sums
extra_change_agg = aggregate(diff_area_clean_ha ~ reg_name + own_name, diff_area_calc, FUN=sum)
names(extra_change_agg)[ncol(extra_change_agg)] <- "extra_change_area_ha"
diff_area_calc = merge(diff_area_calc, extra_change_agg, by=c("reg_name", "own_name"), all.x=TRUE)
diff_area_calc = diff_area_calc[order(diff_area_calc$lcat_code),]
diff_area_calc$extra_change_area_ha = replace(diff_area_calc$extra_change_area_ha, is.na(diff_area_calc$extra_change_area_ha), 0)
diff_area_calc$diff_area_adj_ha = diff_area_calc$diff_area_clean_ha - diff_area_calc$extra_change_area_ha * diff_area_calc$start_area_ha / diff_area_calc$start_reg_own_area_ha
diff_area_adj_sum_ha = sum(diff_area_calc$diff_area_adj_ha)
# annualize difference
diff_area_calc$ann_diff_area_ha = diff_area_calc$diff_area_adj_ha / diff_years
ann_diff_area_sum_ha = sum(diff_area_calc$ann_diff_area_ha)
# drop the zero area cats and add fresh marsh to the Delta region
# make sure that caland checks for land cat existence for each management
# fresh marsh comes out of cultivated only, so only add the initial cultivated ownerships
# fresh marsh has been assigned a land type code of 120
zero_area = diff_area_calc[diff_area_calc$start_area_ha == 0,]
diff_area_calc = diff_area_calc[diff_area_calc$start_area_ha > 0,]
delta_cult = diff_area_calc[diff_area_calc$lt_code == cult_code & diff_area_calc$reg_code == delta_code,]
delta_cult$lcat_code = delta_cult$reg_code * 100000 + fresh_marsh_code * 100 + delta_cult$own_code
delta_cult$lt_name = fresh_marsh_name
delta_cult$lt_code = fresh_marsh_code
delta_cult[,c(8:ncol(delta_cult))] = 0.00
diff_area_calc = rbind(diff_area_calc, delta_cult)
diff_area_calc = diff_area_calc[order(diff_area_calc$lcat_code),]
# total land area
land_area_sum_ha = sum(diff_area_calc$start_area_ha)
# add seagrass
seagrass_df = diff_area_calc[1,]
seagrass_df$reg_name = seagrass_reg_name
seagrass_df$lt_name = seagrass_lt_name
seagrass_df$own_name = seagrass_own_name
seagrass_df$reg_code = seagrass_reg_code
seagrass_df$lt_code = seagrass_lt_code
seagrass_df$own_code = seagrass_own_code
seagrass_df$lcat_code = seagrass_reg_code * 100000 + seagrass_lt_code * 100 + seagrass_own_code
seagrass_df$start_area_ha = seagrass_start_area_ha
seagrass_df[,c(9:ncol(seagrass_df))] = 0.00
diff_area_calc = rbind(diff_area_calc, seagrass_df)
##############################################################################################
########################### assign area calcs to out_scen_df_list ###########################
##############################################################################################
# scen area tables
# initial area
out_scen_df_list[[1]] = data.frame(Land_Cat_ID=diff_area_calc$lcat_code, Region=diff_area_calc$reg_name, Land_Type=diff_area_calc$lt_name,
Ownership=diff_area_calc$own_name, Area_ha=diff_area_calc$start_area_ha)
# annaul area change
out_scen_df_list[[2]] = data.frame(Land_Cat_ID=diff_area_calc$lcat_code, Region=diff_area_calc$reg_name, Land_Type=diff_area_calc$lt_name,
Ownership=diff_area_calc$own_name, Area_change_ha=diff_area_calc$ann_diff_area_ha)
##############################################################################################
# if using new landuse change methods replace diff_area_calc$ann_diff_area_ha with new values
##############################################################################################
if (land_change_method == "Landuse_Avg_Annual") {
# read new area changes csv (units = m2)
new_area_changes <- read.csv(paste0(in_dir, area_gis_files_new), header=TRUE)
# add column for average annual area changes
new_area_changes$Area_change_ha <- rowMeans(new_area_changes[,2:ncol(new_area_changes)])
# convert all areas to ha
new_area_changes[,2:ncol(new_area_changes)] <- new_area_changes[,2:ncol(new_area_changes)] / 10000
# replace "m2" from end of column names with "ha"
colnames(new_area_changes)[2:(ncol(new_area_changes)-1)] <- sub("_m2", "_ha", colnames(new_area_changes[,2:(ncol(new_area_changes)-1)]))
colnames(new_area_changes)[2:(ncol(new_area_changes)-1)] <- sub("Area_", "Area_change_", colnames(new_area_changes[,2:(ncol(new_area_changes)-1)]))
# replace area change with average of individual annual area changes
# match avg area changes column with column in out_scen_df_list[[2]] by landcat
out_scen_df_list[[2]][match(new_area_changes$Landcat,out_scen_df_list[[2]][,"Land_Cat_ID"]), "Area_change_ha"] <- new_area_changes$Area_change_ha
# replace ice and water landtype area changes with 0
out_scen_df_list[[2]][out_scen_df_list[[2]]["Land_Type"] =="Ice" | out_scen_df_list[[2]][,"Land_Type"] == "Water", "Area_change_ha"] <- 0.0
}
# checks if matching is done correctly. Should equal -0.4541131. Good!
out_scen_df_list[[2]][out_scen_df_list[[2]]["Land_Cat_ID"] == 107001, "Area_change_ha"]
###### scen wildfire area table
# the available land types are forest, woodland, savanna, shrubland, grassland; but these are determined annually in caland
# caland can take multiple year columns
# assume that the intensities are: High, Medium, Low
# do not use the trends for BAU!!! They are negative!!!
fire_area_in = read.csv(paste0(in_dir,fire_area_file), stringsAsFactors = FALSE)
# first set up the table
# use all region-ownerships and all three severities so that the table is consistent across years and complete
# remove ocean and aggregate to region-ownership
burn_avail_reg_own = aggregate(Area_ha ~ Region + Ownership, out_scen_df_list[[1]][out_scen_df_list[[1]]$Region != "Ocean",], FUN=sum)
names(burn_avail_reg_own)[ncol(burn_avail_reg_own)] <- "reg_own_area_ha"
burn_avail_reg_own$lcat_code = -1
burn_avail_reg_own$lt_name = "Unspecified"
high_burn = burn_avail_reg_own
high_burn$Severity = "High"
low_burn = burn_avail_reg_own
low_burn$Severity = "Low"
med_burn = burn_avail_reg_own
med_burn$Severity = "Medium"
burn_area_ann = rbind(high_burn, low_burn, med_burn)
burn_area_ann = merge(burn_area_ann, fire_area_in, by = c("Region", "Ownership"), all.x = TRUE)
burn_area_ann = burn_area_ann[order(burn_area_ann$lcat_code, burn_area_ann$Region, burn_area_ann$Ownership, burn_area_ann$Severity),]
# get the initial year burn area by averaging
burn_area_ann$initial_area_tot = 0.0
for (y in initial_fire_start:initial_fire_end) {
# get the column index for this year
col_ind = which(names(burn_area_ann) == paste0("X", y, "_ha"))
burn_area_ann$initial_area_tot = burn_area_ann$initial_area_tot + burn_area_ann[,col_ind] / num_initial_fire_years
} # end for y loop to set initial burn area
burn_area_ann[is.na(burn_area_ann)] <- 0.0
# calc the temporal trend slopes
# fits[1,] are the intercepts by row
# fits[2,] are the slopes by row
#start_col = which(names(burn_area_ann) == paste0("X", initial_fire_start, "_ha"))
#end_col = which(names(burn_area_ann) == paste0("X", initial_fire_end, "_ha"))
#fits = apply(burn_area_ann[,start_col:end_col], 1, function(x) lm(x~c(initial_fire_start:initial_fire_end), na.action=na.exclude)$coefficients)
# set initial year burn severity proportions and column name
burn_area_ann$initial_area_sev = 0.0
burn_area_ann[burn_area_ann$Severity == "High", "initial_area_sev"] = high_fire_frac * burn_area_ann[burn_area_ann$Severity == "High", "initial_area_tot"]
burn_area_ann[burn_area_ann$Severity == "Low", "initial_area_sev"] = low_fire_frac * burn_area_ann[burn_area_ann$Severity == "Low", "initial_area_tot"]
burn_area_ann[burn_area_ann$Severity == "Medium", "initial_area_sev"] = med_fire_frac * burn_area_ann[burn_area_ann$Severity == "Medium", "initial_area_tot"]
# start building the output table
out_scen_df_list[[4]] = data.frame(Land_Cat_ID=burn_area_ann$lcat_code, Region=burn_area_ann$Region, Land_Type=burn_area_ann$lt_name, Ownership=burn_area_ann$Ownership, Severity=burn_area_ann$Severity, start_ha=burn_area_ann$initial_area_sev)
names(out_scen_df_list[[4]])[ncol(out_scen_df_list[[4]])] <- paste0(start_year,"_ha")
# loop through the remaining years
for (y in (start_year+1):(end_year-1)) {
# get the column index for this year
col_ind = which(names(burn_area_ann) == paste0("X", y, "_ha"))
# adjust the fire severity proportions based on BAU; make sure they add to 1
hff_cur = high_fire_frac + (hs_fire_trend * (y-start_year))
lff_cur = low_fire_frac - low_fire_frac / (low_fire_frac + med_fire_frac) * (hs_fire_trend * (y-start_year))
mff_cur = med_fire_frac - low_fire_frac / (low_fire_frac + med_fire_frac) * (hs_fire_trend * (y-start_year))
a = (1 - hff_cur) / (lff_cur + mff_cur)
lff_cur = a * lff_cur
mff_cur = a * mff_cur
# if projected climate, then assign the annual years
# else if BAU/historical then apply initial year burn area
if (CLIMATE != "HIST") {
# projected: set burn severity proportions and column name
burn_area_ann[burn_area_ann$Severity == "High", col_ind] = hff_cur * burn_area_ann[burn_area_ann$Severity == "High", col_ind]
burn_area_ann[burn_area_ann$Severity == "Low", col_ind] = lff_cur * burn_area_ann[burn_area_ann$Severity == "Low", col_ind]
burn_area_ann[burn_area_ann$Severity == "Medium", col_ind] = mff_cur * burn_area_ann[burn_area_ann$Severity == "Medium", col_ind]
} else {
# BAU/historical: first set the annual area, then burn severity
burn_area_ann[,col_ind] = burn_area_ann$initial_area_tot
burn_area_ann[burn_area_ann$Severity == "High", col_ind] = hff_cur * burn_area_ann[burn_area_ann$Severity == "High", col_ind]
burn_area_ann[burn_area_ann$Severity == "Low", col_ind] = lff_cur * burn_area_ann[burn_area_ann$Severity == "Low", col_ind]
burn_area_ann[burn_area_ann$Severity == "Medium", col_ind] = mff_cur * burn_area_ann[burn_area_ann$Severity == "Medium", col_ind]
}
# set negative values to zero
burn_area_ann[,-c(1:5)][burn_area_ann[,-c(1:5)] < 0] == 0
# store this year
out_scen_df_list[[4]] = cbind(out_scen_df_list[[4]], burn_area_ann[,col_ind])
names(out_scen_df_list[[4]])[ncol(out_scen_df_list[[4]])] <- paste0(y,"_ha")
} # end for y loop for remaining years
###### scen mortality table
# this is a multiple year table
# currently, the only land types with explicit mortality are:
# shrubland, savanna, woodland, forest, developed_all
# this is because other land types do not have above gound carbon accumulation (i.e. mortality is implicit)
# read in the mortality values
# all valid land types must be explicitly specified
mort_in = read.csv(paste0(in_dir,mortality_file), stringsAsFactors = FALSE)
# initial year
mortality_types = out_scen_df_list[[1]][out_scen_df_list[[1]]$Land_Type == "Forest" | out_scen_df_list[[1]]$Land_Type == "Woodland" | out_scen_df_list[[1]]$Land_Type == "Savanna" | out_scen_df_list[[1]]$Land_Type == "Shrubland" | out_scen_df_list[[1]]$Land_Type == "Developed_all",]
names(mortality_types)[ncol(mortality_types)] <- "start_frac"
year_col_start = ncol(mortality_types)
mortality_types$start_frac = mortality_default
# create table of the mortality values by location specified
mort_in_landcat = mort_in[mort_in$Region != "All" & mort_in$Ownership != "All",]
mort_in_allreg = mort_in[mort_in$Region == "All" & mort_in$Ownership != "All",]
mort_in_allown = mort_in[mort_in$Region != "All" & mort_in$Ownership == "All",]
mort_in_allregown = mort_in[mort_in$Region == "All" & mort_in$Ownership == "All",]
# land cat
if (nrow(mort_in_landcat) > 0) {
mortality_types = merge(mortality_types, mort_in_landcat, by = c("Region", "Land_Type", "Ownership"), all.x = TRUE)
mortality_types$start_frac[!is.na(mortality_types$Mortality_frac)] = mortality_types$Mortality_frac[!is.na(mortality_types$Mortality_frac)]
mortality_types$Mortality_frac = NULL
}
# all region
if (nrow(mort_in_allreg) > 0) {
mortality_types = merge(mortality_types, mort_in_allreg[,c("Land_Type", "Ownership", "Mortality_frac")], by = c("Land_Type", "Ownership"), all.x = TRUE)
mortality_types$start_frac[!is.na(mortality_types$Mortality_frac)] = mortality_types$Mortality_frac[!is.na(mortality_types$Mortality_frac)]
mortality_types$Mortality_frac = NULL
}
# all ownership
if (nrow(mort_in_allown) > 0) {
mortality_types = merge(mortality_types, mort_in_allown[,c("Region", "Land_Type", "Mortality_frac")], by = c("Region", "Land_Type"), all.x = TRUE)
mortality_types$start_frac[!is.na(mortality_types$Mortality_frac)] = mortality_types$Mortality_frac[!is.na(mortality_types$Mortality_frac)]
mortality_types$Mortality_frac = NULL
}
# all reigon and all ownership
if (nrow(mort_in_allregown) > 0) {
mortality_types = merge(mortality_types, mort_in_allregown[,c("Land_Type", "Mortality_frac")], by = c("Land_Type"), all.x = TRUE)
mortality_types$start_frac[!is.na(mortality_types$Mortality_frac)] = mortality_types$Mortality_frac[!is.na(mortality_types$Mortality_frac)]
mortality_types$Mortality_frac = NULL
}
mortality_types = mortality_types[,c("Land_Cat_ID", "Region", "Land_Type", "Ownership", "start_frac")]
if(forest_mort_fact != 1) {
# increased forest mortality
# need to add the years surround the jumps so that each period has a flat mortality rate
mortality_types$year2_frac = mortality_types$start_frac
mortality_types$year3_frac = mortality_types$start_frac
mortality_types$year3_frac[mortality_types$Land_Type == "Forest"] = forest_mort_fact * mortality_types$year3_frac[mortality_types$Land_Type == "Forest"]
mortality_types$year3_frac[mortality_types$Land_Type == "Forest" & mortality_types$year3_frac >= 1] = 1
mortality_types$year4_frac = mortality_types$year3_frac
mortality_types$year5_frac = mortality_types$start_frac
mortality_types = mortality_types[order(mortality_types$Land_Cat_ID),]
names(mortality_types)[year_col_start+1] <- paste0(forest_mort_adj_first-1,"_frac")
names(mortality_types)[year_col_start+2] <- paste0(forest_mort_adj_first,"_frac")
names(mortality_types)[year_col_start+3] <- paste0(forest_mort_adj_last,"_frac")
names(mortality_types)[year_col_start+4] <- paste0(forest_mort_adj_last+1,"_frac")
}
names(mortality_types)[year_col_start] <- paste0(start_year,"_frac")
out_scen_df_list[[5]] = mortality_types
###### scen climate scalars tables
# these values are the direct scalars
# vegetation first, then soil
climate_c_in = read.csv(paste0(in_dir,climate_c_file), stringsAsFactors = FALSE)
clim_start_col = which(names(climate_c_in) == paste0("X", start_year))
clim_end_col = ncol(climate_c_in)
# veg
# start building the output table
out_scen_df_list[[6]] = out_scen_df_list[[1]][,c(1:4)]
climate_c_veg = out_scen_df_list[[1]][,c(1:4)]
climate_c_veg = merge(climate_c_veg, climate_c_in[climate_c_in$Component == "Vegetation", c(1:3,clim_start_col:clim_end_col)], by = c("Region", "Land_Type", "Ownership"), all.x = TRUE)
climate_c_veg = climate_c_veg[order(climate_c_veg$Land_Cat_ID),]
climate_c_veg$Component = NULL
climate_c_veg[is.na(climate_c_veg)] = 1
if(CLIMATE == "HIST") { climate_c_veg[,-c(1:4)] = 1 }
# add the year data to the table
out_scen_df_list[[6]] = cbind(out_scen_df_list[[6]], climate_c_veg[,c(5:ncol(climate_c_veg))])
# set the column names
for (i in c(5:ncol(out_scen_df_list[[6]]))) {
names(out_scen_df_list[[6]])[i] <- paste0(start_year+i-5)
}
# soil
# start building the output table
out_scen_df_list[[7]] = out_scen_df_list[[1]][,c(1:4)]
climate_c_soil = out_scen_df_list[[1]][,c(1:4)]
climate_c_soil = merge(climate_c_soil, climate_c_in[climate_c_in$Component == "Soil", c(1:3,clim_start_col:clim_end_col)], by = c("Region", "Land_Type", "Ownership"), all.x = TRUE)
climate_c_soil = climate_c_soil[order(climate_c_soil$Land_Cat_ID),]
climate_c_soil$Component = NULL
climate_c_soil[is.na(climate_c_soil)] = 1
if(CLIMATE == "HIST") { climate_c_soil[,-c(1:4)] = 1 }
# add the year data to the table
out_scen_df_list[[7]] = cbind(out_scen_df_list[[7]], climate_c_soil[,c(5:ncol(climate_c_soil))])
# set the column names
for (i in c(5:ncol(out_scen_df_list[[7]]))) {
names(out_scen_df_list[[7]])[i] <- paste0(start_year+i-5)
}
###### read the scenario definition file
scenin_wrkbk = loadWorkbook(paste0(in_dir,scenarios_file))
# worksheet/table names
scenin_sheets = getSheets(scenin_wrkbk)
num_scenin_sheets = length(scenin_sheets)
# NA values need to be converted to numeric
# the warnings thrown by readWorksheet below are ok because they just state that the NA string can't be converted a number so it is
# converted to NA value
scenin_df_list <- list()
for (i in 1:num_scenin_sheets) {
scenin_df_list[[i]] <- readWorksheet(scenin_wrkbk, i, startRow = start_row, colTypes = c(rep("character",4), rep("numeric",50)), forceConversion = TRUE)
# convert management acres to hectares as needed
if (units_scenario=="ac") {
scenin_df_list[[i]][,c("start_area", "end_area")] = scenin_df_list[[i]][,c("start_area", "end_area")] * ac2ha
}
} # end for i loop to read in scenarios
###### read the scenario headers file for the outputs
scenhead_wrkbk = loadWorkbook(scen_head_file)
# worksheet/table names
scenhead_sheets = getSheets(scenhead_wrkbk)
num_scenhead_sheets = length(scenhead_sheets)
scen_head_df_list <- list()
for (i in 1:num_scenhead_sheets) {
scen_head_df_list[[i]] <- readWorksheet(scenhead_wrkbk, i, startRow = 1, header = FALSE)
}
###### read the parameter file
param_wrkbk = loadWorkbook(paste0(in_dir,parameter_file))
# worksheet/table names
param_sheets = getSheets(param_wrkbk)
num_param_sheets = length(param_sheets)
# NA values need to be converted to numeric
# the warnings thrown by readWorksheet below are ok because they just state that the NA string can't be converted a number so it is
# converted to NA value
c_col_types1 = c("character", rep("numeric",50))
c_col_types2 = c("character", "character", rep("numeric",50))
c_col_types3 = c("character", "character", "character", rep("numeric",50))
c_col_types4 = c("character", "character", "character", "character", rep("numeric",50))
# Load the param worksheets into a list of data frames
param_df_list <- list()
param_head_list <- list()
for (i in 1:2) { # vegc_uptake and soilc_accum
param_head_list[[i]] <- readWorksheet(param_wrkbk, i, startRow = 1, endRow = last_head_row, header=FALSE)
param_df_list[[i]] <- readWorksheet(param_wrkbk, i, startRow = start_row, colTypes = c_col_types3, forceConversion = TRUE)
}
for (i in 3:3) { # conversion2ag_urban
param_head_list[[i]] <- readWorksheet(param_wrkbk, i, startRow = 1, endRow = last_head_row, header=FALSE)
param_df_list[[i]] <- readWorksheet(param_wrkbk, i, startRow = start_row, colTypes = c_col_types1, forceConversion = TRUE)
}
# forest_manage
i = 4
param_head_list[[i]] <- readWorksheet(param_wrkbk, i, startRow = 1, endRow = last_head_row, header=FALSE)
param_df_list[[i]] <- readWorksheet(param_wrkbk, i, startRow = start_row, colTypes = c_col_types4, forceConversion = TRUE)
for (i in 5:6) { # dev_manage to grass_manage
param_head_list[[i]] <- readWorksheet(param_wrkbk, i, startRow = 1, endRow = last_head_row, header=FALSE)
param_df_list[[i]] <- readWorksheet(param_wrkbk, i, startRow = start_row, colTypes = c_col_types2, forceConversion = TRUE)
}
# ag_manage
i = 7
param_head_list[[i]] <- readWorksheet(param_wrkbk, i, startRow = 1, endRow = last_head_row, header=FALSE)
param_df_list[[i]] <- readWorksheet(param_wrkbk, i, startRow = start_row, colTypes = c_col_types3, forceConversion = TRUE)
# wildfire
i=8
param_head_list[[i]] <- readWorksheet(param_wrkbk, i, startRow = 1, endRow = last_head_row, header=FALSE)
param_df_list[[i]] <- readWorksheet(param_wrkbk, i, startRow = start_row, colTypes = c_col_types2, forceConversion = TRUE)
###### read the parameter headers file for the outputs
paramhead_wrkbk = loadWorkbook(param_head_file)
# worksheet/table names
paramhead_sheets = getSheets(paramhead_wrkbk)
num_paramhead_sheets = length(paramhead_sheets)
param_head_df_list <- list()
for (i in 1:num_paramhead_sheets) {
param_head_df_list[[i]] <- readWorksheet(paramhead_wrkbk, i, startRow = 1, header = FALSE)
}
if (control_wildfire_lulcc) {
# read in the instruction file for input scenario
control_wildfire_lulcc_df <-read.csv(paste0(in_dir,control_wildfire_lulcc_file), header = TRUE)
}