diff --git a/src/analysis/makieplots.jl b/src/analysis/makieplots.jl
index d03c32eb2d7db33d1074f136e3d53581cf871033..4e53f381f679d1e54fd82375d44a627e179d5efe 100644
--- a/src/analysis/makieplots.jl
+++ b/src/analysis/makieplots.jl
@@ -3,6 +3,8 @@
 ### This file visualises model output using Makie.
 ###
 
+const MIN_AXISLIMIT::Float64 = 1e-5
+
 """
     visualisemap(model, date, landcover)
 
@@ -59,7 +61,7 @@ function populationtrends(model::SimulationModel)
     update_theme!(palette=(color=cgrad(:seaborn_bright, ncolors),), cycle=[:color])
     f = Figure()
     dates = @param(core.startdate):@param(core.enddate)
-    axlimits = (1, length(dates), 0, maximum(pops[!,:Abundance])*1.1)
+    axlimits = (1, length(dates), 0, max(MIN_AXISLIMIT, maximum(pops[!,:Abundance]; init=0)*1.1))
     ax = Axis(f[1,1], xlabel="Date", ylabel="Population size",
               limits=axlimits, xticks=datetickmarks(dates))
     for s in @param(nature.targetspecies)
@@ -82,7 +84,7 @@ function skylarkpopulation(model::SimulationModel)
     pops = @data("skylark_abundance")
     f = Figure()
     dates = @param(core.startdate):@param(core.enddate)
-    axlimits = (1, length(dates), 0, maximum(pops[!,:TotalAbundance]))
+    axlimits = (1, length(dates), 0, max(MIN_AXISLIMIT, maximum(pops[!,:TotalAbundance]; init=0)))
     ax = Axis(f[1,1], xlabel="Date", ylabel="Population size",
               limits=axlimits, xticks=datetickmarks(dates))
     lines!(f[1,1], Vector{Float32}(pops.TotalAbundance), linewidth=3, label="Total abundance")
@@ -108,20 +110,28 @@ function skylarkstats(model::SimulationModel)
     # nesting habitat
     habitattype = unique(nestingdata.Landcover)
     htfrequency = [count(x -> x==h, nestingdata.Landcover) for h in habitattype]
-    barplot(f[1,1], 1:length(habitattype), htfrequency,
-             axis = (xticks=(1:length(habitattype), habitattype), title="Nesting habitat"))
+    if ! isempty(habitattype)
+        barplot(f[1,1], 1:length(habitattype), htfrequency,
+                axis = (xticks=(1:length(habitattype), habitattype), title="Nesting habitat"))
+    end
     croptype = unique(nestingdata.Crop)
     ctfrequency = [count(x -> x==c, nestingdata.Crop) for c in croptype]
-    barplot(f[2,1], 1:length(croptype), ctfrequency,
-             axis = (xticks=(1:length(croptype), croptype), title="Nesting crop type"))
+    if ! isempty(croptype)
+        barplot(f[2,1], 1:length(croptype), ctfrequency,
+                axis = (xticks=(1:length(croptype), croptype), title="Nesting crop type"))
+    end
     # mortality
     causes = unique(mortalitydata.Cause)
     mtfrequency = [count(x -> x==c, mortalitydata.Cause) for c in causes]
-    barplot(f[1,2], 1:length(causes), mtfrequency,
-            axis = (xticks=(1:length(causes), causes), title="Causes of Mortality"))
+    if ! isempty(causes)
+        barplot(f[1,2], 1:length(causes), mtfrequency,
+                axis = (xticks=(1:length(causes), causes), title="Causes of Mortality"))
+    end
     # territory size
-    hist(f[2,2], nestingdata.Territory, bins=20,
-         axis=(ylabel="Frequency", title="Territory size (ha)"))
+    if ! isempty(nestingdata.Territory)
+        hist(f[2,2], nestingdata.Territory, bins=20,
+             axis=(ylabel="Frequency", title="Territory size (ha)"))
+    end
     f
 end
 
@@ -136,17 +146,19 @@ function croptrends(model::SimulationModel)
     f = Figure(size=(1200,1000))
     dates = @param(core.startdate):@param(core.enddate)
     # plot cropped area over time
-    axlimitsarea = (1, length(dates), 0, maximum(cropdata[!,:Area])*1.1)
+    axlimitsarea = (1, length(dates), 0, max(MIN_AXISLIMIT, maximum(cropdata[!,:Area]; init=0)*1.1))
     ax1 = Axis(f[1,1], xlabel="Date", ylabel="Cropped area (ha)",
-              limits=axlimitsarea, xticks=datetickmarks(dates))
+               limits=axlimitsarea, xticks=datetickmarks(dates))
     for c in unique(cropdata.Crop)
         points = @select!(@subset(cropdata, :Crop .== c), :Area)
         (iszero(size(points)[1]) || iszero(sum(points.Area))) && continue
         lines!(f[1,1], Vector{Float32}(points.Area), linewidth=3, label=c)
     end
-    axislegend("Crop"; position=:rt)
+    if ! isempty(cropdata.Crop)
+        axislegend("Crop"; position=:rt)
+    end
     # plot average crop height over time
-    axlimitsheight = (1, length(dates), 0, maximum(cropdata[!,:Height])*1.1)
+    axlimitsheight = (1, length(dates), 0, max(MIN_AXISLIMIT, maximum(cropdata[!,:Height]; init=0)*1.1))
     ax2 = Axis(f[2,1], xlabel="Date", ylabel="Average plant height (cm)",
               limits=axlimitsheight, xticks=datetickmarks(dates))
     for c in unique(cropdata.Crop)
@@ -154,7 +166,9 @@ function croptrends(model::SimulationModel)
         (iszero(size(points)[1]) || iszero(sum(points.Height))) && c != "no growth" && continue
         lines!(f[2,1], Vector{Float32}(points.Height), linewidth=3, label=c)
     end
-    axislegend("Crop"; position=:rt)
+    if ! isempty(cropdata.Crop)
+        axislegend("Crop"; position=:rt)
+    end
     f
 end