library(stargazer)
library(tidyverse)
library(sf)
library(terra)
library(scales)
library(knitr)
library(kableExtra)
library(raster)
library(here)
library(rnaturalearth) # for global spatial boundaries
library(bcmaps) # for BC spatial data (including boundaries)

Project Overview

This report documents a district-level land cover analysis for British Columbia (BC), Canada using Google Earth Engine (GEE) and R.

Question 1 — ESA WorldCover (2021)

1a) Data Product Description

European Space Agency WorldCover 10 m v200 is a global land-cover map with global spatial extent and a temporal extent of 2021. Its spatial resolution is 10 m, and it is an annual snapshot for the reference year 2021, meaning that its temporal resolution is annual (Van De Kerchove et al., 2022; Zanaga et al., 2021).

It was generated using a combination of Sentinel-1 synthetic aperture radar (SAR) data and Sentinel-2 multispectral optical imagery. Sentinel-1 data were acquired in Interferometric Wide (IW) mode, with an approximate spatial resolution of 5 × 20m and a temporal resolution of 12-days per satellite (NASA Earthdata, n.d.). Sentinel-2 imagery contributes data at 10m and 20m spatial resolutions, depending on the spectral band. It has a temporal resolution of ~5–6 days when both Sentinel-2 satellites are considered. These inputs were harmonized and processed to produce a global land-cover classification at 10 m spatial resolution (Van De Kerchove et al., 2022).

The overall accuracy of the ESA WorldCover product was validated by Wageningen University for statistical accuracy and IIASA for spatial accuracy, where it reached an overall accuracy of 76.7% (Tsendbazar et al., 2022). The table below includes a break down of overall accuracy per continent.

worldcover_accuracy <- tibble::tribble(
  ~Region, ~`Overall accuracy`,
  "Global", "76.7 ± 0.5",
  "Africa", "76.5 ± 1.3",
  "Europe", "77.9 ± 1.0",
  "Eurasia", "72.5 ± 1.3",
  "Asia", "82.1 ± 1.0",
  "Oceania and Australia", "72.5 ± 1.3",
  "North America", "74.6 ± 1.2",
  "South America", "77.9 ± 1.1"
)

stargazer_table_html(
  worldcover_accuracy,
  title = "Overall Accuracy for the WorldCover 2021 v200 Data Product"
)
Overall Accuracy for the WorldCover 2021 v200 Data Product
Region Overall accuracy
Global 76.7 ± 0.5
Africa 76.5 ± 1.3
Europe 77.9 ± 1.0
Eurasia 72.5 ± 1.3
Asia 82.1 ± 1.0
Oceania and Australia 72.5 ± 1.3
North America 74.6 ± 1.2
South America 77.9 ± 1.1


1b) Land Cover Categories

The ESA WorldCover product includes the following 11 classes (Van De Kerchove et al., 2022):

product_classes <- tribble(
  ~`Map Code`, ~`Land Cover Class`,
  10, "Tree cover",
  20, "Shrubland",
  30, "Grassland",
  40, "Cropland",
  50, "Built-up",
  60, "Bare / sparse vegetation",
  70, "Snow and ice",
  80, "Permanent water bodies",
  90, "Herbaceous wetland",
  95, "Mangroves",
  100, "Moss and lichen"
)

stargazer_table_html(
  product_classes,
  title = "ESA WorldCover Land Cover Classes"
)
ESA WorldCover Land Cover Classes
Map Code Land Cover Class
10 Tree cover
20 Shrubland
30 Grassland
40 Cropland
50 Built-up
60 Bare / sparse vegetation
70 Snow and ice
80 Permanent water bodies
90 Herbaceous wetland
95 Mangroves
100 Moss and lichen


1c) Fraser Valley Visualization

(i) Map Visualization with Legend

The land cover image for Fraser Valley District was exported from Google Earth Engine and visualized below with a custom legend showing the 11 ESA WorldCover classes.

# load the Fraser Valley visualized image
fv_rast <- rast("data/FV_transparent.tif")

# ESA WorldCover color scheme and labels for the legend
landcover_classes <- tribble(
  ~Value, ~Color, ~Description,
  10, "#006400", "Tree cover",
  20, "#ffbb22", "Shrubland",
  30, "#ffff4c", "Grassland",
  40, "#f096ff", "Cropland",
  50, "#fa0000", "Built-up",
  60, "#b4b4b4", "Bare / sparse vegetation",
  70, "#f0f0f0", "Snow and ice",
  80, "#0064c8", "Permanent water bodies",
  90, "#0096a0", "Herbaceous wetland",
  95, "#00cf75", "Mangroves",
  100, "#fae6a0", "Moss and lichen"
)

# drop black pixels
fv_rast[fv_rast == 0] <- NA
# create two-panel layout: map on left, legend on right
layout(matrix(c(1, 2), nrow = 1), widths = c(2, 1.5))

# plot the map
# par(mar = c(1, 1, 3, 1))
plot(fv_rast, 
     col = landcover_classes$Color,
     main = "Fraser Valley District - ESA WorldCover")

# create legend panel
# par(mar = c(1, 0, 3, 1))
plot.new()
legend(x = 0.05, y = 0.95, 
       legend = landcover_classes$Description,
       fill = landcover_classes$Color,
       border = "black",
       title = "Land Cover Classes",
       cex = 1.2,
       bty = "o")

(ii) Area Calculation by Land Cover Type in the Fraser Valley

# load worldcover CSV data
worldcover_data <- read_csv((here("data", "worldcover.csv")), show_col_types = FALSE) %>%
  rename(District = district) %>%
  dplyr::select(District, everything())

worldcover_data[is.na(worldcover_data)] <- 0

# extract Fraser Valley data
fraser_valley <- worldcover_data %>%
  filter(District == "Fraser Valley")

# Display Fraser Valley data with scrollable table
fraser_valley %>%
  kable(format = "html", 
        caption = "Land Cover Area (km²) by Type in the Fraser Valley",
        digits = 1) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = TRUE,
                position = "center",
                font_size = 11) %>%
  column_spec(1, bold = TRUE, width = "2.5cm") %>%
  scroll_box(width = "100%", height = "300px")
Land Cover Area (km²) by Type in the Fraser Valley
District Tree cover Moss and lichen Shrubland Grassland Cropland Built-up Bare / sparse vegetation Snow and ice Permanent water bodies Herbaceous wetland Mangroves
Fraser Valley 10165 338 2 1557 123 91 746 328 549 0 0

The numbers obtained through GEE’s export makes sense when compared with the visualization. Tree cover dominates other land cover types at 10,165 km\(^2\), which is evident by the large dark green area in the visualization. This also makes sense realistically given that the Fraser Valley is largely forested. Permanent water bodies, cropland, and built-up is also consistent with the visualization as well as the values extracted; we see that water covers 549 km\(^2\) of land cover in the Fraser Valley, consistent with BC’s overall abundance in hydro resources. We can visually see the distribution of cropland and built-up along the lower lip of the Fraser Valley, consistent with the urban centers and agricultural areas in the district.

1d) British Columbia - Cover Type Distribution

The CSV file for BC’s land cover data was exported from Google Earth Engine; link to the code is included in the beginning of the report.

# display data
display_data <- worldcover_data %>%
  head() # only first few rows

# customization for html output
out <- capture.output(stargazer(
  display_data,
  type = "html",
  summary = FALSE,
  rownames = FALSE,
  title = "Land Cover Area (km<sup>2</sup>) by District in British Columbia"
))
html <- paste(out, collapse = "\n")
html <- sub(
  "<table",
  "<table class='bc-equal-cols' style='margin-left:auto;margin-right:auto;'",
  html
)
knitr::asis_output(paste0(
  "<style>
    /* local: force equal column widths for this one table */
    table.bc-equal-cols { table-layout: fixed; width: 100%; font-size: 0.8em; }
    table.bc-equal-cols th, table.bc-equal-cols td {
      width: 1%;
      /* wrap between words, not within a word */
      white-space: normal;
      word-break: keep-all;
      overflow-wrap: normal;
    }
  </style>",
  "<div style='text-align:center;'>",
  html,
  "</div><br>"
))
Land Cover Area (km2) by District in British Columbia
District Tree cover Moss and lichen Shrubland Grassland Cropland Built-up Bare / sparse vegetation Snow and ice Permanent water bodies Herbaceous wetland Mangroves
Comox Valley 1356 13 0 295 5 16 11 1 41 0 0
Strathcona 13315 441 1 1363 1 12 1155 1918 444 0 0
Alberni-Clayoquot 6043 45 0 389 1 11 60 12 330 0 0
Bulkley-Nechako 60987 2126 6 8054 129 38 1528 526 4806 4 0
Capital 2085 0 0 170 3 56 4 0 52 0 0
Cariboo 60946 2317 13 12549 86 47 3250 1355 2376 11 0


Question 2 — Canada AAFC Annual Crop Inventory

2a) Data Product Description

The Canada AAFC Annual Crop Inventory is produced by Agriculture and Agri-Food Canada (AAFC)’s Earth Observation Team. The spatial extent covers most of Canada, focusing on agricultural regions. The temporal extent spans from 2009 to the present, with annual updates making the temporal resolution annual. The spatial resolution is 30m, derived from optical and radar based satellite imagery including Landsat-5, AWiFS, DMC, and RADARSAT-2 based data (Agriculture and Agri-Food Canada, 2023).

The underlying satellite data sources include: - Landsat-5: 30m spatial resolution, 16-day temporal resolution (NASA, n.d.)

  • AWiFS (Advanced Wide Field Sensor): 56m spatial resolution, 5-day revisit period (United States Geological Survey, n.d.)

  • DMC (Disaster Monitoring Constellation): 4m spatial resolution, daily revisit time (eoPortal, n.d.a)

  • RADARSAT-2: 3-100m spatial resolution, 4-day temporal resolution (eoPortal, n.d.b; Huntley et al., 2021)

The overall accuracy of this data product is at minimum 85% at the final spatial resolution of 30m (Agriculture and Agri-Food Canada, 2023).

2b) Land Cover Categories

The AAFC Annual Crop Inventory includes detailed agricultural classifications plus broader land cover classes. Land Cover categories include:

aafc_data <- read_csv(here("data", "aafc.csv"), show_col_types = FALSE) %>%
  rename(District = district) %>%
  dplyr::select(District, everything())

aafc_data[is.na(aafc_data)] <- 0

aafc_columns <- names(aafc_data)[names(aafc_data) != "District"]

aafc_categories <- tibble(
  `Land Cover Class` = aafc_columns
)

stargazer_table_html(
  aafc_categories,
  title = "AAFC Annual Crop Inventory Land Cover Classes"
)
AAFC Annual Crop Inventory Land Cover Classes
Land Cover Class
Cloud
Grassland
Agriculture (undifferentiated)
Pasture and Forages
Too Wet to be Seeded
Fallow
Cereals
Barley
Other Grains
Millet
Oats
Rye
Spelt
Triticale
Wheat
Switchgrass
Sorghum
Quinoa
Winter Wheat
Spring Wheat
Corn
Tobacco
Ginseng
Oilseeds
Borage
Camelina
Canola and Rapeseed
Flaxseed
Mustard
Safflower
Sunflower
Soybeans
Pulses
Other Pulses
Peas
Chickpeas
Beans
Fababeans
Lentils
Vegetables
Tomatoes
Potatoes
Sugarbeets
Other Vegetables
Fruits
Berries
Blueberry
Cranberry
Other Berry
Orchards
Other Fruits
Vineyards
Hops
Sod
Herbs
Nursery
Buckwheat
Canaryseed
Hemp
Vetch
Other Crops
Water
Forest (undifferentiated)
Coniferous
Broadleaf
Mixedwood
Exposed Land and Barren
Urban and Developed
Greenhouses
Shrubland
Wetland
Peatland


2c) Fraser Valley Visualization

(i) Map Visualization with Legend

The land cover image for Fraser Valley District was exported from Google Earth Engine using the AAFC Annual Crop Inventory dataset.

# Load the Fraser Valley AAFC visualized image
fv_aafc_rast <- rast("data/FV_unVisualized_AAFC.tif")

# AAFC color scheme - all land cover types present in BC with official colors, ordered by importance
aafc_classes <- tribble(
  ~Code, ~Color, ~Description,
  210, "#006600", "Coniferous",
  30, "#996666", "Exposed Land and Barren", 
  50, "#ffff00", "Shrubland",
  220, "#00cc00", "Broadleaf",
  80, "#993399", "Wetland",
  230, "#cc9900", "Mixedwood",
  110, "#cccc00", "Grassland",
  20, "#3333ff", "Water",
  34, "#cc6699", "Urban and Developed",
  122, "#ffcc33", "Pasture and Forages",
  120, "#cc6600", "Agriculture (undifferentiated)",
  153, "#d6ff70", "Canola and Rapeseed",
  133, "#dae31d", "Barley",
  146, "#92a55b", "Spring Wheat",
  162, "#8f6c3d", "Peas",
  147, "#ffff99", "Corn",
  188, "#ff6666", "Orchards",
  182, "#d20000", "Blueberry",
  131, "#ff9900", "Fallow",
  136, "#d1d52b", "Oats",
  190, "#7442bd", "Vineyards",
  183, "#cc0000", "Cranberry",
  167, "#82654a", "Beans",
  185, "#dc3200", "Other Berry",
  137, "#cacd32", "Rye",
  35, "#e1e1e1", "Greenhouses",
  145, "#809769", "Winter Wheat",
  177, "#ffcccc", "Potatoes",
  179, "#ffccff", "Other Vegetables",
  192, "#b5fb05", "Sod",
  194, "#07f98c", "Nursery",
  175, "#b74b15", "Vegetables",
  139, "#b9bc44", "Triticale",
  130, "#7899f6", "Too Wet to be Seeded",
  199, "#749a66", "Other Crops",
  189, "#c5453b", "Other Fruits",
  168, "#a39069", "Fababeans",
  197, "#8e7672", "Hemp",
  191, "#ffcc99", "Hops",
  174, "#b85900", "Lentils",
  160, "#896e43", "Pulses",
  155, "#d6cc00", "Mustard",
  154, "#8c8cff", "Flaxseed"
)

# Handle no-data values
fv_aafc_rast[fv_aafc_rast == 0] <- NA
# create two-panel layout: map on left, legend on right
layout(matrix(c(1, 2), nrow = 1), widths = c(2, 1.5))

plot(fv_aafc_rast, 
     main = "Fraser Valley District - AAFC Annual Crop Inventory",
     col = aafc_classes$Color,
     legend = FALSE)

plot.new()
par(mar = c(1, 1, 1, 1))

legend(x = 0.05, y = 0.98, 
       legend = aafc_classes$Description,
       fill = aafc_classes$Color,
       border = "black",
       title = "AAFC Land Cover Classes (All Present in BC)",
       cex = 1.1,
       bty = "o",
       xjust = 0,
       yjust = 1,
       box.lwd = 1.5,
       box.col = "black")

(ii) Area Calculation by Land Cover Type in the Fraser Valley

fraser_valley_aafc <- aafc_data %>%
  filter(District == "Fraser Valley") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  mutate(across(-District, ~ round(. / 1000000, 2)))

# Display Fraser Valley AAFC data with scrollable table
fraser_valley_aafc %>%
  kable(format = "html", 
        caption = "Land Cover Area (km²) by Type in the Fraser Valley - AAFC Data",
        digits = 2) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = TRUE,
                position = "center",
                font_size = 11) %>%
  column_spec(1, bold = TRUE, width = "2.5cm") %>%
  scroll_box(width = "100%", height = "300px")
Land Cover Area (km²) by Type in the Fraser Valley - AAFC Data
District Grassland Agriculture (undifferentiated) Pasture and Forages Fallow Barley Corn Beans Potatoes Other Vegetables Blueberry Cranberry Other Berry Orchards Other Fruits Vineyards Hops Sod Nursery Other Crops Water Coniferous Broadleaf Mixedwood Exposed Land and Barren Urban and Developed Greenhouses Shrubland Wetland
Fraser Valley 43.9 0.07 94.59 0.34 0.03 108.78 0.08 1.79 2.76 28.24 1.02 5.9 0.47 0.02 0.66 0.06 2.44 1.15 0.62 512 9343.79 787.1 619.12 1136.48 462.45 2.4 703.65 39.11

The AAFC data shows a different perspective on Fraser Valley land cover compared to ESA WorldCover. Rather than looking at land cover type, it is now separated into agricultural commodities. The data corroborates the visualization, though specifics are hard to pick out visually due to the level of detail provided by the AAFC dataset. We can see that there is much diversity in the bottom left border of the Fraser Valley.

2d) British Columbia - AAFC Cover Type Distribution

(i) Missing Data

Based on visual inspection of the data as well as through GEE’s inspector UI, several districts in British Columbia have incomplete information due to the product’s focus on agricultural areas. These include:

  • Central Coast - Minimal agricultural coverage
  • Mount Waddington - Limited to small agricultural pockets
  • Northern Rockies - Limited agricultural development
  • Skeena-Queen Charlotte - Very limited agricultural areas
  • Stikine - Minimal agricultural areas

Through GEE’s inspector UI, I also identified Kitimat-Stikine to have sparse agricultural coverage.

missing_aafc <- aafc_data %>%
  filter(`Agriculture (undifferentiated)` == 0) %>%
  filter(Potatoes == 0) %>%
  filter(`Pasture and Forages` == 0)

print(missing_aafc$District)
[1] "Central Coast"          "Mount Waddington"       "Northern Rockies"      
[4] "Skeena-Queen Charlotte" "Stikine"

(ii) Total Land Cover Area

The CSV file for BC’s AAFC land cover data was exported from Google Earth Engine; link to the code is included at the beginning of the report. Note that given the amount of categories available, I’ve filtered out any columns that are completely 0, as in not produced in BC as a whole.

# display data
aafc_nozero_data <- aafc_data %>%
  dplyr::select(where(~ any(. != 0))) %>%
  mutate(across(-District, ~ round(. / 1000000, 2)))

aafc_display_data <- aafc_nozero_data %>%
  head()

# Display first few districts with scrollable table for better readability
aafc_display_data %>%
  kable(format = "html", 
        caption = "Land Cover Area (km²) by District in British Columbia - AAFC Data (First 6 Districts)",
        digits = 2) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = TRUE,
                position = "center",
                font_size = 11) %>%
  column_spec(1, bold = TRUE, width = "2.5cm") %>%
  scroll_box(width = "100%", height = "400px")
Land Cover Area (km²) by District in British Columbia - AAFC Data (First 6 Districts)
District Grassland Agriculture (undifferentiated) Pasture and Forages Too Wet to be Seeded Fallow Barley Oats Rye Triticale Winter Wheat Spring Wheat Corn Canola and Rapeseed Flaxseed Mustard Pulses Peas Beans Fababeans Lentils Vegetables Potatoes Other Vegetables Blueberry Cranberry Other Berry Orchards Other Fruits Vineyards Hops Sod Nursery Hemp Other Crops Water Coniferous Broadleaf Mixedwood Exposed Land and Barren Urban and Developed Greenhouses Shrubland Wetland
Comox Valley 1.12 36.56 0.00 0 0 0.00 0.00 0 0 0 0 0 0 0 0 0.00 0 0 0 0 0.00 0 0 0 0 0 0 0 0 0 0 0 0 0 32.05 1203.95 130.58 0.00 37.23 129.60 0 141.59 26.35
Strathcona 275.56 1.69 0.00 0 0 0.00 0.00 0 0 0 0 0 0 0 0 0.00 0 0 0 0 0.00 0 0 0 0 0 0 0 0 0 0 0 0 0 313.29 11561.85 602.20 0.49 3513.86 64.60 0 771.12 132.91
Alberni-Clayoquot 2.42 6.87 0.00 0 0 0.00 0.00 0 0 0 0 0 0 0 0 0.00 0 0 0 0 0.00 0 0 0 0 0 0 0 0 0 0 0 0 0 287.49 5604.56 173.35 2.18 239.19 60.39 0 480.77 32.94
Bulkley-Nechako 724.40 585.20 317.66 0 0 0.00 0.37 0 0 0 0 0 0 0 0 0.00 0 0 0 0 2.39 0 0 0 0 0 0 0 0 0 0 0 0 0 4888.28 51976.44 2731.91 4786.37 4388.29 158.94 0 3700.06 2096.01
Capital 0.04 29.37 0.00 0 0 0.00 0.00 0 0 0 0 0 0 0 0 0.00 0 0 0 0 0.00 0 0 0 0 0 0 0 0 0 0 0 0 0 35.35 1790.79 95.78 0.13 12.48 303.94 0 91.04 12.45
Cariboo 2988.32 72.22 387.22 0 0 0.72 0.00 0 0 0 0 0 0 0 0 0.04 0 0 0 0 0.12 0 0 0 0 0 0 0 0 0 0 0 0 0 2284.44 56853.92 1209.73 691.85 10644.63 399.25 0 4228.83 3188.85

2e) Cropland and Grassland Comparison

To compare cropland coverage between AAFC and WorldCover, I made the assumption that any land cover that produces crops, which includes greenhouses and nurseries - with this assumption, I hypothesize that it will lead to discrepancies with the WorldCover classification as they would likely fall under Built-up.

fraser_aafc <- aafc_data %>%
  filter(District == "Fraser Valley") %>%
  dplyr::select(-Peatland, -Wetland, -Shrubland, -`Urban and Developed`, -`Exposed Land and Barren`, -`Mixedwood`, -`Broadleaf`, -`Coniferous`, -`Forest (undifferentiated)`, -Water, -Sod, -`Pasture and Forages`, -`Too Wet to be Seeded`)

fraser_valley_total_crop <- fraser_aafc %>%
  summarise(
    total_crop = sum(rowSums(
        dplyr::select(., where(is.numeric), -Grassland),
        na.rm = TRUE)))

Comparison between WorldCover and AAFC:

ESA WorldCover provides data at 10m resolution while AAFC data is at 30m resolution, however, GEE automatically converted using ee.Image.pixelArea() which outputs data in m\(^2\) units. For my ESA data output, I’ve pre-rounded to km\(^2\) using var area = ee.Number(areaDict.get('sum')).divide(1e6).round();, so to make the comparisons I will need to do the unit conversions for AAFC into \(km^2\).

The ESA WorldCover data product showed that the district of Fraser Valley houses 123 km\(^2\) of cropland and 1557 km\(^2\) of grassland. Meanwhile, AAFC’s Annual Crop Inventory shows a sum of 154 km\(^2\) of cropland and 44 km\(^2\) of Grassland. I believe these differences are due to the following factors:

  1. Cropland As hypothesized, AAFC’s cropland estimate should be larger than ESA WorldCover as I included nurseries and greenhouses as part of “cropland”’s definition, while WorldCover will likely classify them as buildings. Regardless, the estimates are quite similar, which brings confidence to the estimate.

  2. Grassland: There is a very large difference between AAFC and ESA WorldCover’s estimate on grassland, where AAFC showed only 44 km\(^2\) compared to ESA WorldCover’s 1557 km\(^2\). As a sanity check, I added up the total area for both data products for the district of Fraser Valley.

sanity_ESA <- fraser_valley %>%
  summarise(
    total_area = sum(rowSums(
        dplyr::select(., where(is.numeric)),
        na.rm = TRUE)))

sanity_aafc <- fraser_valley_aafc %>%
  summarise(
    total_area = sum(rowSums(
        dplyr::select(., where(is.numeric)),
        na.rm = TRUE)))

print(sanity_aafc$total_area)
[1] 13899.02
print(sanity_ESA$total_area)
[1] 13899

The two data products returned near identical areas, which suggests there is no unit conversion / missing data issues.

It is possible that AAFC has more detailed classifications of natural terrain, such as wetland, peatland, Exposed Land and Barren, Sod etc; most importantly, AAFC’s Annual Crop Inventory has pasture and forages as its own category which is very likely to be classified as Grassland in the ESA WorldCover data. Overall, the granularity in AAFC’s data product likely led to this difference as ESA WorldCover included many of these types of terrain together.

Question 3 — AAFC Land Cover Data Visualization

(i) Most Dominant Crop Type by District

In the AAFC data, there is a large category titled Agriculture (undifferentiated). Since I’m interested in the specific crop type, I will filter out this category while acknowleding the limitations of my visualizations. I’ve also filtered out Fallow, which is land currently resting between seedings.

dominant_crop <- aafc_data %>%
  dplyr::select(-Peatland, -Wetland, -Shrubland, -`Urban and Developed`, -`Exposed Land and Barren`, -`Mixedwood`, -`Broadleaf`, -`Coniferous`, -`Forest (undifferentiated)`, -Water, -Sod, -`Pasture and Forages`, -`Too Wet to be Seeded`, -Cloud, -Grassland, -`Agriculture (undifferentiated)`, -Fallow) 

dominant_crop$max_column <- colnames(dominant_crop)[apply(dominant_crop, 1, which.max)]

dominant_crop <- dominant_crop %>%
  dplyr::select(District, max_column)

bc_admin2 <- regional_districts()

bc_admin2 <- bc_admin2 %>%
  mutate(District_clean = str_remove_all(ADMIN_AREA_NAME, "Regional District of |Regional District| Region \\(Unincorporated\\)"),
         District_clean = str_trim(District_clean),
         District_clean = case_when(
           District_clean == "Columbia Shuswap" ~ "Columbia-Shuswap",
           District_clean == "Metro Vancouver" ~ "Greater Vancouver",
           District_clean == "qathet" ~ "Powell River",
           TRUE ~ District_clean
         )) %>%
  dplyr::filter(ADMIN_AREA_NAME != "North Coast Regional District")

bc_crops <- bc_admin2 %>%
  left_join(dominant_crop, by = c("District_clean" = "District"))

ggplot() +
  geom_sf(data = bc_crops, aes(fill = max_column), color = "black", linewidth = 0.5) +
  scale_fill_brewer(palette = "Paired") +
  theme_minimal() +
  labs(fill = "Dominant Crop Type",
       title = "Dominant Crop Types by Regional District in British Columbia")

Based on the plot of how dominant crop types vary across districts in BC, we can see that geographical elements likely dictate the type of commodity produced. Specifically looking at the Okanagan, we see dominant land use for Vineyards while canola and rapeseed dominates along the eastern border of BC. We also see corn being produced in the Fraser Valley, likely correlated with the concentration of dairy production in the region. We also see the western coast of BC favouring cereal production due to soil and moisture favourability. However, purely looking at the most produced commodity at the district level foregoes analysis of the amount of commodity output by land conditions, which I will examine in the following plot.

(ii) District Commodity Output

To better understand the scale of agricultural production across BC districts, I’ll visualize the total commodity output by district. This will show not just what is grown, but how much is produced in each area.

# Calculate total crop production by district (excluding non-crop categories)
crop_totals <- aafc_data %>%
  dplyr::select(-Peatland, -Wetland, -Shrubland, -`Urban and Developed`, 
                -`Exposed Land and Barren`, -`Mixedwood`, -`Broadleaf`, 
                -`Coniferous`, -`Forest (undifferentiated)`, -Water, -Sod, 
                -`Pasture and Forages`, -`Too Wet to be Seeded`, -Cloud, 
                -Grassland, -`Agriculture (undifferentiated)`, -Fallow) %>%
  mutate(total_crop_area = rowSums(dplyr::select(., where(is.numeric)), na.rm = TRUE)) %>%
  dplyr::select(District, total_crop_area) %>%
  mutate(total_crop_area_km2 = total_crop_area / 1000000)

bc_crop_totals <- bc_admin2 %>%
  left_join(crop_totals, by = c("District_clean" = "District"))

ggplot() +
  geom_sf(data = bc_crop_totals, aes(fill = total_crop_area_km2), linewidth = 0.3) +
  scale_fill_gradient(low = "lightblue", high = "darkgreen", name = "Total Crop Area km^2") +
  theme_minimal() +
  labs(title = "Total Cropland Area by Districts in British Columbia")

# Analyze crop diversity by calculating number of different crops grown per district
crop_diversity <- aafc_data %>%
  dplyr::select(-Peatland, -Wetland, -Shrubland, -`Urban and Developed`, 
                -`Exposed Land and Barren`, -`Mixedwood`, -`Broadleaf`, 
                -`Coniferous`, -`Forest (undifferentiated)`, -Water, -Sod, 
                -`Pasture and Forages`, -`Too Wet to be Seeded`, -Cloud, 
                -Grassland, -`Agriculture (undifferentiated)`, -Fallow) %>%
  mutate(
    crop_diversity = rowSums(dplyr::select(., where(is.numeric)) > 0, na.rm = TRUE),
    total_crop_area = rowSums(dplyr::select(., where(is.numeric)), na.rm = TRUE)
  ) %>%
  dplyr::select(District, crop_diversity, total_crop_area) %>%
  mutate(total_crop_area_km2 = total_crop_area / 1000000)

# Join with spatial data
bc_crop_diversity <- bc_admin2 %>%
  left_join(crop_diversity, by = c("District_clean" = "District"))

# Create the diversity map
ggplot() +
  geom_sf(data = bc_crop_diversity, aes(fill = crop_diversity), linewidth = 0.3) +
  scale_fill_gradient(low = "lightyellow", high = "darkred", 
                      name = "Number of Crop Types") +
  theme_minimal() +
  labs(title = "Crop Diversity by Districts in British Columbia")

Based on the commodity output analysis, the Peace River district emerges as the district with the largest crop area at over 3000 km\(^2\), followed by the Northern Rockies and Fort Nelson. This concentration in northeastern BC reflects the region’s extensive grain production, particularly canola and wheat, which benefit from the prairie-like conditions and fertile soils. The maps clearly shows that agricultural production is concentrated in specific regions - primarily the Peace River region in the northeast, the Fraser Valley in the southwest, and scattered production in the interior valleys. The coastal and mountainous regions show minimal agricultural activity. The diversity of crop types reveals that while districts such as Peace River is home to high cropland coverage, they may focus on fewer crop types while districts like the Fraser Valley and Okanagan show higher crop diversity despite smaller total areas. This reflects different agricultural strategies, such as extensive grain farming in the northeast versus intensive, diversified farming in southern valleys.

AI Use Disclaimer

AI was used in formatting the html styling and graphical coding. All code logic and interpretation are mine.

References

Agriculture and Agri-Food Canada. (2023). Annual crop inventory. https://open.canada.ca/data/en/dataset/ba2645d5-4458-414d-b196-6303ac06c1c9.
eoPortal. (n.d.a). DMC-3 (disaster monitoring constellation-3). https://www.eoportal.org/satellite-missions/dmc-3.
eoPortal. (n.d.b). RADARSAT-2. https://www.eoportal.org/satellite-missions/radarsat-2.
Huntley, D., Rotheram-Clarke, D., Pon, A., Tomaszewicz, A., Leighton, J., Cocking, R., & Joseph, J. (2021). Benchmarked RADARSAT-2, SENTINEL-1 and RADARSAT constellation mission change-detection monitoring at north slide, thompson river valley, british columbia: Ensuring a landslide-resilient national railway network. Canadian Journal of Remote Sensing, 47(4), 635–656. https://doi.org/10.1080/07038992.2021.1937968
NASA. (n.d.). Landsat 5. https://science.nasa.gov/mission/landsat-5/.
NASA Earthdata. (n.d.). Sentinel-1. https://www.earthdata.nasa.gov/data/platforms/space-based-platforms/sentinel-1.
Tsendbazar, N., Xu, P., Herold, M., Lesiv, M., & Duerauer, M. (2022). ESA WorldCover product validation report (v2.0) (WorldCover_PVR_v2.0). ESA WorldCover Consortium (VITO/Wageningen University/IIASA). https://worldcover2021.esa.int/data/docs/WorldCover_PVR_V2.0.pdf
United States Geological Survey. (n.d.). USGS EROS archive - ISRO resourcesat 1 and resourcesat 2 - AWiFS). https://www.usgs.gov/centers/eros/science/usgs-eros-archive-isro-resourcesat-1-and-resourcesat-2-awifs.
Van De Kerchove, R., Zanaga, D., Xu, P., Tsendbazar, N., & Lesiv, M. (2022). ESA WorldCover product user manual (v2.0) (WorldCover_PUM_v2.0). ESA WorldCover Consortium (VITO/Wageningen University/IIASA). https://worldcover2021.esa.int/data/docs/WorldCover_PUM_V2.0.pdf
Zanaga, D., Van De Kerchove, R., Daems, D., De Keersmaecker, W., Brockmann, C., Kirches, G., Wevers, J., Cartus, O., Santoro, M., Fritz, S., Lesiv, M., Herold, M., Tsendbazar, N. E., Xu, P., Ramoino, F., & Arino, O. (2021). ESA WorldCover 10 m 2021 v200 [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.7254221