Skip to contents

Introduction

The two main function families in pliman —shapefile_*() and mosaic_*()— serve as the foundation for high-throughput phenotyping (HTP) analysis. The shapefile_build() function is central for constructing shapefiles, offering customizable layouts and export options that allow users to tailor shapefiles to diverse experimental designs and field layouts, enhancing the relevance and reproducibility of the data collected.

Following shapefile creation, the mosaic_*() family is pivotal for performing mosaic analysis —a process that includes calculating vegetation indices, segmenting plant canopies, and visualizing plant-level data. The mosaic_analyze() function, in particular, automates the extraction of detailed phenotypic information, making it easy to capture metrics on plant health, spatial distribution, and structure across fields. Through practical examples, users learn to leverage this function for automating plant-level data extraction, thus boosting the efficiency and depth of HTP workflows.

All shapefiles and orthomosaics used in this vignette are publicly available here, ensuring reproducibility and allowing users to follow along with each example. Since the data is accessed remotely, an active internet connection is required to run the examples. Alternatively, you can download the entire repository and modify the code to import the files from a local source, making it more convenient for offline use and adaptable to different working environments.

R packages

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(ggplot2)
library(pliman)
#> ╭ Welcome to pliman version "3.1.0"! ──────────────────────────────╮
#> │                                                                  │
#> │   Developed collaboratively by NEPEM <https://nepemufsc.com>
#> │   Group lead: Prof. Tiago Olivoto                                │
#> │   For citation, type `citation('pliman')`                        │
#> │   We welcome your feedback and suggestions!                      │
#> │                                                                  │
#> ╰────────────── Simplifying high-throughput plant phenotyping in R
#> 
#> Attaching package: 'pliman'
#> The following object is masked from 'package:dplyr':
#> 
#>     %>%

Shapefiles

Building

Shapefiles are a widely used format in geographic information systems (GIS) for representing vector data such as points, lines, and polygons. They are essential in spatial analysis and can store information about geographical features and their attributes. There are three main geometry types.

  • Points: Represent specific locations (e.g., control points, cities).
  • Lines: Represent linear features (e.g., roads, rivers).
  • Polygons: Represent areas (e.g., boundaries of regions or land parcels).

Below, we explore how theshapefile_build() function works for constructing shapefiles. By default, calling shapefile_build(mosaic, ...) allows the creation of either rectangular grids (defined by rows and columns) or custom, free-form shapes, providing flexibility for different experimental designs.

To create a rectangular grid, you define the corners of your region of interest (ROI) in this sequence: top left → top right → bottom right → bottom left → back to top left to close the polygon. If more than four points are provided, the function will interpret them as vertices of a free-form shape, constructing a customized polygon that precisely fits your specified boundaries. This flexibility makes shapefile_build() a versatile tool, allowing users to adapt their shapefile structures to match complex field layouts or specific plot arrangements.

For the purpose of this example, we will use a predefined shapefile containing control points to demonstrate the process, ensuring reproducibility and allowing you to follow along with the output.

url <- "https://github.com/TiagoOlivoto/images/raw/refs/heads/master/pliman/ortho/"
mosaic <-    mosaic_input(paste0(url, "orthosmall.tif"), info = FALSE)
cpoint <- shapefile_input(paste0(url, "controlpoints.rds"), info = FALSE)

mosaic_plot_rgb(mosaic)
shapefile_plot(cpoint, add = TRUE, lwd = 5)


# Create a basemap for further plots
bm <- mosaic_view(mosaic, r = 1, g = 2, b = 3)
#> Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
#> deprecated. It might return a CRS with a non-EPSG compliant axis order. Further
#> messages of this type will be suppressed.

shp <- shapefile_build(
  mosaic = mosaic,        # the raster file
  controlpoints = cpoint, # control points (optional)
  basemap = bm,           # basemap (optional)
  nrow = 5,               # number of rows
  ncol = 3,               # number of columns
  layout = "tbrl",        # layout definition
  serpentine = FALSE      # serpentine layout?
  
)
#>  Building the mosaic Mosaic built [13ms]
#>  Cropping the mosaic Mosaic cropped [425ms]
#>  Creating the shapes Shapes created [143ms]
#>  Finishing the shapefile Shapefile built [49ms] 

# see key aspects of the created shapefiles
shapefile_measures(shp)
#> Simple feature collection with 15 features and 11 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 734315.6 ymin: 4488975 xmax: 734327 ymax: 4488979
#> Projected CRS: WGS 72BE / UTM zone 14N
#> First 10 features:
#>    unique_id block plot_id row column   xcoord  ycoord     area perimeter width
#> 1          1   B01   P0011   1      1 734317.5 4488978 2.943949  9.127652 3.786
#> 2          2   B01   P0012   2      1 734317.5 4488978 2.943949  9.127652 3.786
#> 3          3   B01   P0013   3      1 734317.5 4488977 2.943949  9.127652 3.786
#> 4          4   B01   P0014   4      1 734317.5 4488976 2.943949  9.127652 3.786
#> 5          5   B01   P0015   5      1 734317.6 4488975 2.943949  9.127652 3.786
#> 6          6   B01   P0006   1      2 734321.2 4488978 2.943949  9.127652 3.786
#> 7          7   B01   P0007   2      2 734321.3 4488978 2.943949  9.127652 3.786
#> 8          8   B01   P0008   3      2 734321.3 4488977 2.943949  9.127652 3.786
#> 9          9   B01   P0009   4      2 734321.3 4488976 2.943949  9.127652 3.786
#> 10        10   B01   P0010   5      2 734321.3 4488975 2.943949  9.127652 3.786
#>    height                       geometry
#> 1   0.778 POLYGON ((734315.6 4488979,...
#> 2   0.778 POLYGON ((734315.6 4488978,...
#> 3   0.778 POLYGON ((734315.6 4488977,...
#> 4   0.778 POLYGON ((734315.6 4488976,...
#> 5   0.778 POLYGON ((734315.7 4488976,...
#> 6   0.778 POLYGON ((734319.3 4488979,...
#> 7   0.778 POLYGON ((734319.4 4488978,...
#> 8   0.778 POLYGON ((734319.4 4488977,...
#> 9   0.778 POLYGON ((734319.4 4488976,...
#> 10  0.778 POLYGON ((734319.4 4488976,...

The function shapefile_build() returns a Simple Feature sf object , which is a data structure used to store spatial objects (points, lines, polygons) along with associated attributes, as follows:

Attribute Description
unique_id A unique identifier for each feature (e.g., plot or area).
block The block grouping of the feature, often used in experimental designs (e.g., “B01” represents Block 1).
plot_id The identifier for the specific plot within the block (e.g., “P0001” represents Plot 1).
row The row number within the plot layout (spatial positioning of the plot within a block).
column The column number within the plot layout (spatial positioning of the plot within a block).
geometry The spatial data representing the polygon boundaries of each feature POLYGON ((x y, ...)).

shapefile_build() returns a grid layout that by default goes from left to right and top to bottom (layout = "lrtb")

bm + shapefile_view(shp, attribute = "plot_id")

By combining the layout and serpentine arguments, you can generate a total of 16 distinct layouts. The layout argument controls the primary arrangement of items, while the serpentine argument introduces an optional serpentine pattern, which alters the direction of item placement in alternating rows or columns.

The layout argument specifies the orientation of the layout and is a character string. You can choose from the following options:

  • ‘tblr’: Top to Bottom, Left to Right
  • ‘tbrl’: Top to Bottom, Right to Left
  • ‘btlr’: Bottom to Top, Left to Right
  • ‘btrl’: Bottom to Top, Right to Left
  • ‘lrtb’: Left to Right, Top to Bottom
  • ‘lrbt’: Left to Right, Bottom to Top
  • ‘rltb’: Right to Left, Top to Bottom
  • ‘rlbt’: Right to Left, Bottom to Top

The serpentine argument determines whether a serpentine layout is applied. When set to TRUE, items in alternating rows or columns will be placed in reverse order, creating a “zig-zag” pattern. By default, serpentine is set to FALSE, which means the layout follows the specified direction without altering the order in alternating rows or columns. Copy and run the following code to build a Shiny app to demonstrate how these two arguments interact and affect the layout:

library(shiny)
library(bs4Dash)
library(pliman)
library(leaflet)

# Define the UI
ui <- bs4DashPage(
  sidebar = bs4DashSidebar(disable = TRUE), 
  body = bs4DashBody(
    fluidRow(
      column(
        width = 4, # Controls will be in a 3-column layout
        title = "Controls",
        selectInput(
          inputId = "layout",
          label = "Select Layout Orientation:",
          choices = c('tblr', 'tbrl', 'btlr', 'btrl', 'lrtb', 'lrbt', 'rltb', 'rlbt'),
          selected = 'tblr'
        ),
        checkboxInput(
          inputId = "serpentine",
          label = "Apply Serpentine Layout?",
          value = FALSE
        ),
        numericInput(
          inputId = "nrow",
          label = "Number of Rows:",
          value = 5,
          min = 1,
          max = 10
        ),
        numericInput(
          inputId = "ncol",
          label = "Number of Columns:",
          value = 3,
          min = 1,
          max = 10
        ),
        numericInput(
          inputId = "pwidth",
          label = "Plot width (optional):",
          value = NULL,
          min = 0,
          max = Inf
        ),
        numericInput(
          inputId = "pheight",
          label = "Plot height (optional):",
          value = NULL,
          min = 0,
          max = Inf
        )
      ),
      column(
        width = 8, # Map plot will take 9 columns
        title = "Map",
        leafletOutput("map", height = "640px")
      )
    )
  ),
  header = bs4DashNavbar(
    title = dashboardBrand(
      title = "Live demonstration",
      color = "white",
      opacity = 0.8
    ),
    status = "white",
    fixed = TRUE
  )
)


# Define the server logic
server <- function(input, output, session) {
  url <- "https://github.com/TiagoOlivoto/images/raw/refs/heads/master/pliman/ortho/"
  mosaic <-    mosaic_input(paste0(url, "orthosmall.tif"), info = FALSE)
  cpoint <- shapefile_input(paste0(url, "controlpoints.rds"), info = FALSE)
  bm <- mosaic_view(mosaic, r = 1, g = 2, b = 3)
  # Build the shapefile
  map <- reactive({
    if(!is.na(input$pwidth) && !is.na(input$pheight)){
      pwidth <- input$pwidth
      pheight <- input$pheight
    } else{
      pwidth <- NULL
      pheight <- NULL
    }
    shp <- shapefile_build(mosaic = mosaic,
                           controlpoints = cpoint,
                           basemap = bm,
                           nrow = input$nrow,
                           ncol = input$ncol,
                           layout = input$layout,
                           serpentine = input$serpentine,
                           plot_width = pwidth,
                           plot_height = pheight,
                           verbose = FALSE)
    req(shp)
    (bm + shapefile_view(shp, attribute = "plot_id"))@map
  })
  
  output$map <- renderLeaflet({
    map()
  })
}

# Run the application 
shinyApp(ui = ui, server = server)

Exporting

When working with spatial data in R, two common file formats are .shp (shapefiles) and .rds (R serialized files). A shapefile is a standard format in Geographic Information Systems (GIS) for storing spatial vector data (points, lines, polygons) and can be imported in any GIS software, like QGIS. Despite the name, a shapefile is not a single file but a collection of related files:

  • .shp: Contains the geometry (shapes) of the features (points, lines, or polygons).
  • .cpg: Contains the character encoding used to interpret the text data in the .dbf file (defaults to UTF-8).
  • .shx: An index file that speeds up data access.
  • .dbf: Stores the attributes or properties of each feature (like ID, name, etc.).
  • .prj: Defines the Coordinate Reference System (CRS), ensuring correct placement on the Earth’s surface.

The .rds file (suggested to work with pliman) is a format specific to R, used for saving single R objects (including spatial data) in a serialized form. It’s ideal for saving R objects in their native format and later loading them back into R exactly as they were saved.

The function shapefile_export() can be used to export a shapefile created with shapefile_build() or any other SpatVector or sf object.

# export to a .rds file
shapefile_export(shp, "shape_rds.rds")
# export to a .shp file
shapefile_export(shp, "shape_shp.shp")

Importing shapefiles

You can import previously saved shapefiles using the shapefile_input() function. This function supports both .rds files and .shp files, whether they were exported using shapefile_export() or created in other software.

However, since all functions in the pliman package are designed to work with shapefiles generated within pliman, it’s crucial to ensure that specific fields—such as unique_id, block, plot_id, row, and column—are present in the shapefile. If any of these required fields are missing, unexpected errors may occur during processing.

shp <- shapefile_input("shape_rds.rds")

Exploring the mosaic_analyze() function

mosaic_analyze() is the cornerstone function in pliman for high-throughput phenotyping. It enables users to efficiently process orthomosaics and extract a wealth of data from satellite or drone imagery with just a few lines of code. In most cases, all you need is an orthomosaic (or even a .jpg image from a cellphone) and the right function parameters to unlock its full potential.

Case study

In the example below, mosaic_analyze() is used to count, measure, and extract image indices at the block, plot, and individual levels in a lettuce trial. This process is based on an orthomosaic image, as described in this paper.

A big thank you to the authors for providing the full-resolution .tif file, which enabled me to advance several functionalities in pliman, including high-throughput image analysis and data extraction at multiple levels. This kind of data sharing is invaluable for driving further innovation and tool development.

The trial was conducted using a randomized complete block design with four blocks. The researchers tested the effects of Aspergillus niger application (six different levels, combining both concentration and formulation) and three levels of phosphorus (0%, 50%, 100%) on lettuce growth.

Field Experiment
Field Experiment

In the plimans shapefile, each plot within the four blocks is represented by a unique plot_id, such as “P0001,” “P0002,” etc. These correspond to the following treatments:

Plot ID Inoculant Phosphorus (%)
P0001 NI 0
P0002 NI 50
P0003 NI 100
P0004 TS 0
P0005 TS 50
P0006 TS 100
P0007 GR2 0
P0008 GR2 50
P0009 GR2 100
P0010 GR6 0
P0011 GR6 50
P0012 GR6 100
P0013 SC2 0
P0014 SC2 50
P0015 SC2 100
P0016 SC6 0
P0017 SC6 50
P0018 SC6 100

Importing the needed files

The mosaic_input() function is used to load the mosaic of a lettuce field, and the shapefile_input() function is used to load the corresponding shapefile that delineates the plots. You can also create a shapefile with shapefile_build() (as in the previous section) or simply define the nrow and ncol arguments in mosaic_analyze().

In this example, a basemap is created using a mosaic image to serve as the foundation for further visualizations. While creating a basemap is not mandatory, it can significantly speed up the process, as functions like mosaic_analyze(), shapefile_build(), and shapefile_edit() will automatically render a leaflet map if one is not provided. By pre-creating the basemap, you avoid the overhead of rendering multiple maps, making the workflow more efficient.

Additionally, a shapefile layer is overlaid on top of the basemap to display the levels of the inoculante factor.

url <- "https://github.com/TiagoOlivoto/images/raw/refs/heads/master/pliman/lettuce/"
mos <-    mosaic_input(paste0(url, "lettuce.tif"), info = FALSE)
shp <- shapefile_input(paste0(url, "lettuce.rds"), info = FALSE)

# create a basemap
bm <- mosaic_view(mos, r = 1, g = 2, b = 3) # defaults is 1e6.. so here, a bit higher resolution is used
#>  Using `downsample = 2` to match the max_pixels
#> constraint.
bm + shapefile_view(shp, attribute = "p", color_regions = ggplot_color(3))