R/plot_distmap.R

In rSymbiota: Interface to 'Symbiota' Portals for Accessing Multi-Organismal Biodiversity Data

### Code adapted from karthik/gbifmap2.R:
### https://gist.github.com/karthik/4254076

#' Plot distribution data on map
#'
#' @param x an object of class "\code{records}", see \link{symbiota}
#' @param mapdatabase The map database to use in mapping. What you choose here
#' 		determines what you can choose in the region parameter. One of: county,
#' 		state, usa, world, world2, france, italy, or nz.
#' @param region The region of the world to map. From the maps package, run
#' 		\code{sort(unique(map_data("world")$region))} to see region names for the
#' 		world database layer, or e.g., \code{sort(unique(map_data("state")$region))}
#' 		for the state layer.
#' @param legend logical
#' @param interactive logical, if TRUE map will be plotted using function \code{mapview::mapview}
#' @param panel plots panels for each species for species above threshold supplied to \code{panel}, e.g., 1000
#' @param jitter If you use jitter, the amount by which to jitter
#' 		points in width, height, or both.
#' @param plot logical
#' @param gazetter logical if Gazetter should be used to add long lat for USA counties where coordinates are missing; default = TRUE

#' @return Map (using \link{ggplot2} package) of points or tiles on a world map.
#' @import ggplot2 mapview sf sp
#' @importFrom utils read.csv
#' @importFrom stats complete.cases
#'
#' @author Franz-Sebastian Krah
#'
#' @examples
#' \dontrun{
#' am.dist <- symbiota(taxon = "Amanita muscaria", db = "MyCoPortal")
#' head(recordsTable(am.dist))
#' plot_distmap(am.dist, mapdatabase = "state", interactive = TRUE)
#' plot_distmap(am.dist, mapdatabase = "state", interactive = FALSE)
#' }
#' @export

plot_distmap <- function(x,
                         mapdatabase = "world",
                         region = ".",
                         legend = FALSE,
                         panel = FALSE,
                         interactive = TRUE,
                         jitter = position_jitter(width = 0, height = 0),
                         plot = TRUE,
                         gazetter = TRUE) {


  ## Add lat long Gazeter file
  if(gazetter){
    ipt <- x@records
    fpath <- system.file("extdata", "2017_Gaz_counties_national.txt", package="rSymbiota")
    county.coord <- read.csv(fpath, sep = "\t")
    # can be found and downloaded here:
    # https://www.census.gov/geo/maps-data/data/gazetteer2017.html

    county.coord$NAME <- gsub(" County", "", county.coord$NAME)
    ipt$County <- gsub(" Co\\.", "", ipt$County)
    ipt <- data.frame(ipt, county.coord[match(ipt$County, county.coord$NAME),][,c("INTPTLAT","INTPTLONG")])

    ipt$lat[is.na(ipt$lat)] <- ipt$INTPTLAT[is.na(ipt$lat)]
    ipt$lon[is.na(ipt$lon)] <- ipt$INTPTLONG[is.na(ipt$lon)]
    x@records <- ipt
  }

  if(!interactive){
    ipt <- x@records

    tomap <- ipt[complete.cases(ipt$lat, ipt$lon), ]
    tomap <- tomap[-(which(tomap$lat <=90 || tomap$lon <=180)), ]
    world <- map_data(map=mapdatabase, region=region) # get world map data

    tomap <- tomap[(tomap$lat >= min(world$lat) & tomap$lat <= max(world$lat)) &
                     (tomap$lon >= min(world$long) & tomap$lon <= max(world$long)),]

    message(paste("Rendering map...plotting ", nrow(tomap), " points. Not all records have coordinates.", sep=""))

    if(panel){
      nr.rec <- table(tomap$spec)
      nr.rec <- names(nr.rec)[nr.rec>panel]
      tomap <- tomap[tomap$spec %in% nr.rec, ]
    }

    p <- ggplot(world, aes(world$long, world$lat)) + # make the plot
      geom_polygon(aes(group=world$group), fill="white", color="gray40", size=0.2) +
      geom_point(data=tomap, aes(tomap$lon, tomap$lat, color = tomap$species),
                 alpha = 0.4, size = 3, position = jitter) +
      labs(x="", y="") + theme_bw(base_size = 14) +
      ggtitle(paste("Distribution map for", x@query$taxon)) +
      coord_fixed(1.3)

    if(!legend)
      p <- p + theme(legend.position = "none")

    if(panel){
      p <- p + facet_wrap(~spec)
    }
    if(plot)
      plot(p)
    return(p)
  }
  if(interactive){
    x@records <- x@records[!is.na(x@records$lat), ]
    x@records <- x@records[!is.na(x@records$lon), ]
    my.sf.point <- st_as_sf(x = x@records,
                            coords = c("lon", "lat"),
                            crs = "+proj=longlat +datum=WGS84")
    mapview(my.sf.point)
  }
}