Missing Links: Density of Bigfoot Sightings in North America

October 14, 2018
By

(This article was first published on Weird Data Science, and kindly contributed to R-bloggers)

Sightings of homo sapiens cognatus, or homo sylvestris, are well-recorded, with a particular prevalence in North America. Bigfoot, otherwise known as the Sasquatch, is one of the most widely-reported cryptids in the world; it is the subject of documentaries, film, music, and countless books.

The Bigfoot Field Research Organisation has compiled a detailed database of Bigfoot sightings going back to the 1920s. Each sighting is dated, located to the nearest town or road, and contains a full description of the sighting. In many cases, sightings are accompanied by a follow-up report from the organisation itself.

As previously with UFO sightings and paranormal manifestations, our first step is to retrieve the data and parse it for analysis. Thankfully, the bfro.net dataset is relatively well-formatted; reports are broken down by region, with each report following a mainly standard format.

As before, we rely on the rvest package in R to explore and scrape the website. In this case, the key elements were to retrieve each state’s set of reports from the top level page, and retrieve the link for each report. Conveniently, these are in a standard format; the website also allows a printer-friendly mode that greatly simplifies scraping.

The scraping code is given here:

Show scraping code.
library(tidyverse)
library(progress)
library(rvest)

# Base URLs for scraping
index_url <- "https://www.bfro.net/GDB/"
base_url <- "https://www.bfro.net"
report_base_url_pattern <- "https:\\/\\/www.bfro.net\\/GDB\\/show_report.asp\\?id="

# Retrieve list of state-level county report pages
get_state_listing_urls <- function( url ) {

	read_html( url ) %>%
		html_nodes( 'a' ) %>%
		html_attr( 'href' ) %>%
		str_match( '.*state_listing.asp\\?state=.*' ) %>%		
		na.omit %>%
		unique %>%
		{ paste0( base_url, . ) }

}

# Return all URLs pointing to a county-level list of reports 
get_county_report_urls <- function( url ) {

	read_html( url ) %>%
		html_nodes( 'a' ) %>%
		html_attr( 'href' ) %>%
		str_match( '.*show_county_reports.asp\\?state=.*' ) %>%		
		na.omit %>%
		unique %>%
		{ paste0( index_url, . ) }

}

# Return all URLs pointing to a report in this page.
get_report_urls <- function( url ) {
	
	read_html( url ) %>%
		html_nodes( 'a' ) %>%
		html_attr( 'href' ) %>%
		str_match( '.*show_report.asp\\?id=[0-9]+' ) %>%
		na.omit %>%
		unique %>%
		{ paste0( index_url, . ) }

}

progressive_get_county_report_urls <- function( url, progress_bar ) {

	progress_bar$tick()$print()
	cat( paste(url, "\n") , file="work/progress.log", append=TRUE )
	chime()
	get_county_report_urls( url )

}

progressive_get_report_urls <- function( url, progress_bar ) {

	progress_bar$tick()$print()
	cat( paste(url, "\n") , file="work/progress.log", append=TRUE )
	get_report_urls( url )

}

# Pull the report listing from a page.
store_report <- function( url ) {

	report_id <-
		url %>%
		str_replace( report_base_url_pattern, "" ) %>%
		str_replace( "/", "-" ) %>%
		str_replace( ".html", "" )


	#message( paste0("Processing report ", report_id, "... " ), appendLF=FALSE )

	# Check if this report has already been stored.
	if( file.exists( paste0( "data/", report_id, ".rds" ) ) ) {
		message( paste0( "Report already retrieved: ", report_id ) )
		return()
	}

	url_pf <- paste0( url, "&PrinterFriendly=True" )

	report <- 
		tryCatch( 
					{
						# Fetch and parse HTML of current page.
						read_html( url_pf ) %>%
							as.character
					},
					error=function( cond ) {
						message( paste( "URL caused an error:", url ))
						message( "Error message:")
						message( cond )
						return( NULL )
					},
					warning=function( cond ) {
						message( paste( "URL caused a warning:", url ))
						message( "Warning message:")
						message( cond )
						return( NULL )
					},
					finally={
					}
					)

	# Write report result to disk
	saveRDS( report, paste0("data/", report_id, ".rds") )

}

# Progressive version of store_report
progressive_store_report <- function( url, progress_bar ) {

	progress_bar$tick()$print()
	cat( paste(url, "\n") , file="work/progress.log", append=TRUE )
	store_report( url )

}

# The key to this is that links of the form 
#  
# link to each state's data. At the top level, then, we can just get all links matching that form.

# At each sub-page, the links are:
# 
# We can just pull all of those for each state.
# NOTE: There are also non-US reports linked like this directly from the main GDB page.

# Finally, it seems that adding "&PrinterFriendly=True" will strip a lot of
# unnecessary formatting.

# From the top-level page, get every URL that matches 'show_county_reports.asp?state=', make
# the list unique, then download the pages.

# Get all non-US state indices
message("Getting state report lists...", appendLF=FALSE )
bfro_state_indices <- get_state_listing_urls( index_url )
message("done." )
saveRDS( bfro_state_indices, "work/bfro_state_indices.rds" )

# Get all US county-level indices
message("Getting county-level report urls... " )
pb <- progress_estimated( length( bfro_state_indices ) )
bfro_county_urls <- 
	bfro_state_indices %>%
	map( progressive_get_county_report_urls, pb ) %>%
	unlist %>%
	unique
saveRDS( bfro_county_urls, "work/bfro_county_urls.rds" )

# Get URLs of US reports from county pages.
# (Final subset line handles cases where pages are empty of reports, but
# contain other links such as media articles.)
pb <- progress_estimated( length( bfro_county_urls ) )
bfro_report_urls <- 
	bfro_county_urls %>%
	map( progressive_get_report_urls, pb ) %>%
	unlist %>%
	unique %>%
	str_subset( "GDB\\/." )
saveRDS( bfro_report_urls, "work/bfro_report_urls.rds" )

# Store all US reports
pb <- progress_estimated( length( bfro_report_urls ) )
bfro_report_urls %>%
	map( progressive_store_report, pb )

# Get all non-US indices
message("Getting non-US top-level report lists...", appendLF=FALSE )
bfro_nonus_indices <- 
	get_county_report_urls( index_url ) %>%
	str_replace( "GDB\\/\\/", "\\/" )
message("done." )
saveRDS( bfro_nonus_indices, "work/bfro_nonus_indices.rds" )

# Get URLs of US reports from county pages
message("Getting non-US report lists...", appendLF=FALSE )
pb <- progress_estimated( length( bfro_nonus_indices ) )
bfro_nonus_report_urls <- 
	bfro_nonus_indices %>%
	map( progressive_get_report_urls, pb ) %>%
	unlist %>%
	unique  %>%
	str_subset( "GDB\\/." )
saveRDS( bfro_nonus_report_urls, "work/bfro_nonus_report_urls.rds" )

# Store all US reports
pb <- progress_estimated( length( bfro_nonus_report_urls ) )
bfro_nonus_report_urls %>%
	map( progressive_store_report, pb )

With each page retrieved, we step through and parse each report. Again, each page is fairly well-formatted, and uses a standard set of tags for date, location, and similar. The report parsing code is given here:

Show parsing code.
# NOTES

# Process entire list to check for standardised headers. These are in capitals, and located in  tags.

# Report starts with:
# 
#  (Look up what these mean.)
# Some following  types contain general summary information, eg:
# " Submitted  by  witness   on Thursday, November 1, 2007."
# Following fields are in span tags separated by paragraph tags, eg:
# 

YEAR: 2007

# STATE and COUNTY fields are typically links; we should pull their text. (I can't see a good reason to parse links to anything other than the link text for our purposes.) library(tidyverse) library(magrittr) library(progress) library(rvest) library(lubridate) # Exploratory functions # List all capitalised fields seen in any file list_all_fields <- function() { filenames <- list.files( "data", pattern="*.rds", full.names=TRUE) # In each file, select the classes, match the uppercase field names, and extract text. all_fields <- filenames %>% map( list_report_fields ) %>% unlist %>% unique saveRDS( all_fields, "work/fields.rds" ) return( all_fields ) } fields_to_colnames <- function( fields ) { # In total there are 18 fields reported in the data, with the final one # being an apparently miscoded date and place from a report. # Format the text appropriately for dataframe column names fields %>% head( 17 ) %>% str_remove( ":" ) %>% str_replace_all( " ", "_" ) %>% str_to_lower() } list_report_fields <- function( report ) { to_process <- readRDS( report ) bfro_fields <- NULL if( !is.null( to_process ) ) { bfro_fields <- to_process %>% read_html %>% html_nodes( "span" ) %>% html_nodes(xpath = '//*[@class="field"]') %>% html_text %>% str_subset( "[A-Z]+:" ) } return( bfro_fields ) } # Extract node date parse_report_data <- function( report_data ) { # report_data should be an xml_nodeset # Metadata metadata_list <- report_data %>% html_text %>% str_subset( "^[A-Z ]+: " ) %>% str_split_fixed( ": ", 2 ) %>% as.tibble %>% spread( key=V1, value=V2 ) %>% set_colnames( fields_to_colnames( colnames(.) ) ) # Extra details extra_text <- report_data %>% html_text %>% str_remove( "^[A-Z ]+:.*" ) %>% str_flatten( " " ) %>% str_trim # Add extra details as a column metadata_list$extra <- extra_text # Note whether date is rough or accurately reported metadata_list$rough_date <- FALSE # "Fix" missing date or month columns if( "date" %in% colnames( metadata_list ) == FALSE ) { metadata_list$date <- "1" metadata_list$rough_date <- TRUE } if( "month" %in% colnames( metadata_list ) == FALSE ) { metadata_list$month <- "January" metadata_list$rough_date <- TRUE } # Combine date columns (year, month, date) to a true date. metadata_list <- metadata_list %>% mutate( year = str_replace( year, ".*([0-9]{4}).*", "\\1" ) ) %>% mutate( date = paste0( year, "-", month, "-", date ), year=NULL, month=NULL ) %>% mutate( date = as.POSIXct( date, format="%Y-%B-%d" ) ) } # Read a file and pass to `post_get_all_thread` process_file <- function( filename ) { # Read stored file to_process <- readRDS( filename ) # Don't process null files if( is.null( to_process ) ) return( NULL ) if( length( to_process ) == 0 ) return( NULL ) # Produce an xml_nodeset for parsing xml_thread <- to_process %>% read_html %>% html_nodes( "p" ) parse_report_data( xml_thread ) } # Progressive version of process_file progressive_process_file <- function( filename, progress_bar ) { progress_bar$tick()$print() cat( paste(filename, "\n") , file="progress.log", append=TRUE ) report <- tryCatch( { process_file( filename ) }, error=function( cond ) { message( paste( "File caused an error:", filename )) message( "Error message:") message( cond ) return( NULL ) }, warning=function( cond ) { message( paste( "URL caused a warning:", url )) message( "Warning message:") message( cond ) return( NULL ) }, finally={ } ) return( report ) } ## Processing starts here # Read all .rds data files and process into a thread filenames <- list.files("data", pattern="*.rds", full.names=TRUE) pb <- progress_estimated( length( filenames ) ) # Begin logging cat( "Processing... ", file="progress.log", append=FALSE ) bfro_tbl <- filenames %>% map( progressive_process_file, pb ) %>% bind_rows saveRDS( bfro_tbl, "data/bfro_processed.rds" )

With each report parsed into a form suitable for analysis, the final step in scraping the site is to geolocate the reports. As in previous posts, we rely on Google’s geolocation API. For each report, we extract an appropriate address and parse it into a set of latitude and longitude coordinates. For the purposes of this initial scrape we restrict ourselves to North America, which compromises a large majority of reports on `bfro.net`. Geolocation code is included below. (Note that a Google Geolocation API key is required for this code to run.)

Show geolocation code.
library(googleway)
library(progress)
library(tidyverse)
library(magrittr)

# weird.data.science Google API key (Geocoding API enabled)
key <- 

# Load bfro tibble
bfro_data <- readRDS( "data/bfro_processed.rds" )

# Geocode entries

# BFRO entries contain some, all, or none of:
# country, province, state, county, nearest_town

# Country is only used for Canada, in which case a province is given.
# Country and province are NA for the US.

# Best plan, then is to create a string of (nearest_town, province, country) for Canadian results, and (nearest_town, county, state, "US") for US results.

# Bound the request to be only in North America.
# (Used: )
# (SW->NE latitude and longitude.)
bounding_box <- list( c( -170.3, 	24.4),
							 c(  -52.3,		83.3) )

form_location_string <- function( country, province, state, county, nearest_town ) {

	location_string <- NA

	# US case, then Canadian
	if( is.na( country ) )
		location_string <- paste0( nearest_town, ", ", county, ", ", state, ", US" )
	else 
		location_string <- paste0( nearest_town, ", ", province, ", ", ", Canada" )

	# Strip double commas caused by NA values and remove bracketed clauses.
	location_string %<>%
		str_remove_all( "\\([^\\)]*\\)" ) %>%
		str_replace_all( ", , ", ", " ) %>%
		str_replace_all( " ,", "," ) %>%
		str_remove_all( "NA, " )

}

# Create a safe version of google_geocode
safe_geocode <- safely( google_geocode )

# Wrap google_geocode with a progress bar call.
# Also, optionally remove any bracketed substrings entirely (strip_brackets)
progressive_geocode <- function( location_string, progress_bar ) {

	# Write status to log file.
	progress_bar$tick()$print()

	cat( paste0( location_string, "\n" ), file="progress.log", append=TRUE )

	result <- 
		safe_geocode( 	
						 address = location_string,
						 key = key, 
						 bounds = bounding_box, 
						 simplify = TRUE ) 

	# Note that this can be NULL
	return( result$result )

}

# Logging and output
## Clear the screen first
cat(c("\033[2J","\033[0;0H"))
cat("Geolocating entries...\n")
cat( "Geolocating entries...\n", file="progress.log", append=FALSE )

# Create the location string
bfro_data %<>%
	mutate( location = form_location_string( country, province, state, county, nearest_town ) )

# Geolocate each location.
pb <- progress_estimated(nrow(bfro_data))
bfro_data$geolocation <- 
	map( bfro_data$location, progressive_geocode, progress_bar = pb )

cat("\nComplete.\n")

## With all values scraped, save geolocated data.
saveRDS( bfro_data, file="work/bfro_data_geolocated.rds")

With geolocated data in hand, we can now venture into the wilds. In which areas of North America are Sasquatch most commonly seen to roam? The plot below shows the overall density of Bigfoot sightings, with individual reports marked.

Density plot of Bigfoot sightings in North America

Density of Bigfoot sightings in North America. (PDF Version)

There are particular clusters on the Great Lakes, particularly in Southern Ontario; as well as in the Pacific Northwest. Smaller notable clusters exist in Florida, centered around Orlando. As with most report-based datasets, sightings are skewed towards areas of high population density.

The obvious first question to ask of such data is which, if any, environmental features correlate with these sightings. Other analyses of Bigfoot sightings, such as the seminal work of Lozier et al.% filter( !is.na( lat ) & !is.na( lng ) ) %>% filter( not( lat == modal( lat ) & lng == modal( lng ) ) ) # Convert the bfro dataframe to a spatial dataframe that contains # explicit longitude and latitude projected appropriately for plotting. coordinates( bfro_tbl ) <- ~lng+lat proj4string( bfro_tbl ) <- CRS("+init=epsg:4326 +lon_wrap=170") bfro_tbl_spatial <- spTransform(bfro_tbl,CRS(proj4string(bfro_tbl))) # Restrict bfro_tbl to those points in the polygons defined by world_subset bfro_tbl_rows <- bfro_tbl_spatial %>% over( world_subset ) %>% is.na() %>% not() %>% rowSums() %>% `!=`(0) %>% which bfro_tbl <- as.tibble( bfro_tbl_spatial[ bfro_tbl_rows, ] ) # Create window for spatial analysis bfro_owin <- as.owin.SpatialPolygons(world_subset) # Function to plot density of a specific manifestation type. # plot_resolution is for the density raster, and is mainly used for quick prototyping of the output. density_plot <- function( density_res ) { cat( paste0( "Plotting density: ... " ) ) bfro_ppp <- ppp( x=coordinates(bfro_tbl_spatial)[,1], y=coordinates(bfro_tbl_spatial)[,2], window = bfro_owin ) # This discards 'illegal' points outside of the window bfro_ppp <- as.ppp(bfro_ppp) bfro_density <- density( bfro_ppp, diggle=T, sigma=2, dimyx = c( density_res, density_res ) ) # Make density image object usable by ggplot as a raster bfro_density_raster <- raster( bfro_density ) bfro_raster_tbl <- as.tibble( rasterToPoints( bfro_density_raster ) ) # Show the map gp <- ggplot() # Add density of sightings as raster. gp <- gp + geom_raster( data = bfro_raster_tbl, aes( x=x, y=y, fill = layer ), alpha=0.8, show.legend=FALSE ) + scale_fill_viridis( option = "cividis", direction=1 ) # Add sightings as points. gp <- gp + geom_point( data = as.tibble( bfro_ppp ), aes( x=x, y=y ), colour="yellow", size=0.1, alpha=0.5, show.legend=FALSE ) gp <- gp + geom_map( data = world_df, aes( map_id=id ), colour="#3c3f4a", fill = "transparent", size = 0.4, map = world_df ) # Theming gp <- gp + theme_map() gp <- gp + theme( plot.background = element_rect(fill = "transparent", colour = "transparent"), panel.border = element_blank(), ) gp <- gp + guides( fill = guide_colourbar( title.position="top", direction="horizontal", barwidth=6, barheight=0.4 ) ) cat("done.\n" ) return(gp) } # Add bioclimatic points to an existing plot add_bioclimatic_variables <- function( gp, bioclim = seq( 1, 16 ) ) { # This gets climatic variables for global range. # (res=0.5 requires latitude and longitude to define a tile to download, # but 2.5 minutes of a degree is ~5 miles, so probably good enough.) #wc <- getData( 'worldclim', res=2.5, var='bio' ) # Value Label # 0 Water # 1 Evergreen Needleleaf forest # 2 Evergreen Broadleaf forest # 3 Deciduous Needleleaf forest # 4 Deciduous Broadleaf forest # 5 Mixed forest # 6 Closed shrublands # 7 Open shrublands # 8 Woody savannas # 9 Savannas # 10 Grasslands # 11 Permanent wetlands # 12 Croplands # 13 Urban and built-up # 14 Cropland/Natural vegetation mosaic # 15 Snow and ice # 16 Barren or sparsely vegetated # 254 Unclassified # 255 Fill Value glcf <- raster( "~/opt/datasets/gis/landcover/glcf/LC_5min_global_2012.tif" ) # We can't just reproject a raster to wrap differently. # This splits and merges the two edges. # (This is fragile, unfortunately, but works here.) # (Reprojecting is best done as the last step, apparently.) glcf_west <- crop(glcf, extent(-180, 0, 18, 84)) glcf_east <- crop(glcf, extent(0, 180, 18, 84)) extent(glcf_west) <- c(180, 360, 18, 84) glcf <- merge(glcf_west, glcf_east ) glcf <- crop( glcf, extent( world_subset ) ) glcf <- projectRaster( glcf, crs = CRS("+init=epsg:4326") ) # Make density image object usable by ggplot as a raster glcf_raster_tbl <- as.tibble( rasterToPoints( glcf ) ) colnames( glcf_raster_tbl ) <- c ("x", "y", "layer" ) # Add bioclimatic variables gp <- gp + geom_raster( data = glcf_raster_tbl[ which( glcf_raster_tbl$layer %in% bioclim ), ], alpha=0.8, aes( x=x, y=y ), fill="#7cfc00", show.legend=FALSE ) gp } # Cowplot for full-panel plotting. (Parchment background). theme_set(theme_cowplot(font_size=4, font_family = "mapfont" ) ) # Only calculate the density plot if it isn't already stored on disk. if( not( file.exists( "work/density_plot.rds" ) ) ) { gp <- density_plot( 1024 ) saveRDS( gp, "work/density_plot.rds" ) } else { gp <- readRDS( "work/density_plot.rds" ) } # Create a plot including the bioclimatic variables. bioclim_gp <- add_bioclimatic_variables( gp, seq( 1, 5 ) ) # Construct full plot, with title and backdrop. title <- ggdraw() + draw_label("Bigfoot Sightings in North America", fontfamily="mapfont", colour = "#3c3f4a", size=20, hjust=0, vjust=1, x=0.02, y=0.88) + draw_label("http://www.weirddatascience.net | @WeirdDataSci", fontfamily="mapfont", colour = "#3c3f4a", size=12, hjust=0, vjust=1, x=0.02, y=0.40) data_label <- ggdraw() + draw_label("Data: http://www.bfro.net", fontfamily="mapfont", colour = "#3c3f4a", size=12, hjust=1, x=0.98 ) tgp <- plot_grid(title, gp, data_label, ncol=1, rel_heights=c(0.1, 1, 0.1)) parchment_plot <- ggdraw() + draw_image("img/parchment.jpg", scale=1.4 ) + draw_plot(tgp) save_plot("output/bigfoot-density.pdf", parchment_plot, base_width = 16, base_height = 9, base_aspect_ratio = 1.78 ) # Add bioclimatic variables to plot. title <- ggdraw() + draw_label("Bigfoot Sightings in North America Showing Forest Cover", fontfamily="mapfont", colour = "#3c3f4a", size=20, hjust=0, vjust=1, x=0.02, y=0.88) + draw_label("http://www.weirddatascience.net | @WeirdDataSci", fontfamily="mapfont", colour = "#3c3f4a", size=12, hjust=0, vjust=1, x=0.02, y=0.40) data_label <- ggdraw() + draw_label("Data: http://www.bfro.net", fontfamily="mapfont", colour = "#3c3f4a", size=12, hjust=1, x=0.98 ) tgp <- plot_grid(title, bioclim_gp, data_label, ncol=1, rel_heights=c(0.1, 1, 0.1)) parchment_plot <- ggdraw() + draw_image("img/parchment.jpg", scale=1.4 ) + draw_plot(tgp) save_plot("output/bigfoot-density-forest.pdf", parchment_plot, base_width = 16, base_height = 9, base_aspect_ratio = 1.78 )

From this initial plot we can see that, whilst tree cover is certainly not a bad predictor of Bigfoot sightings, it is far from a definitive correlation. The largest cluster, around the US-Canada border near Toronto, is principally lakes; whilst the secondary cluster in Florida is neither significantly forested or even close to the Everglades, which might have been expected. From the other perspective, there are significant forested areas for which sightings are reasonably rare.

The mystery of Bigfoot’s natural habitat and preferences is, therefore, very much unanswered from our initial analysis. With a broad range of variables still to explore — climate, altitude, food sources — future posts will attempt to determine what conditions lend themselves to strange survivals of pre-human primate activity. Perhaps changing conditions have pushed our far-distant cousins to previously unsuspected regions.

Until then, we keep following these trails into data unknown.

References

To leave a comment for the author, please follow the link and comment on their blog: Weird Data Science.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers

Most visited articles of the week

  1. How to write the first for loop in R
  2. Learning R: The Ultimate Introduction (incl. Machine Learning!)
  3. Modern Data Science with R: A review
  4. R Studio Shortcuts and Tips – part 2
  5. 5 Ways to Subset a Data Frame in R
  6. How to interactively examine any R code - 4 ways to not just read the code, but delve into it step-by-step
  7. R – Sorting a data frame by the contents of a column
  8. Create Animation in R : Learn by Examples
  9. Using apply, sapply, lapply in R

Sponsors

RSS Jobs for R users

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)