Goal

This document is intended as an example of how to use the apsimx package to generate the inputs needed to run an APSIM (Next Gen) simulation. (With very minor changes it can be used for Classic as well). For this you will need:

Disclaimer

This document is not intended as advise or recommendation. The only goal is to illustrate the functionality in the apsimx package.

Preliminaries

Make sure that the apsimx package can find APSIM Next Gen. It should be able if installed in a (more or less) standard location. If you have problems. I suggest running the following in the R console

library(apsimx)
apsimx_options(warn.versions = FALSE) 
## This disables warnings when multiple versions are present
packageVersion("apsimx") ## apsimx package version
## [1] '2.3.1'
apsim_version() ## Check which versions of APSIM you have installed
## OS: Darwin 
## 
## 
## |APSIM           |Version.1     |
## |:---------------|:-------------|
## |Classic         |NA            |
## |Next Generation |APSIM6990.app |
wheat <- apsimx_example("Wheat") ## Check that examples can be run

Example for Meknes, Morocco

Weather data

Here I will illustrate how to construct a weather file (‘met’) for use in APSIM. However, this might be useful even if you do not work with APSIM. First, I will get data from NASA-POWER

library(apsimx) ## Load the pacakge if it was not loaded before
## After obtaining coordinates we can retrive the data
pwr1 <- get_power_apsim_met(lonlat = c(-5.5498, 33.9131), 
                            dates = c("2010-01-01", "2020-12-31"))
## The summary function might be useful for a preliminary look at the data
summary(pwr1, years = 2015:2020)
##   year months days high_maxt high_mint avg_maxt avg_mint low_maxt low_mint
## 1 2015   1:12 1:31     44.09     27.65    26.94    13.44     7.30     0.40
## 2 2016   1:12 1:31     45.11     29.64    27.23    13.55     9.30     0.24
## 3 2017   1:12 1:31     45.86     29.22    28.22    13.85     9.27    -0.16
## 4 2018   1:12 1:31     43.41     25.96    24.35    12.19     7.54     1.01
## 5 2019   1:12 1:31     43.98     25.67    26.73    12.78    10.82     1.94
## 6 2020   1:12 1:31     45.21     26.73    27.90    13.87    11.73     2.70
##   rain_sum radn_sum radn_avg
## 1   353.14  7001.12    19.18
## 2   616.68  7124.70    19.47
## 3   357.78  7272.85    19.93
## 4  1080.54  6691.01    18.33
## 5   352.97  7421.82    20.33
## 6   553.94  7219.99    19.73
## Simple graph of temperature, however hard to interpret it
plot(pwr1, years = 2015:2020)

In the previous step we obtained data from NASA-POWER, which is a gridded product, and performed a simple summary and visualization of the maximum and minimum temperature. In order to compare the different years, however, it might be easier to look at cumulative temperature and also the climatology.

## For temperature
plot(pwr1, years = 2015:2020, 
     cumulative = TRUE, climatology = TRUE)

## For precipitation
plot(pwr1, met.var = "rain", years = 2015:2020, 
     cumulative = TRUE, climatology = TRUE)

The previous graphs make it easier to see the separation for the different years. However, how about getting station data for this location? We can do this.

gsd1 <- get_gsod_apsim_met(lonlat = c(-5.5498, 33.9131), 
                           dates = c("2010-01-01", "2020-12-31"))
check_apsim_met(gsd1)
## Warning in check_apsim_met(gsd1): 38 date discontinuities found. Consider using
## napad_apsim_met.
## noname.met 
## site = station-ID 601500-99999 
## latitude = 33.879 
## longitude = -5.515 
## tav = 18.3140075376885 
## amp = 12.9108793969849 
## year day radn maxt mint rain rh windspeed 
## () () (MJ/m2/day) (oC) (oC) (mm) (%) (m/s) 
##   year day radn maxt mint rain   rh windspeed
## 1 2010   1   NA 15.0  8.0 6.10 72.3       2.9
## 2 2010   2   NA 18.7  6.4 0.00 63.1       2.5
## 3 2010   3   NA 20.1  7.5 0.00 59.7       2.6
## 4 2010   4   NA 18.5 10.6 6.10 64.0       6.8
## 5 2010   5   NA 14.2 10.0 0.51 76.0       7.9
## 6 2010   6   NA 15.4  9.5 1.02 74.1       3.3
## Warning in check_apsim_met(gsd1): Missing values found for solar radiation
## Warning in min(met[["radn"]], na.rm = TRUE): no non-missing arguments to min;
## returning Inf
## Warning in max(met[["radn"]], na.rm = TRUE): no non-missing arguments to max;
## returning -Inf
## noname.met 
## site = station-ID 601500-99999 
## latitude = 33.879 
## longitude = -5.515 
## tav = 18.3140075376885 
## amp = 12.9108793969849 
## year day radn maxt mint rain rh windspeed 
## () () (MJ/m2/day) (oC) (oC) (mm) (%) (m/s) 
##     year day radn maxt mint rain   rh windspeed
## 110 2010 111   NA 23.0 15.0   NA 70.7       3.2
## 120 2010 121   NA 24.4 14.0   NA 74.9       2.3
## 140 2010 141   NA 36.8 17.5   NA 33.3       5.5
## 254 2010 259   NA 30.0 19.9   NA 46.3       3.1
## 275 2010 282   NA 27.0 16.0   NA 66.4       2.9
## 315 2010 322   NA 14.0 10.0   NA 88.8       1.5
## Warning in check_apsim_met(gsd1): Missing values found for precipitation
## Warning in check_apsim_met(gsd1): Rain is possibly too high

The warning messages mean that there is missing data. There is no solar radiation at all and there are some gaps in the data. The function napad_apsim_met attempts to fill in these gaps.

gsd2 <- napad_apsim_met(gsd1) ## Fills in gaps 
gsd2$radn <- pwr1$radn ## replaces solar radiation using NASA-POWER
gsd3 <- impute_apsim_met(gsd2) ## Imputes missing values using linear interpolation
check_apsim_met(gsd3) ## One instance of high precipitation. Does not seem to be a problem.
## Warning in check_apsim_met(gsd3): Rain is possibly too high
## Using the plot function again
plot(gsd3, met.var = "rain",
     climatology = TRUE,
     cumulative = TRUE)

## Precipitation is not identical to NASA-POWER
## I assume it is more accurate for this location
## But I do not have other source of information to 
## cross-check. If you do, let me know.

If you want to save these file for APSIM use the write_apsim_met function.

write_apsim_met(gsd3, wrt.dir = ".", filename = "Meknes-2010-2020.met")

Soil data

Obtaining data from SoilGrids (https://soilgrids.org/) is easy. However, it might not be appropriate for point simulation and there is uncertainty associated with these estimates, but we are not taking this into account.

## Get the data
sgrd1 <- get_isric_soil_profile(lonlat = c(-5.5498, 33.9131))
## Visualize some of the data
plot(sgrd1, property = "water")

plot(sgrd1, property = "BD")

plot(sgrd1, property = "Carbon")

## We need to make changes as APSIM will not accept these values
sgrd1$soil$SAT <- sgrd1$soil$SAT * 0.9

Running an APSIM simulation

In order to run an APSIM simulation, I will use the wheat example which is distributed with the APSIM software. First, we copy that file to our current directory.

## Find where the Examples are located
ex.dir <- auto_detect_apsimx_examples()
## Copy Wheat file to current directory
file.copy(file.path(ex.dir, "Wheat.apsimx"), ".") 
## [1] TRUE

Now we inspect the weather, clock and manager.

inspect_apsimx("Wheat.apsimx", node = "Weather")
## Met file: %root%\Examples\WeatherFiles\Dalby.met
inspect_apsimx("Wheat.apsimx", node = "Clock")
## Start: 1900-01-01T00:00:00 
## End: 2000-12-31T00:00:00
inspect_apsimx("Wheat.apsimx", node = "Manager",
               parm = list("SowingRule", NA))
## Management Scripts:  SowingFertiliser Harvest SowingRule1 
## 
## Name:  SowingRule1 
## 
## 
## |parm         |value  |
## |:------------|:------|
## |StartDate    |1-may  |
## |EndDate      |10-jul |
## |MinESW       |100    |
## |MinRain      |25     |
## |RainDays     |7      |
## |CultivarName |Hartog |
## |SowingDepth  |30     |
## |RowSpacing   |250    |
## |Population   |120    |

We need to edit these options before we can run it for our example.

## Edit the path for the weather
edit_apsimx("Wheat.apsimx", 
            node = "Weather",
            value = "Meknes-2010-2020.met",
            overwrite = TRUE)
## Edited (node):  Weather 
## Edited (child):  none 
## Edited parameters:  
## New values:  Meknes-2010-2020.met 
## Created:  ./Wheat.apsimx
## Edit the Clock
edit_apsimx("Wheat.apsimx", 
            node = "Clock",
            parm = c("Start", "End"),
            value = c("2010-01-01T00:00:00", "2020-12-31T00:00:00"),
            overwrite = TRUE)
## Edited (node):  Clock 
## Edited (child):  none 
## Edited parameters:  Start End 
## New values:  2010-01-01T00:00:00 2020-12-31T00:00:00 
## Created:  ./Wheat.apsimx
## Change the sowing rule for when rain is available
edit_apsimx("Wheat.apsimx", 
            node = "Manager",
            manager.child = "SowingRule1",
            parm = "StartDate", ## This is for start date
            value = "1-oct",
            overwrite = TRUE)
## Edited (node):  Manager 
## Edited (child):  SowingRule1 
## Edited parameters:  StartDate 
## New values:  1-oct 
## Created:  ./Wheat.apsimx
edit_apsimx("Wheat.apsimx", 
            node = "Manager",
            manager.child = "SowingRule1",
            parm = "EndDate", ## This is for end date
            value = "1-nov",
            overwrite = TRUE)
## Edited (node):  Manager 
## Edited (child):  SowingRule1 
## Edited parameters:  EndDate 
## New values:  1-nov 
## Created:  ./Wheat.apsimx

Now let’s replace the soil profile.

edit_apsimx_replace_soil_profile("Wheat.apsimx",
                                 soil.profile = sgrd1,
                                 overwrite = TRUE)
## Created:  ./Wheat.apsimx

Finally we can run the Wheat example for this location.

sim <- apsimx("Wheat.apsimx")
head(sim)
##   CheckpointID SimulationID  Zone         Clock.Today
## 1            1            1 Field 2011-04-20 12:00:00
## 2            1            1 Field 2012-05-21 12:00:00
## 3            1            1 Field 2013-05-24 12:00:00
## 4            1            1 Field 2016-03-24 12:00:00
## 5            1            1 Field 2017-04-21 12:00:00
## 6            1            1 Field 2019-04-05 12:00:00
##   Wheat.Phenology.Zadok.Stage Wheat.Phenology.CurrentStageName
## 1                          90                      HarvestRipe
## 2                          90                      HarvestRipe
## 3                          90                      HarvestRipe
## 4                          90                      HarvestRipe
## 5                          90                      HarvestRipe
## 6                          90                      HarvestRipe
##   Wheat.AboveGround.Wt Wheat.AboveGround.N    Yield Wheat.Grain.Protein
## 1             1.548387          0.06193548    0.000             0.00000
## 2          1371.397666         26.68122763 5288.970            15.20962
## 3          1400.396364         15.19175661 4821.460            11.95603
## 4             1.548387          0.06193548    0.000            11.95603
## 5             1.548387          0.06193548    0.000            11.95603
## 6          1607.102436         29.34962706 5644.164            14.98801
##   Wheat.Grain.Size Wheat.Grain.Number Wheat.Grain.Total.Wt Wheat.Grain.Total.N
## 1       0.00000000               0.00               0.0000             0.00000
## 2       0.02448865           21597.64             528.8970            14.08813
## 3       0.02811418           17149.57             482.1460            10.09554
## 4       0.00000000               0.00               0.0000             0.00000
## 5       0.00000000               0.00               0.0000             0.00000
## 6       0.04055596           13916.98             564.4164            14.81519
##   Wheat.Total.Wt       Date
## 1       1.778722 2011-04-20
## 2    1570.075823 2012-05-21
## 3    1611.942972 2013-05-24
## 4       1.797956 2016-03-24
## 5       1.779879 2017-04-21
## 6    1843.013439 2019-04-05
## Visualize the data using ggplot
library(ggplot2)
ggplot(data = sim, aes(Date, Yield)) + 
  geom_point()

It looks like the conditions are not ideal for wheat growth as we do not get consistent yield.

I hope you found this tutorial useful. Let me know if you have any questions.