Introduction

This is a skeleton document illustration several type of plot. Maps of Leicestershire are included along with a set of data for the Green Carpet moth in the same county.

The data format is fussy, but other templates are in the pipeline to assist with processing spreadsheet data into the format used here.

Your text can be inserted anywhere between the code block and formatted with markdown.

See https://www.markdownguide.org/

Markdown is designed to be simple and is worth mastering becuase it makes it possible to automate the production of complex documents. These days authors of books often use markdown because of its simplicity.

Code is inserted into the special blocks and can be run in stages within R Studio. Code can also be adjusted by hand. Interactive help is coming soon.

My aim is to automate as much as possible within the template.

Code setup

# Set code chunk options for all chunks
# These can be overridden at the chunk level, but setting global options 
# ensures consistency of chunk behaviour.
# To print a version of this document without code set echo = FALSE. 
#Include false will ignore the code chunk completely!
knitr::opts_chunk$set(include = TRUE, 
                      echo = FALSE, 
                      warning=FALSE, 
                      message=FALSE, 
                      error=FALSE, 
                      results='asis',
                      out.width='\\textwidth') 
# Globally set the figure width here in inches.
# Load some libraries that are needed but not yet transferred into the package.
library(sp)
library(raster)
library(ggmap)
library(tidyverse)
library(broom)
library(ggpolypath)
library(PrepareDataForETL)
library(reshape2)
library(fossil)

Setup map

data("leics") # This is already part of the package so no path needed.
# This code makes a polygon stencil to mask the edge of the map.
# Masking is needed to hide edge effects when desity distributions
# are used
temp.bb <- unlist(attr(leics,"bb"))
# Note make sure the coordinates are in the correct order
temp.coords <- cbind(
  temp.bb[c(1,3,3,1)],
  temp.bb[c(2,2,4,4)]
)

temp.sp <- SpatialPolygons(
  list(Polygons(list(Polygon(temp.coords)), "id")), 
  proj4string = CRS(proj4string(leics)))
temp.sp_diff <- erase(temp.sp, leics)
sp_diff_df <- fortify(temp.sp_diff) 
# This is a stencil with a leicestershire shape hole.

# Try this test plot
my.map <- ggplot() + 
   coord_map(projection = "mercator") + # Set a mercator projection
  geom_polypath(
         aes( long, lat, group=group),
         data = sp_diff_df,
         fill="red",
         alpha=1
       ) 


# Plot the map

my.map

# Clean up
rm(list = ls(pattern = "temp*"))

Data Preparation

The data has been cleaned using the PrepareDataForETL package written by author PJP. More templates are on the way to help convert Excel spreadsheets and other data into the required format for analysis.

Records for Green Carpet Moth

 # These are already part of the package so no path needed.
data("Green.Carpet.Data")

# This is already part of the package so no path needed.
data("leics") 

# This package also saves images in an image subdirectory
dir.create("images", showWarnings = FALSE)

The data comprises of r nrow(Green.Carpet.Data) observations. Several calculated fields have been added:

\newpage

Distribution all records

# A trial plot of all data
# Note that as the data preparation is now complete, 
# this code is all about visualisation.
data.to.plot <- Green.Carpet.Data
my.map <- ggplot(data = leics) + 
  geom_path( aes(x = long, y = lat, group = group)) +
  coord_map(projection = "mercator") + 
  geom_point(data = data.to.plot, aes(x = Long, y = Lat, colour = "black")) +
  theme(axis.line=element_blank(),
      axis.text.x=element_blank(),
      axis.text.y=element_blank(),
      axis.ticks=element_blank(),
      axis.title.x=element_blank(),
      axis.title.y=element_blank(),
      legend.position="none",
      panel.background=element_blank(),
      panel.border=element_blank(),
      panel.grid.major=element_blank(),
      panel.grid.minor=element_blank(),
      plot.background=element_blank())

my.map

# This line saves the plot
ggsave("images/all_leics.png", device = "png")

# Clean up

rm(my.map, data.to.plot)
rm(list = ls(pattern = "temp*"))

This map shows all r nrow(Green.Carpet.Data) observations marked on the map of Leicestershire with underlying geographical areas marked.

\newpage

Phenology plot of all records using ggplot

 temp.my.plot <- 
  ggplot(Green.Carpet.Data,aes(x=as.numeric(Green.Carpet.Data$DOY) ) ) + 
  geom_density(fill="Red", bw = 4) +
   labs(title="Phenology",x="Day of Year", y = "Density") +
  # Uncomment these lines to add vertical lines through the brood peaks
  #geom_vline(aes(xintercept=Brood.peaks[1]), color="blue") + 
  #geom_vline(aes(xintercept=Brood.peaks[2]), color="blue") +
  #geom_vline(aes(xintercept=Brood.peaks[3]), color="blue") +
  #geom_vline(aes(xintercept=Brood.peaks[4]), color="blue") +
   xlim(1, 365) # Set the x-axis limits.

temp.my.plot

# This line saves the plot
ggsave("images/Phenology-density.png", device = "png")

# Clean up
rm(list = ls(pattern = "temp*"))

Phenology by Year

This plot uses a heatmap type approach coupled with filtering

# Now to try an create a summary.
# Confine the date range from 1960

temp.filter.year <- 1959
temp.data.to.plot <- Green.Carpet.Data[Green.Carpet.Data$YYYY > temp.filter.year,]

temp.data.to.plot$Observations <- 1 
# Dummy numeric place holder for observations. One row equals one observation.


# Now coerce to numeric types to plot nicely.
temp.data.to.plot$YYYY <- temp.data.to.plot$YYYY %>% as.integer()
temp.data.to.plot$DOY <- temp.data.to.plot$DOY %>% as.integer()



# We are now ready to make a heatmap using geom_tile to stack the slices.
# The x scale is DOY so to plot all the data use 1, 365.
# The y-scale is year,
# If data is outside of these ranges you will get warning along side your plots.
temp.my.plot <- ggplot(data = temp.data.to.plot ) + 
   scale_x_continuous(limits = c(1, 365), expand = c(0, 0)) +
  scale_y_continuous(limits = c(1960, 2015), expand = c(0, 0)) + 

  labs(fill='Observations') +
  geom_tile(aes(x = DOY, 
                        y = YYYY, 
                        group = DOY, 
                        fill= "red")) +
  ggtitle("Phenology chart by DOY number and year") +

  # geom_vline(aes(xintercept=Brood.peaks[1]), color="blue") +
  # geom_vline(aes(xintercept=Brood.peaks[2]), color="blue") +
  # geom_vline(aes(xintercept=Brood.peaks[3]), color="blue") +
  # geom_vline(aes(xintercept=Brood.peaks[4]), color="blue") +
 #xlim(1,52,  expand = c(0, 0)) + ylim(1960, 2015,  expand = c(0, 0)) +
  theme(
          legend.position="none",
          plot.background=element_blank()
          )

temp.my.plot

ggsave("images/Phenology-heatmap-kmeans-doy.png", device = "png")

# Clean up
rm(list = ls(pattern = "temp*"))

\newpage

Latitute plot

temp.data.to.plot <- Green.Carpet.Data
temp.data.to.plot$Observations <- 1 # Numeric place holder for observations


# Now coerce to types to play nicely.
temp.data.to.plot$YYYY <- temp.data.to.plot$YYYY %>% as.integer()
temp.data.to.plot$DOY <- temp.data.to.plot$DOY %>% as.integer()
temp.data.to.plot$Lat <- temp.data.to.plot$Lat %>% as.double()
temp.data.to.plot$Long <- temp.data.to.plot$Long %>% as.double()

# We are now ready to make a heatmap.

temp.my.plot <- ggplot() + 
  geom_point(data = temp.data.to.plot, aes(x = Lat, 
                                           y = YYYY
                                           )) +

  geom_tile(alpha = 0.8) +
 # scale_fill_gradient(low = "grey", high= "black") +
   ggtitle("Lattitude Chart") +

  theme(
          legend.position="right",
          plot.background=element_blank()
          )

temp.my.plot

ggsave("images/Lattitude-chart.png", device = "png")

# Clean up
rm(list = ls(pattern = "temp*"))

\newpage

Plot by Longitude

temp.data.to.plot <- Green.Carpet.Data
temp.data.to.plot$Observations <- 1 # Numeric place holder for observations


# Now coerce to types to play nicely.
temp.data.to.plot$YYYY <- temp.data.to.plot$YYYY %>% as.integer()
temp.data.to.plot$DOY <- temp.data.to.plot$DOY %>% as.integer()
temp.data.to.plot$Lat <- temp.data.to.plot$Lat %>% as.double()
temp.data.to.plot$Long <- temp.data.to.plot$Long %>% as.double()

# We are now ready to make a heatmap.

temp.my.plot <- ggplot() + 
  geom_point(data = temp.data.to.plot, aes(x = Long, 
                                           y = YYYY
                                           )) +

  geom_tile(alpha = 0.8) +
 # scale_fill_gradient(low = "grey", high= "black") +
   ggtitle("Longitude Chart") +

  theme(
          legend.position="right",
          plot.background=element_blank()
          )

temp.my.plot

ggsave("images/Longitude-chart.png", device = "png")

# Clean up
rm(list = ls(pattern = "temp*"))

\newpage

Density distribution

data.to.plot <- Green.Carpet.Data

    # Build the map in layers
     my.map <- ggplot() +  ggtitle("Density distribution") + 
       geom_path(data =leics, aes(x = long, y = lat, group = group)) + 

       # Observation density
       stat_density2d(  aes( x = Long, 
                             y = Lat, 
                             fill = ..level.., 
                             alpha = ..level.. ,
                            colour = "black" ), 
                            size = .5, bins = 8, 
                            data = data.to.plot, 
                            geom = "polygon", 
                            contour = TRUE, 
                            show.legend = FALSE, 
                           inherit.aes = FALSE) +
      geom_point(data = data.to.plot, 
                 aes(x = Long, y = Lat,  alpha=.8)) +
      coord_map(projection = "mercator") + 
          theme(axis.line=element_blank(),
          axis.text.x=element_blank(),
          axis.text.y=element_blank(),
          axis.ticks=element_blank(),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          legend.position="none",
          panel.background=element_blank(),
          panel.border=element_blank(),
          panel.grid.major=element_blank(),
          panel.grid.minor=element_blank(),
          plot.background=element_blank()) +
        geom_polypath(
         aes( long, lat, group=group),
         sp_diff_df,
         fill="white",
         alpha=1
       )  


   print(my.map)

# Clean up
     rm(data.to.plot, my.map, area.name)
     rm(list = ls(pattern = "temp*"))


enpjp/PrepareDataForETL documentation built on Oct. 3, 2020, 3:01 a.m.