dev/notebook/example.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# GrowthMaps.jl example\n\n_Rafael Schouten and James Maino_\n\n\nIn this example we will calculate the expected population growth rates of\nSpotted Wing Drosophila (SWD) _D. suzukii_, for each month of the year on a 9km grid\naccross North America.\n\nWe'll experiment with running the model with a number of different datasets.\nGrowthMaps.jl, via Rasters.jl faciliates using a wide range of input sources.\nFirst we will use aggregated data in tiff files, then aggregate them to in-memory\nfiles to run the inferface. Finally if you can download it, we'll use the\n(huge) SMAP dataset, which is in a custom HDF5 format.\n\nIf you are using this as a .jmd file, it's best in\nto use atom with the \"uber-juno\" and \"language-weave\" plugins.\n\n## Load some required packages\n\nThese packages take care of loading and plotting data, and handling sci units and dates."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "using GrowthMaps, Plots, Unitful, UnitfulRecipes, Dates\nusing Rasters, ArchGDAL, NCDatasets\nusing CSV, DataFrames\nusing Unitful: °C, K, cal, mol\nbasedir = joinpath(dirname(@__FILE__), \"../\")"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Define model components\n\nFirst we'll define the growth model using `SchoolfieldIntrinsicGrowth`, based on\nSchoolfield (1981).\n\nWhen defining model components, the first parameter is a `:symbol` for the\nrequired raster layer in the source data."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "p = 3e-01\nΔH_A = 3e4cal/mol\nΔH_L = -1e5cal/mol\nΔH_H = 3e5cal/mol\nThalf_L = 2e2K\nThalf_H = 3e2K\nT_ref = K(25.0°C)\ngrowthmodel = SchoolfieldIntrinsicGrowth(p, ΔH_A, ΔH_L, Thalf_L, ΔH_H, Thalf_H, T_ref)\ngrowth = Layer(:surface_temp, K, growthmodel)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "If these are only estimated parameters, We can fit the model to a a dataset of growth\nrate and temperature. First extract our independent and dependent variables from the example CSV:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "mkpath(joinpath(basedir, \"build\"))\nobsdata = CSV.File(joinpath(basedir, \"swd_ecophys_data.csv\"), select=[:x_value, :y_value, :y_key]) |> DataFrame\nobsdata = filter(d -> d.y_key == \"r_m\", obsdata) |> dropmissing"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Then extract the required data colummns, and convert temperature values from\nunitless Celcius to explit Kelvins, using Unitful.jl:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "obsrate = obsdata.y_value\nobstemp = obsdata.x_value * °C .|> K\nobs = collect(zip(obstemp, obsrate))"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now we can fit the model. The `fitrate` function provides an easy way to fit the\nmodels in GrowthMaps or your own custom `RateModel`s, using the LsqFit.jl\npackage:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "fittedgrowth = fit(growth, obs)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now plot the fit against the data:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "temprange = (270.0:0.1:310.0)K\np = plot(x -> GrowthMaps.rate(fittedgrowth, x), temprange; label=\"fitted\")\nscatter!(p, obs; label=\"observed \")"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "We can also try tweaking the fitting the model manually in a user interface.\nModel components are immutable (for performance reasons), so we wrap the model\nin a mutable wraper so we can use the results.\nWe parametrise the model over the same temperature range that we are plotting,\nusing the :surface_temp key that the model requires:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "model = Model(fittedgrowth)\ntempdata=(surface_temp=temprange,)\nmanualfit!(model, tempdata; obs=obs)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "If you are happy with the result, you we can update extract the manual fit\nto use to generate our growth rate maps:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "fittedgrowth = parent(model)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Note that `manualfit!` will also work for a tuple of model components\nthat use the same source data, like `(growth, heatstress, coldstress)`.\n\n\n## Load spatial data\n\nLater we can use real SMAP datasets using Rasters.jl SMAPseries loader.\nBut downloading the dataset takes too long for an example.\nInstead we will download and unzip some lower resolution monthly\ndata to use in the model:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "dataurl = \"https://media.githubusercontent.com/media/cesaraustralia/GrowthMaps.jl/data/SMAP_aggregated27km.zip\"\nzipfilepath = joinpath(basedir, \"SMAP_aggregated27km.zip\")\nunzippedfolder = joinpath(basedir, \"SMAP_aggregated27km\")\nisfile(zipfilepath) || download(dataurl, zipfilepath)\nrun(`unzip -o $zipfilepath -d $basedir`);\n\n# Get the paths for to all the wilting and surface temp files\nfilenames = readdir(unzippedfolder)\nwilting_filenames = filter(fn -> occursin(r\"land_fraction_wilting\", fn), filenames)\nsurface_temp_filenames = filter(fn -> occursin(r\"surface_temp\", fn), filenames)\nwilting_paths = joinpath.(Ref(unzippedfolder), wilting_filenames)\nsurface_temp_paths = map(p -> joinpath(unzippedfolder, p), surface_temp_filenames)\n\n# Get the dates covered in the data from the `surface_temp` files list using regex\ndf = DateFormat(\"yyyymmddTHHMMSS\");\ndates = DateTime.(replace.(surface_temp_paths, Ref(r\".*_(\\d+T\\d+).tif\" => s\"\\1\")), Ref(df))"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "We know the `land_fraction_wilting` files are for the same dates.\n\nNow we have the files and date seies, we can put together a series of\nRasters.jl stacks to load lazily from disk while running `mapgrowth`.\n\nThe beauty of this approach is that we can use a lot of different source file\ntypes and folder configurations without converting them or running out of RAM."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "function makestack(i)\n    GDALstack((land_fraction_wilting=wilting_paths[i],\n               surface_temp=surface_temp_paths[i]);\n        window=(Band(1),), # lazy view applied when the file is loaded.\n        childkwargs=(usercrs=EPSG(4326),), # lets us use lat/lon without knowing the underlying projection\n    )\nend\nstacks = [makestack(i) for i in 1:length(surface_temp_paths)]\ntimedim = Ti(dates; mode=Sampled(span=Regular(Hour(3))))\ntiffseries = RasterSeries(stacks, (timedim,))"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "We can plot a layer from a file at some date in the series:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "tiffseries[Ti(Near(DateTime(2016, 2)))][:surface_temp] |> plot"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Set the time period to a month, and set the length of the subsample period to\nthe available times over one day:\n\n\n## Run a model over the spatial data\n\nDefine a timespan range to run the model over:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "tspan=DateTime(2016, 1):Month(1):DateTime(2016, 12)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Then to start, we'll run a simple model that only calculates growth rate."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "output = mapgrowth(fittedgrowth;\n  series=tiffseries,\n  tspan=tspan,\n)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Then plot the results:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "output[Ti(1)] |> plot"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "It doesn't really capture realistic population growth: there are no growth rates\nbelow zero. We need to add some stress models. Stress models model processes\nthat produce negative deathrates in the population, as oposed to rate models,\nwhere the minimum growth rate is zero.\n\nStress models have a threshold in K and mortality\nrate per degree K. `LowerStress` models stress below the given threshold, while\n`UpperStress` models stress induced above a threshold.\n\nWe will define stress models for cold, heat and wilt stress.\nAs SWD live on and around plants, we use the proportion of plants\nwilting as an indicater of stress induced by lack of moisture."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "coldthresh = 7.0°C |> K  # Enriquez2017\ncoldmort = -log(1.00) * K^-1\ncoldstress = Layer(:surface_temp, K, ColdStress(coldthresh, coldmort))\n\nheatthresh = 30.0°C |> K # Kimura2004\nheatmort = -log(1.15) * K^-1\nheatstress = Layer(:surface_temp, K, HeatStress(heatthresh, heatmort))\n\nwiltthresh = 0.5 # default?\nwiltmort = -log(1.1);\nwiltstress = Layer(:land_fraction_wilting, WiltStress(wiltthresh, wiltmort));"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "To build a more complex model, we can chain components together\nin a tuple, and again, run and plot them:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "output = mapgrowth(fittedgrowth, heatstress;\n    series=tiffseries,\n    tspan=tspan,\n)\noutput[Ti(1)] |> plot"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "That looks better. There are negative growth rates in hot regions.\nNow lets run a full model with growth and three stressors:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "output = mapgrowth(fittedgrowth, coldstress, heatstress, wiltstress;\n    series=tiffseries,\n    tspan=tspan,\n)\noutput[Ti(1)] |> plot"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Compare with observation data\n\nTo compare out simulation with observations data,\nwe'll first load them from a CSV file:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "csvurl = \"https://raw.githubusercontent.com/cesaraustralia/GrowthMaps.jl/data/Oersted_occurrence.csv\"\ncsvfilename = joinpath(basedir, \"Oersted_occurrence.csv\")\nisfile(csvfilename) || download(csvurl, csvfilename)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Then scatter them on a map:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "obs = CSV.File(csvfilename)\noccurrence = collect(zip(obs.Longitude, obs.Latitude))\np = plot(output[Ti(1)])\nscatter!(p, occurrence; markersize=2.0, markercolor=:white, markershape=:circle, label=\"obs\")"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Parametrising models using interactive maps\n\nIf you need to adjust the model based on the distribution, this\ncan be done live in the interface, as with the manual fit.\n\nBut parametrising maps on large datasets is processor intensive, which inhibits\ninteractive fedback. To reduce processing, we can aggregate the spatial data to\na more manageable size.\n\nYou can experiment with the `agg` size to compromise between quality and render time.\nLarge values will look pixelated but will run fast.\n\n`Center()` simply takes the central cell. `Statistics.mean` would take the mean value."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "agg = 8\naggseries = Rasters.aggregate(Center(), tiffseries, (Lon(agg), Lat(agg)))"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "As the combination of models is additive, we can pre-build parts of the model\nwe don't want to fit manually, which simplifies the interfaces and helps performance.\nSeeing we allready fit the growth response to empiracle data, lets just fit\nthe stress responses to the map:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "modelkwargs = (series=aggseries, tspan=tspan)\nprecomputed = mapgrowth(fittedgrowth; modelkwargs...)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Then fit the other components. `throttle` will dictate how fast the interface\nupdates. Make it larger on a slow machine, smaller on a faster one.\nThe `window` argument lets us select a window of the output to view, using\nDimensionalData.jl/Rasters.jl `Dimension`s and `Selector`s. Here we only plot\nthe first time: but you can select any/multiple, and select them with\n`Near(DateTime(2016,2))` or similar."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "model = Model(wiltstress, coldstress, heatstress)\nthrottle = 0.2\ninterface = mapfit!(model, modelkwargs;\n    occurrence=occurrence,\n    precomputed=precomputed,\n    throttle=throttle,\n    markershape=:cross,\n    markercolor=:lightblue,\n    markeropacity=0.4,\n    window=(Ti(1),),\n)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "In atom, this will show the interface in the plot pane:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "display(interface)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "To use the it in a desktop app, use Blink.jl:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "using Blink\nw = Blink.Window()\nbody!(w, interface)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "And get the updated model components from the model:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "wiltstress, coldstress, heatstress = parent(model)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Now we will put together decent population growth maps\nusing the higher resolutions data, and a monthly tiestep:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "output = mapgrowth(fittedgrowth, wiltstress, coldstress, heatstress;\n    series=tiffseries,\n    tspan=tspan,\n);\nplot(output[Ti(1:3:12)]; axis=false)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "This output is ready to use for projecting growthrates, or in a dispersal\nsimulation using Dispersal.jl. Lets save it as a NetCDF file using Rasters.jl:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "savedir = joinpath(basedir, \"build\") |> mkpath\nsavepath = joinpath(savedir, \"growthrates.ncd\")\nwrite(savepath, output)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "We can load it again with:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "rast = Raster(savepath)\nplot(rast[Ti(1)])"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Use the SMAP dataset\n\nUsing aggregated data loses some of the impact of the stress\nmodels, which respond to extremes, not averages. If you need to use\nthe model for a serious application, run it on the real SMAP dataset.\n\nUnfortunately we can't just download the data for you with a script,\nas you need to log in to an account to have access.\n\nBut Rasters.jl has a `SMAPseries()` constructor that will automate the\nwhole process of loading SMAP HDF5 files from a folder once you have\nthe data.\n\nWe can also reduce the numer of days processed to 8 each month (or whatever)\nusing the `Where` selector from DimensionalData.jl."
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "using HDF5\n# Set this to your own SMAP data folder\nsmappath = \"/home/raf/Storage/Data/SMAP/SMAP_L4_SM_gph_v4\"\n# Choose the dayes of month to use\ndays = (1, 4, 8, 12, 16, 20, 23, 27) # 8 days per month\n#days = 1:31 # all days\n# Use the `Where` selector to choose dates in the time index by day of month\nsmapseries = SMAPseries(smappath)[Ti(Where(d -> dayofmonth(d) in days))]"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "Run the model:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "@time output = mapgrowth(fittedgrowth, wiltstress, coldstress, heatstress;\n    series=smapseries, tspan=tspan,\n)\nplot(output[Ti(1:3:12)]; axis=false)"
      ],
      "metadata": {},
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "source": [
        "And save a netcdf:"
      ],
      "metadata": {}
    },
    {
      "outputs": [],
      "cell_type": "code",
      "source": [
        "write(joinpath(basedir, \"growthrates.ncd\"), output)"
      ],
      "metadata": {},
      "execution_count": null
    }
  ],
  "nbformat_minor": 2,
  "metadata": {
    "language_info": {
      "file_extension": ".jl",
      "mimetype": "application/julia",
      "name": "julia",
      "version": "1.6.3"
    },
    "kernelspec": {
      "name": "julia-1.6",
      "display_name": "Julia 1.6.3",
      "language": "julia"
    }
  },
  "nbformat": 4
}