PleistoDistR: Calculating Island Geomorphologies Over Pleistocene Time

Documented in pleistodist_centroid pleistodist_euclidean pleistodist_leastcost pleistodist_leastshore pleistodist_meanshore pleistodist_netmig pleistodist_visibility

#' Calculate Euclidean distance between points
#'
#' This function calculates the straight-line distance ("as the crow flies) between user-specified points.
#' Because these points are invariant across space and time and independent of sea level change, this calculation
#' is only performed once instead of repeatedly across different sea level or time intervals.
#'
#' @param points A user-generated multi-point shapefile (.shp) containing at least two points. If the shapefile attribute table contains a column
#' labelled 'Name', the output distance matrix will use the identifiers in this column to label each pairwise comparison. Otherwise, the
#' output distance matrix will use the attribute FID values instead.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param outdir Output directory for the Euclidean distance matrix (island_euclideandist.csv).
#' If the specified output directory doesn't already exist, PleistoDist will create the output directory.
#' @return This function outputs a pairwise distance matrix of inter-point Euclidean distances in long format.
#' @examples
#' #create temp directory
#' path <- file.path(tempdir())
#' #create points shapefile
#' points <- sf::st_multipoint(rbind(c(179.41256,-17.79426), c(179.27600,-17.97850)))
#' #convert points to feature geometry
#' points <- sf::st_sfc(points)
#' points <- sf::st_cast(points, "POINT")
#' #set default projection (WGS84)
#' sf::st_crs(points) <- 4326
#' #save shapefile
#' sf::write_sf(points,paste0(path,"/points.shp"))
#' #Project input points using EPSG:3141 and calculate Euclidean distances
#' pleistodist_euclidean(points=paste0(path,"/points.shp"),epsg=3141,outdir=path)
#' @export
pleistodist_euclidean <- function(points,epsg,outdir) {

  if(base::dir.exists(outdir)==FALSE) {
    base::dir.create(outdir)
  }

  points = sf::st_read(points,fid_column_name="FID")

  #check projection on input points, set to WGS84 (EPSG:4326) if undefined
  if (is.na(sf::st_crs(points)$input)) {
    sf::st_crs(points) = 4326
  }

  #reproject input points to target projection
  points_transformed <- sf::st_transform(points,epsg)

  #check to see if the points contains a column of unique names, otherwise use FID numbers as identifiers
  if ("Name" %in% colnames(points_transformed) && anyDuplicated(points_transformed$Name) == 0) {
    p1_names <- points_transformed$Name
    p2_names <- points_transformed$Name

  } else {
    p1_names <- points_transformed$FID
    p2_names <- points_transformed$FID
    message("No 'Name' column with unique identifiers detected in points file, defaulting to FID values instead")
  }

  #create dataframe for storing Euclidean distance matrix
  euclideandist <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names,include.equals=T)[,2],
    Island2 = expand.grid(p1_names,p2_names,include.equals=T)[,1]
  )

  euclidean <- sf::st_distance(points_transformed)
  blankmatrix <- expand.grid(1:nrow(points_transformed),1:nrow(points_transformed))
  blankmatrix$interval0 <- tidyr::gather(as.data.frame(euclidean))$value
  #blankmatrix <- blankmatrix %>%
  #ilter(Var1>=Var2)
  euclideandist$interval0 <- as.numeric(blankmatrix$interval0)
  utils::write.csv(euclideandist,base::paste0(outdir,"/island_euclideandist.csv"))
}

#' Calculate inter-island centroid-to-centroid distance over time
#'
#' This function calculates the centroid-to-centroid distance between island landmasses over Pleistocene time. For each pairwise combination
#' of source points provided by the user, this function selects the landmasses that correspond to each point, and calculates the distance between
#' the centroids of the selected landmasses for each time/sea level interval. If both points lie on the same landmass, a distance value of 0 is returned. If one or both islands
#' is/are underwater, then a value of NA is returned. After calculating the pairwise island distances for each interval, this function calculates
#' the time-weighted average centroid-to-centroid distance between each island pair.
#'
#' @param points A user-generated multi-point shapefile (.shp) containing at least two points. If the shapefile attribute table contains a column
#' labelled 'Name', the output distance matrix will use the identifiers in this column to label each pairwise comparison. Otherwise, the
#' output distance matrix will use the attribute FID values instead.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param intervalfile This is the master control file generated using either the getintervals_time() or getintervals_sealvl() functions, or custom generated by the user,
#' that defines the number of intervals, the sea level at each interval, and the duration of each interval.
#' @param mapdir The directory containing map outputs from the makemaps() function
#' (i.e. the directory containing the raster_flat, raster_topo, and shapefile folders).
#' @param outdir Output directory for the inter-island centroid-to-centroid distance matrix file.
#' If the specified output directory doesn't already exist, PleistoDist will create the output directory.
#' @return This function outputs a pairwise distance matrix of pairwise inter-island centroid-to-centroid distances in long format, with
#' one column per interval. This function also generates a file of pairwise relative island widths, which is used for calculating the net inter-island migration ratio.
#' @examples
#' \donttest{
#' #create temp directory
#' path <- file.path(tempdir())
#' #create points shapefile
#' points <- sf::st_multipoint(rbind(c(179.41256,-17.79426), c(179.27600,-17.97850)))
#' #convert points to feature geometry
#' points <- sf::st_sfc(points)
#' points <- sf::st_cast(points, "POINT")
#' #set default projection (WGS84)
#' sf::st_crs(points) <- 4326
#' #save shapefile
#' sf::write_sf(points,paste0(path,"/points.shp"))
#' #load bathymetry file
#' fiji <- system.file("extdata","FJ.asc",package="PleistoDist")
#' #generate interval file for 2 intervals and 20 kya cutoff time, binning by time
#' getintervals_time(time=20,intervals=2,outdir=path)
#' #generate maps based on interval file, projecting map using EPSG:3141
#' makemaps(inputraster=fiji,epsg=3141,intervalfile=paste0(path,"/intervals.csv"),outdir=path)
#' #calculate inter-island centroid-to-centroid distances, projecting points using EPSG:3141
#' pleistodist_centroid(points=paste0(path,"/points.shp"),epsg=3141,
#'     intervalfile=paste0(path,"/intervals.csv"), mapdir=path,outdir=path)
#' }
#' @export
pleistodist_centroid <- function(points,epsg,intervalfile,mapdir,outdir) {
  #check for existence of output directory, create output directory if it doesn't already exist.
  if(base::dir.exists(outdir)==FALSE) {
    base::dir.create(outdir)
  }
  intervalfile = utils::read.csv(intervalfile)
  points = sf::st_read(points,fid_column_name="FID")

  #check projection on input points, set to WGS84 (EPSG:4326) if undefined
  if (is.na(sf::st_crs(points)$input)) {
    sf::st_crs(points) = 4326
  }

  #reproject input points to target projection
  points_transformed <- sf::st_transform(points,epsg)
  numintervals = max(intervalfile$Interval)

  #check to see if the points contains a column of unique names, otherwise use FID numbers as identifiers
  if ("Name" %in% colnames(points_transformed) && anyDuplicated(points_transformed$Name) == 0) {
    p1_names <- points_transformed$Name
    p2_names <- points_transformed$Name

  } else {
    p1_names <- points_transformed$FID
    p2_names <- points_transformed$FID
    message("No 'Name' column with unique identifiers detected in points file, defaulting to FID values instead")
  }

  #create dataframe for storing centroid-to-centroid distance matrix
  centroiddist <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names,include.equals=T)[,2],
    Island2 = expand.grid(p1_names,p2_names,include.equals=T)[,1]
  )

  #create dataframe for storing relative island width matrix
  relativeislandwidth <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names)[,2],
    Island2 = expand.grid(p1_names,p2_names)[,1]
  )

  #for each interval...
  for (i in 0:numintervals) {
    #load interval shapefile
    invector <- sf::st_read(base::paste0(mapdir,"/shapefile/interval",i,".shp")) #load interval shapefile
    intvlname <- base::paste0("interval",i)

    #create empty lists for storing distance values at each interval
    centroidd <- c()
    iwidth <- c()

    message(base::paste0("Analysing island distances for interval ",i,"..."))
    #start of first order for-loop to select island 1
    for (f in 1:nrow(points_transformed)) {

      #start second order for-loop to select island 2
      for (t in 1:nrow(points_transformed)) {

        #extract coordinates for island 1
        p1 <- points_transformed$geometry[f]
        #extract coordinates for island 2
        p2 <- points_transformed$geometry[t]

        message(base::paste0("Calculating distance between points ",p1_names[f]," and ",p2_names[t],"... "))

        if (f == t) { #write NA for self-comparisons
          centroidd <- c(centroidd,NA)
          iwidth <- c(iwidth,NA)
        } else if ((nrow(invector[p1, ,op=sf::st_contains]) || nrow(invector[p2, ,op=sf::st_contains])) == 0) { #check to see if island1 or island2 is submerged under water
          centroidd <- c(centroidd,NA)
          iwidth <- c(iwidth,NA)
          message("One or both islands are underwater during this interval, writing a value of NA")
        } else {
          combinedpoints <- sf::st_sfc(rbind(p1,p2),crs=epsg)
          islandpair <- invector[combinedpoints, ,op=sf::st_contains]
          if (nrow(islandpair) == 1) { #check to see if both points are on the same island, if TRUE write distance of 0 for island-based comparisons
            centroidd <- c(centroidd,0)
            iwidth <- c(iwidth,NA)
            message("Both points are on the same island, writing a value of 0 for centroid-to-centroid distance")
          } else {
            island1 <- invector[p1, ,op=sf::st_contains]
            island2 <- invector[p2, ,op=sf::st_contains]
            island1_centroid <- sf::st_centroid(island1)
            island2_centroid <- sf::st_centroid(island2)
            centroidline <- sf::st_linestring(rbind(sf::st_coordinates(island1_centroid),sf::st_coordinates(island2_centroid)))

            m = -1/(stats::line(centroidline)$coefficients[2])
            c = sf::st_coordinates(island1_centroid)[,2]-(m*sf::st_coordinates(island1_centroid)[,1])

            #calculate 2 points on widthline at edges of bounding box
            x1 = sf::st_bbox(island1)$xmax
            y1 = m*x1+c

            x2 = sf::st_bbox(island1)$xmin
            y2 = m*x2+c

            coord1 <- sf::st_point(c(x1,y1))
            coord2 <- sf::st_point(c(x2,y2))
            widthline <- sf::st_sfc(sf::st_linestring(c(coord1,coord2)),crs=epsg)
            #because complex island shapes can slice widthlines into multiple fragments, convert widthline into coordinates,
            #generate a new line, and calculate length from this new line
            islandwidth<- sf::st_length(sf::st_linestring(sf::st_coordinates(sf::st_intersection(widthline,island1))))
            iwidth = c(iwidth,islandwidth)

            cdistance <- as.numeric(sf::st_length(centroidline))
            centroidd <- c(centroidd,cdistance)
            message(base::paste0("Distance between island centroids is ",cdistance,"m"))
          }
        }
      }
    }
    centroiddist[intvlname] <- centroidd
    relativeislandwidth[intvlname] <- iwidth
  }
  centroiddist$mean <- apply(centroiddist[3:ncol(centroiddist)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  relativeislandwidth$mean <- apply(relativeislandwidth[3:ncol(relativeislandwidth)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  utils::write.csv(centroiddist,paste0(outdir,"/island_centroiddist.csv"))
  utils::write.csv(relativeislandwidth,paste0(outdir,"/island_relativewidth.csv"))
}

#' Calculate least inter-island shore-to-shore distance over time
#'
#' This function calculates the least shore-to-shore distance between island landmasses over Pleistocene time. For each pairwise combination
#' of source points provided by the user, this function selects the landmasses that correspond to each point, and calculates the minimum
#' shore-to-shore distance between the selected landmasses for each time/sea level interval. If both points lie on the same landmass, a distance value of 0 is returned. If one or both islands
#' is/are underwater, then a value of NA is returned. After calculating the pairwise island distances for each interval, this function calculates
#' the time-weighted average least shore-to-shore distance between each island pair.
#'
#' @param points A user-generated multi-point shapefile (.shp) containing at least two points. If the shapefile attribute table contains a column
#' labelled 'Name', the output distance matrix will use the identifiers in this column to label each pairwise comparison. Otherwise, the
#' output distance matrix will use the attribute FID values instead.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param intervalfile This is the master control file generated using either the getintervals_time() or getintervals_sealvl() functions, or custom generated by the user,
#' that defines the number of intervals, the sea level at each interval, and the duration of each interval.
#' @param mapdir The directory containing map outputs from the makemaps() function
#' (i.e. the directory containing the raster_flat, raster_topo, and shapefile folders).
#' @param outdir Output directory for the inter-island least shore-to-shore distance matrix file.
#' If the specified output directory doesn't already exist, PleistoDist will create the output directory.
#' @return This function outputs a pairwise distance matrix of pairwise inter-island least shore-to-shore distances in long format, with
#' one column per interval.
#' @examples
#' \donttest{
#' #create temp directory
#' path <- file.path(tempdir())
#' #create points shapefile
#' points <- sf::st_multipoint(rbind(c(179.41256,-17.79426), c(179.27600,-17.97850)))
#' #convert points to feature geometry
#' points <- sf::st_sfc(points)
#' points <- sf::st_cast(points, "POINT")
#' #set default projection (WGS84)
#' sf::st_crs(points) <- 4326
#' #save shapefile
#' sf::write_sf(points,paste0(path,"/points.shp"))
#' #load bathymetry file
#' fiji <- system.file("extdata","FJ.asc",package="PleistoDist")
#' #generate interval file for 2 intervals and 20 kya cutoff time, binning by time
#' getintervals_time(time=20,intervals=2,outdir=path)
#' #generate maps based on interval file, projecting map using EPSG:3141
#' makemaps(inputraster=fiji,epsg=3141,intervalfile=paste0(path,"/intervals.csv"),outdir=path)
#' #calculate inter-island least shore-to-shore distances, projecting points using EPSG:3141
#' pleistodist_leastshore(points=paste0(path,"/points.shp"),epsg=3141,
#'     intervalfile=paste0(path,"/intervals.csv"),mapdir=path,outdir=path)
#' }
#' @export
pleistodist_leastshore <- function(points,epsg,intervalfile,mapdir,outdir) {
  #check for existence of output directory, create output directory if it doesn't already exist.
  if(base::dir.exists(outdir)==FALSE) {
    base::dir.create(outdir)
  }
  intervalfile = utils::read.csv(intervalfile)
  points = sf::st_read(points,fid_column_name="FID")

  #check projection on input points, set to WGS84 (EPSG:4326) if undefined
  if (is.na(sf::st_crs(points)$input)) {
    sf::st_crs(points) = 4326
  }

  #reproject input points to target projection
  points_transformed <- sf::st_transform(points,epsg)
  numintervals = max(intervalfile$Interval)

  #check to see if the points contains a column of unique names, otherwise use FID numbers as identifiers
  if ("Name" %in% colnames(points_transformed) && anyDuplicated(points_transformed$Name) == 0) {
    p1_names <- points_transformed$Name
    p2_names <- points_transformed$Name

  } else {
    p1_names <- points_transformed$FID
    p2_names <- points_transformed$FID
    message("No 'Name' column with unique identifiers detected in points file, defaulting to FID values instead")
  }

  #create dataframe for storing shore-to-shore distance matrix
  shoredist <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names,include.equals=T)[,2],
    Island2 = expand.grid(p1_names,p2_names,include.equals=T)[,1]
  )

  #create dataframe for storing relative island width matrix
  relativeislandwidth <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names)[,2],
    Island2 = expand.grid(p1_names,p2_names)[,1]
  )

  #for each interval...
  for (i in 0:numintervals) {
    #load interval shapefile
    invector <- sf::st_read(base::paste0(mapdir,"/shapefile/interval",i,".shp")) #load interval shapefile
    intvlname <- base::paste0("interval",i)
    #create empty lists for storing distance values at each interval
    shored <- c()
    iwidth <- c()
    message(base::paste0("Analysing island distances for interval ",i,"..."))
    #start of first order for-loop to select island 1
    for (f in 1:nrow(points_transformed)) {
      #start second order for-loop to select island 2
      for (t in 1:nrow(points_transformed)) {
        #extract coordinates for island 1
        p1 <- points_transformed$geometry[f]
        #extract coordinates for island 2
        p2 <- points_transformed$geometry[t]
        message(base::paste0("Calculating distance between points ",p1_names[f]," and ",p2_names[t],"... "))
        if (f == t) { #write NA for self-comparisons
          shored <- c(shored,NA)
          iwidth <- c(iwidth,NA)
        } else if ((nrow(invector[p1, ,op=sf::st_contains]) || nrow(invector[p2, ,op=sf::st_contains])) == 0) { #check to see if island1 or island2 is submerged under water
          shored <- c(shored,NA)
          iwidth <- c(iwidth,NA)
          message("One or both islands are underwater during this interval, writing a value of NA")
        } else {
          combinedpoints <- sf::st_sfc(rbind(p1,p2),crs=epsg)
          islandpair <- invector[combinedpoints, ,op=sf::st_contains]
          if (nrow(islandpair) == 1) { #check to see if both points are on the same island, if TRUE write distance of 0 for island-based comparisons
            shored <- c(shored,0)
            iwidth <- c(iwidth,0)
            message("Both points are on the same island, writing a value of 0 for least shore-to-shore distance")
          } else {
            island1 <- invector[p1, ,op=sf::st_contains]
            island2 <- invector[p2, ,op=sf::st_contains]
            island1_centroid <- sf::st_centroid(island1)
            island2_centroid <- sf::st_centroid(island2)
            centroidline <- sf::st_linestring(rbind(sf::st_coordinates(island1_centroid),sf::st_coordinates(island2_centroid)))

            m = -1/(stats::line(centroidline)$coefficients[2])
            c = sf::st_coordinates(island1_centroid)[,2]-(m*sf::st_coordinates(island1_centroid)[,1])

            #calculate 2 points on widthline at edges of bounding box
            x1 = sf::st_bbox(island1)$xmax
            y1 = m*x1+c

            x2 = sf::st_bbox(island1)$xmin
            y2 = m*x2+c

            coord1 <- sf::st_point(c(x1,y1))
            coord2 <- sf::st_point(c(x2,y2))
            widthline <- sf::st_sfc(sf::st_linestring(c(coord1,coord2)),crs=epsg)
            #because complex island shapes can slice widthlines into multiple fragments, convert widthline into coordinates,
            #generate a new line, and calculate length from this new line
            islandwidth<- sf::st_length(sf::st_linestring(sf::st_coordinates(sf::st_intersection(widthline,island1))))
            iwidth = c(iwidth,islandwidth)

            sdistance <- as.numeric(sf::st_distance(island1,island2))
            shored <- c(shored,sdistance)
            message(base::paste0("Least shore-to-shore distance between islands is ",sdistance,"m"))
          }
        }
      }
    }
    shoredist[intvlname] <- shored
    relativeislandwidth[intvlname] <- iwidth
  }
  shoredist$mean <- apply(shoredist[3:ncol(shoredist)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  relativeislandwidth$mean <- apply(relativeislandwidth[3:ncol(relativeislandwidth)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  utils::write.csv(shoredist,base::paste0(outdir,"/island_leastshoredist.csv"))
  utils::write.csv(relativeislandwidth,base::paste0(outdir,"/island_relativewidth.csv"))
}

#' Calculate least-cost distance between points over time
#'
#' This function calculates the least-cost overland distance between each pairwise combination of user-specified points. By assigning exposed
#' land a resistance value of 1 and submerged areas a resistance value of 999,999, this function calculates the shortest route between each pair of points
#' that minimises the number of overwater crossings. If both points lie on the same landmass, a distance value of 0 is returned. If one or both islands
#' is/are underwater, then a value of NA is returned. After calculating the pairwise least-cost distances for each interval, this function calculates
#' the time-weighted average least-cost distance between each pair of points.
#'
#' @param points A user-generated multi-point shapefile (.shp) containing at least two points. If the shapefile attribute table contains a column
#' labelled 'Name', the output distance matrix will use the identifiers in this column to label each pairwise comparison. Otherwise, the
#' output distance matrix will use the attribute FID values instead.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param intervalfile This is the master control file generated using either the getintervals_time() or getintervals_sealvl() functions, or custom generated by the user,
#' that defines the number of intervals, the sea level at each interval, and the duration of each interval.
#' @param mapdir The directory containing map outputs from the makemaps() function
#' (i.e. the directory containing the raster_flat, raster_topo, and shapefile folders).
#' @param outdir Output directory for the point-to-point least-cost distance matrix file.
#' If the specified output directory doesn't already exist, PleistoDist will create the output directory.
#' @return This function outputs a pairwise distance matrix of pairwise inter-island least shore-to-shore distances in long format, with
#' one column per interval.
#' @examples
#' \donttest{
#' #create temp directory
#' path <- file.path(tempdir())
#' #create points shapefile
#' points <- sf::st_multipoint(rbind(c(178.68408,-16.82573), c(179.13585,-16.59487)))
#' #convert points to feature geometry
#' points <- sf::st_sfc(points)
#' points <- sf::st_cast(points, "POINT")
#' #set default projection (WGS84)
#' sf::st_crs(points) <- 4326
#' #save shapefile
#' sf::write_sf(points,paste0(path,"/points.shp"))
#' #load bathymetry file
#' fiji <- system.file("extdata","FJ.asc",package="PleistoDist")
#' #generate interval file for 2 intervals and 20 kya cutoff time, binning by time
#' getintervals_time(time=20,intervals=2,outdir=path)
#' #generate maps based on interval file, projecting map using EPSG:3141
#' makemaps(inputraster=fiji,epsg=3141,intervalfile=paste0(path,"/intervals.csv"),outdir=path)
#' #calculate inter-point least cost distances, projecting points using EPSG:3141
#' pleistodist_leastcost(points=paste0(path,"/points.shp"),epsg=3141,
#'     intervalfile=paste0(path,"/intervals.csv"),mapdir=path,outdir=path)
#' }
#' @export
pleistodist_leastcost <- function(points,epsg,intervalfile,mapdir,outdir) {
  #check for existence of output directory, create output directory if it doesn't already exist.
  if(base::dir.exists(outdir)==FALSE) {
    base::dir.create(outdir)
  }
  intervalfile = utils::read.csv(intervalfile)
  points = sf::st_read(points,fid_column_name="FID")

  #check projection on input points, set to WGS84 (EPSG:4326) if undefined
  if (is.na(sf::st_crs(points)$input)) {
    sf::st_crs(points) = 4326
  }

  #reproject input points to target projection
  points_transformed <- sf::st_transform(points,epsg)
  numintervals = max(intervalfile$Interval)

  #check to see if the points contains a column of unique names, otherwise use FID numbers as identifiers
  if ("Name" %in% colnames(points_transformed) && anyDuplicated(points_transformed$Name) == 0) {
    p1_names <- points_transformed$Name
    p2_names <- points_transformed$Name

  } else {
    p1_names <- points_transformed$FID
    p2_names <- points_transformed$FID
    message("No 'Name' column with unique identifiers detected in points file, defaulting to FID values instead")
  }

  #create dataframe for storing point-to-point least cost distance matrix
  leastcostdist <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names,include.equals=T)[,2],
    Island2 = expand.grid(p1_names,p2_names,include.equals=T)[,1]
  )

  #create dataframe for storing relative island width matrix
  relativeislandwidth <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names)[,2],
    Island2 = expand.grid(p1_names,p2_names)[,1]
  )

  #for each interval...
  for (i in 0:numintervals) {
    #load interval shapefile
    invector <- sf::st_read(base::paste0(mapdir,"/shapefile/interval",i,".shp")) #load interval shapefile
    inraster <- raster::raster(base::paste0(mapdir,"/raster_flat/interval",i,".tif")) #load flat raster for least cost path calculation
    inraster[is.na(inraster[])] <- 1/9999999 #convert raster NA values (water) to very low conductance values
    trraster <- gdistance::transition(inraster,mean,directions=8) #convert inraster into a transition matrix
    trraster_corr <- gdistance::geoCorrection(trraster,type="c",scl=FALSE) #apply geographical correction to the transition matrix (for least cost path)
    intvlname <- base::paste0("interval",i)

    #create empty lists for storing distance values at each interval
    leastcostd <- c()
    iwidth <- c()
    message(base::paste0("Analysing island distances for interval ",i,"..."))
    #start of first order for-loop to select island 1
    for (f in 1:nrow(points_transformed)) {
      #start second order for-loop to select island 2
      for (t in 1:nrow(points_transformed)) {
        #extract coordinates for island 1
        p1 <- points_transformed$geometry[f]
        #extract coordinates for island 2
        p2 <- points_transformed$geometry[t]
        message(base::paste0("Calculating distance between points ",p1_names[f]," and ",p2_names[t],"... "))
        if (f == t) { #write NA for self-comparisons
          leastcostd <- c(leastcostd,NA)
          iwidth <- c(iwidth,NA)
        } else if ((nrow(invector[p1, ,op=sf::st_contains]) || nrow(invector[p2, ,op=sf::st_contains])) == 0) { #check to see if island1 or island2 is submerged under water
          leastcostd <- c(leastcostd,NA)
          iwidth <- c(iwidth,NA)
          message("One or both islands are underwater during this interval, writing a value of NA")
        } else {
          combinedpoints <- sf::st_sfc(rbind(p1,p2),crs=epsg)
          islandpair <- invector[combinedpoints, ,op=sf::st_contains]
          if (nrow(islandpair) == 1) { #check to see if both points are on the same island, if TRUE write distance of 0 for island-based comparisons
            lcdistance <- as.numeric(sp::SpatialLinesLengths(gdistance::shortestPath(trraster_corr,sf::st_coordinates(p1),sf::st_coordinates(p2),output="SpatialLines")))
            leastcostd <- c(leastcostd,lcdistance) #for least cost distance, calculate least cost path
            iwidth <- c(iwidth,0)
            message(base::paste0("Least cost distance between points is ",lcdistance,"m"))
          } else {
            island1 <- invector[p1, ,op=sf::st_contains]
            island2 <- invector[p2, ,op=sf::st_contains]
            island1_centroid <- sf::st_centroid(island1)
            island2_centroid <- sf::st_centroid(island2)
            centroidline <- sf::st_linestring(rbind(sf::st_coordinates(island1_centroid),sf::st_coordinates(island2_centroid)))

            m = -1/(stats::line(centroidline)$coefficients[2])
            c = sf::st_coordinates(island1_centroid)[,2]-(m*sf::st_coordinates(island1_centroid)[,1])

            #calculate 2 points on widthline at edges of bounding box
            x1 = sf::st_bbox(island1)$xmax
            y1 = m*x1+c

            x2 = sf::st_bbox(island1)$xmin
            y2 = m*x2+c

            coord1 <- sf::st_point(c(x1,y1))
            coord2 <- sf::st_point(c(x2,y2))
            widthline <- sf::st_sfc(sf::st_linestring(c(coord1,coord2)),crs=epsg)
            #because complex island shapes can slice widthlines into multiple fragments, convert widthline into coordinates,
            #generate a new line, and calculate length from this new line
            islandwidth<- sf::st_length(sf::st_linestring(sf::st_coordinates(sf::st_intersection(widthline,island1))))
            iwidth = c(iwidth,islandwidth)
            lcdistance <- as.numeric(sp::SpatialLinesLengths(gdistance::shortestPath(trraster_corr,sf::st_coordinates(p1),sf::st_coordinates(p2),output="SpatialLines")))
            leastcostd <- c(leastcostd,lcdistance)
            message(base::paste0("Least cost distance between points is ",lcdistance,"m"))
          }
        }
      }
    }
    leastcostdist[intvlname] <- leastcostd
    relativeislandwidth[intvlname] <- iwidth
  }
  leastcostdist$mean <- apply(leastcostdist[3:ncol(leastcostdist)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  relativeislandwidth$mean <- apply(relativeislandwidth[3:ncol(relativeislandwidth)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  utils::write.csv(leastcostdist,base::paste0(outdir,"/island_leastcostdist.csv"))
  utils::write.csv(relativeislandwidth,base::paste0(outdir,"/island_relativewidth.csv"))
}

#' Calculate inter-island mean shore-to-shore distance over time
#'
#' This function calculates the mean shore-to-shore distance between island landmasses over Pleistocene time. This function is slightly different
#' from other PleistoDist functions because the resultant distance matrix is asymmetrical. This is due to the fact that the mean shore-to-shore distance from
#' a large island to a small one will be greater relative to the reverse, since larger islands have longer and more complex shorelines, so the mean dispersal distance
#' from any point along the larger island's shoreline facing the smaller island is likely to be greater than in the reverse case. To calculate the mean
#' shore-to-shore distance, PleistoDist generates a number of points at equal intervals (1 point every 100m by default) along the shoreline of the source island at each time/sea level interval, filters out points not directly
#' facing the receiving island, and calculates the average minimum distance between each point along the source island shoreline and the shoreline of the
#' receiving island. Because the number of shoreline points can run into the thousands for very large islands, thereby increasing computational time, users
#' can set a cap on the maximum number of shoreline points to sample (the default cap is 1000 points). As with other PleistoDist functions, if both points lie on the same landmass, a distance value of 0 is returned. If one or both islands
#' is/are underwater, then a value of NA is returned. After calculating the pairwise island distances for each interval, this function calculates
#' the time-weighted average mean shore-to-shore distance between each island pair.
#'
#' @param points A user-generated multi-point shapefile (.shp) containing at least two points. If the shapefile attribute table contains a column
#' labelled 'Name', the output distance matrix will use the identifiers in this column to label each pairwise comparison. Otherwise, the
#' output distance matrix will use the attribute FID values instead.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param intervalfile This is the master control file generated using either the getintervals_time() or getintervals_sealvl() functions, or custom generated by the user,
#' that defines the number of intervals, the sea level at each interval, and the duration of each interval.
#' @param mapdir The directory containing map outputs from the makemaps() function
#' (i.e. the directory containing the raster_flat, raster_topo, and shapefile folders).
#' @param outdir Output directory for the inter-island mean shore-to-shore distance matrix file.
#' If the specified output directory doesn't already exist, PleistoDist will create the output directory.
#' @param maxsamp Maximum number of points to sample along the donor island shoreline, mainly to limit computational time. If users are
#' ok with long runtimes, then maxsamp can be set to NA. Default number of sample points is 1000.
#' @return This function outputs an asymmetric pairwise distance matrix of mean shore-to-shore distances in long format, with
#' one column per interval. Because this distance matrix is asymmetric, the directionality of the distance calculation matters.
#' @examples
#' \donttest{
#' #create temp directory
#' path <- file.path(tempdir())
#' #create points shapefile
#' points <- sf::st_multipoint(rbind(c(178.68408,-16.82573), c(179.13585,-16.59487)))
#' #convert points to feature geometry
#' points <- sf::st_sfc(points)
#' points <- sf::st_cast(points, "POINT")
#' #set default projection (WGS84)
#' sf::st_crs(points) <- 4326
#' #save shapefile
#' sf::write_sf(points,paste0(path,"/points.shp"))
#' #load bathymetry file
#' fiji <- system.file("extdata","FJ.asc",package="PleistoDist")
#' #generate interval file for 1 interval and 20 kya cutoff time, binning by time
#' getintervals_time(time=20,intervals=1,outdir=path)
#' #generate maps based on interval file, projecting map using EPSG:3141
#' makemaps(inputraster=fiji,epsg=3141,intervalfile=paste0(path,"/intervals.csv"),outdir=path)
#' #calculate inter-island mean shore-to-shore distances, projecting points using EPSG:3141,
#' #with 500 sampling points.
#' pleistodist_meanshore(points=paste0(path,"/points.shp"),epsg=3141,
#'     intervalfile=paste0(path,"/intervals.csv"),mapdir=path,outdir=path,maxsamp=500)
#' }
#' @export
pleistodist_meanshore <- function(points,epsg,intervalfile,mapdir,outdir,maxsamp=1000) {
  #check for existence of output directory, create output directory if it doesn't already exist.
  if(base::dir.exists(outdir)==FALSE) {
    base::dir.create(outdir)
  }
  intervalfile = utils::read.csv(intervalfile)
  points = sf::st_read(points,fid_column_name="FID")

  #check projection on input points, set to WGS84 (EPSG:4326) if undefined
  if (is.na(sf::st_crs(points)$input)) {
    sf::st_crs(points) = 4326
  }

  #reproject input points to target projection
  points_transformed <- sf::st_transform(points,epsg)
  numintervals = max(intervalfile$Interval)

  #check to see if the points contains a column of unique names, otherwise use FID numbers as identifiers
  if ("Name" %in% colnames(points_transformed) && anyDuplicated(points_transformed$Name) == 0) {
    p1_names <- points_transformed$Name
    p2_names <- points_transformed$Name

  } else {
    p1_names <- points_transformed$FID
    p2_names <- points_transformed$FID
    message("No 'Name' column with unique identifiers detected in points file, defaulting to FID values instead")
  }

  #create dataframe for storing mean shore distance matrix
  meanshoredist <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names)[,2],
    Island2 = expand.grid(p1_names,p2_names)[,1]
  )

  #create dataframe for storing relative island width matrix
  relativeislandwidth <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names)[,2],
    Island2 = expand.grid(p1_names,p2_names)[,1]
  )

  #for each interval...
  for (i in 0:numintervals) {
    #load interval shapefile
    invector <- sf::st_read(base::paste0(mapdir,"/shapefile/interval",i,".shp")) #load interval shapefile
    intvlname <- base::paste0("interval",i)
    meanshored <- c()
    iwidth <- c()

    #start of first order for-loop to select island 1
    for (f in 1:nrow(points_transformed)) {

      #start second order for-loop to select island 2
      for (t in 1:nrow(points_transformed)) {

        #extract coordinates for island 1
        p1 <- points_transformed$geometry[f]
        #extract coordinates for island 2
        p2 <- points_transformed$geometry[t]

        message(base::paste0("Calculating distance between points ",p1_names[f]," and ",p2_names[t],"... "))

        if (f == t) { #write NA for self-comparisons
          meanshored <- c(meanshored,NA)
          iwidth <- c(iwidth,NA)
        } else if ((nrow(invector[p1, ,op=sf::st_contains]) || nrow(invector[p2, ,op=sf::st_contains])) == 0) { #check to see if island1 or island2 is submerged under water
          meanshored <- c(meanshored,NA)
          iwidth <- c(iwidth,NA)
          message("One or both islands are underwater during this interval, writing a value of NA")
        } else {
          combinedpoints <- sf::st_sfc(rbind(p1,p2),crs=epsg)
          islandpair <- invector[combinedpoints, ,op=sf::st_contains]
          if (nrow(islandpair) == 1) { #check to see if both points are on the same island, if TRUE write distance of 0 for island-based comparisons
            meanshored <- c(meanshored,0)
            iwidth <- c(iwidth,NA)
            message("Both points are on the same island, writing a value of 0")
          } else {
            island1 <- invector[p1, ,op=sf::st_contains]
            island2 <- invector[p2, ,op=sf::st_contains]
            island1_centroid <- sf::st_centroid(island1)
            island2_centroid <- sf::st_centroid(island2)
            centroidline <- sf::st_linestring(rbind(sf::st_coordinates(island1_centroid),sf::st_coordinates(island2_centroid)))

            m = -1/(stats::line(centroidline)$coefficients[2])
            c = sf::st_coordinates(island1_centroid)[,2]-(m*sf::st_coordinates(island1_centroid)[,1])

            #calculate 2 points on widthline at edges of bounding box
            x1 = sf::st_bbox(island1)$xmax
            y1 = m*x1+c

            x2 = sf::st_bbox(island1)$xmin
            y2 = m*x2+c

            coord1 <- sf::st_point(c(x1,y1))
            coord2 <- sf::st_point(c(x2,y2))
            widthline <- sf::st_sfc(sf::st_linestring(c(coord1,coord2)),crs=epsg)
            #because complex island shapes can slice widthlines into multiple fragments, convert widthline into coordinates,
            #generate a new line, and calculate length from this new line
            islandwidth<- sf::st_length(sf::st_linestring(sf::st_coordinates(sf::st_intersection(widthline,island1))))
            iwidth = c(iwidth,islandwidth)

            island1_shoreline <- sf::st_cast(island1,"LINESTRING")
            #pick the line segment that is closest to the centroid of the receiving island (since there may be inland lakes)
            lineseg <- which.min(sf::st_distance(island1_shoreline,island2))
            island1_shoreline <- island1_shoreline$geometry[lineseg]

            #generate points at 100m intervals along shoreline
            island1_shorepoints <- sf::st_line_sample(island1_shoreline,density=1/100)
            island1_shorepoints <- sf::st_cast(island1_shorepoints,"POINT")

            #downsample the number of shorepoints if the island is very big (to speed up computations)
            if (length(island1_shorepoints) > maxsamp && !is.na(maxsamp)) {
              island1_shorepoints <- sf::st_line_sample(island1_shoreline,n=maxsamp)
              island1_shorepoints <- sf::st_cast(island1_shorepoints,"POINT")
            }

            #draw a line between each shorepoint and the receiving island, and exclude points with lines crossing the shore multiple times (indicating that it faces away from the receiving island)
            island1_shorepoints_filter1 <- sf::st_crosses(sf::st_nearest_points(island1_shorepoints,island2),island1_shoreline,sparse=FALSE)
            #subset points that face receiving island
            island1_shorepoints_filtered <- island1_shorepoints[which(island1_shorepoints_filter1==FALSE)]
            #apply a second filter to eliminate points that cross the island more than once
            island1_shorepoints_filter2 <- sf::st_overlaps(sf::st_nearest_points(island1_shorepoints_filtered,island2),island1_shoreline,sparse=FALSE)
            island1_shorepoints_final <- island1_shorepoints_filtered[which(island1_shorepoints_filter2==FALSE)]

            msd <- as.numeric(mean(sf::st_distance(island1_shorepoints_final,island2)))
            meanshored <- c(meanshored,msd)
            #meanshored_stdev <- c(meanshored_stdev,as.numeric(sd(sf::st_distance(island1_shorepoints,islandpair$geometry[2]))))
            message(paste0("The mean shore-to-shore distance between islands is ",msd,"m"))
          }
        }
      }
    }
    meanshoredist[intvlname] <- meanshored
    relativeislandwidth[intvlname] <- iwidth
  }
  meanshoredist$mean <- apply(meanshoredist[3:ncol(meanshoredist)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  relativeislandwidth$mean <- apply(relativeislandwidth[3:ncol(relativeislandwidth)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  utils::write.csv(meanshoredist,base::paste0(outdir,"/island_meanshoredist.csv"))
  utils::write.csv(relativeislandwidth,base::paste0(outdir,"/island_relativewidth.csv"))
}

#' Calculate width of islands relative to each other
#'
#' This function calculates the width of source islands relative to their recipient islands for every pairwise combination of islands as specified
#' by the user-defined source points. The island width is calculated by drawing a line between the centroids of the two islands, and plotting
#' a line perpendicular to the centroid-centroid line extending from the centroid of the source island to the island edges. This widthline changes
#' relative to the orientation of the receiving island, and is mostly needed for calculating the expected net inter-island migration, as defined
#' in chapter 6 of MacArthur and Wilson (1967).
#' @param points A user-generated multi-point shapefile (.shp) containing at least two points. If the shapefile attribute table contains a column
#' labelled 'Name', the output distance matrix will use the identifiers in this column to label each pairwise comparison. Otherwise, the
#' output distance matrix will use the attribute FID values instead.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param intervalfile This is the master control file generated using either the getintervals_time() or getintervals_sealvl()
#' function that defines the number of intervals, the sea level at each interval, and the duration of each interval. By default, this
#' function will use the "intervals.csv" file stored in the output folder, but users can also specify their own custom interval
#' file (with nice round mean sea level values, for example), although users need to ensure that the same column names are preserved, and
#' be aware that custom interval files may lead to inaccurate weighted mean distance calculations.
#' @param mapdir The directory containing map outputs from the makemaps() function
#' (i.e. the directory containing the raster_flat, raster_topo, and shapefile folders).
#' @param outdir Output directory for the inter-island relative width matrix file.
#' If the specified output directory doesn't already exist, PleistoDist will create the output directory.
#' @return This function outputs an asymmetric pairwise distance matrix of relative island widths in long format, with
#' one column per interval. Because this distance matrix is asymmetric, the directionality of the distance calculation matters.
#' @examples
#' \donttest{
#' #create temp directory
#' path <- file.path(tempdir())
#' #create points shapefile
#' points <- sf::st_multipoint(rbind(c(179.41256,-17.79426), c(179.27600,-17.97850)))
#' #convert points to feature geometry
#' points <- sf::st_sfc(points)
#' points <- sf::st_cast(points, "POINT")
#' #set default projection (WGS84)
#' sf::st_crs(points) <- 4326
#' #save shapefile
#' sf::write_sf(points,paste0(path,"/points.shp"))
#' #load bathymetry file
#' fiji <- system.file("extdata","FJ.asc",package="PleistoDist")
#' #generate interval file for 2 intervals and 20 kya cutoff time, binning by time
#' getintervals_time(time=20,intervals=2,outdir=path)
#' #generate maps based on interval file, projecting map using EPSG:3141
#' makemaps(inputraster=fiji,epsg=3141,intervalfile=paste0(path,"/intervals.csv"),outdir=path)
#' #calculate relative source island width, projecting points using EPSG:3141
#' pleistodist_relativewidth(points=paste0(path,"/points.shp"),epsg=3141,
#'     intervalfile=paste0(path,"/intervals.csv"),mapdir=path,outdir=path)
#' }
#' @export
pleistodist_relativewidth <- function (points, epsg,intervalfile,mapdir,outdir) {
  #check for existence of output directory, create output directory if it doesn't already exist.
  if(base::dir.exists(outdir)==FALSE) {
    base::dir.create(outdir)
  }
  intervalfile = utils::read.csv(intervalfile)
  points = sf::st_read(points,fid_column_name="FID")

  #check projection on input points, set to WGS84 (EPSG:4326) if undefined
  if (is.na(sf::st_crs(points)$input)) {
    sf::st_crs(points) = 4326
  }

  #reproject input points to target projection
  points_transformed <- sf::st_transform(points,epsg)
  numintervals = max(intervalfile$Interval)

  #check to see if the points contains a column of unique names, otherwise use FID numbers as identifiers
  if ("Name" %in% colnames(points_transformed) && anyDuplicated(points_transformed$Name) == 0) {
    p1_names <- points_transformed$Name
    p2_names <- points_transformed$Name
  } else {
    p1_names <- points_transformed$FID
    p2_names <- points_transformed$FID
    message("No 'Name' column with unique identifiers detected in points file, defaulting to FID values instead")
  }

  #create dataframe for storing relative island width matrix
  relativeislandwidth <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names)[,2],
    Island2 = expand.grid(p1_names,p2_names)[,1]
  )

  #for each interval...
  for (i in 0:numintervals) {
    #load interval shapefile
    invector <- sf::st_read(base::paste0(mapdir,"/shapefile/interval",i,".shp")) #load interval shapefile
    intvlname <- base::paste0("interval",i)

    #create empty lists for storing distance values at each interval
    iwidth <- c()

    message(paste0("Analysing island width for interval ",i,"..."))
    #start of first order for-loop to select island 1
    for (f in 1:nrow(points_transformed)) {

      #start second order for-loop to select island 2
      for (t in 1:nrow(points_transformed)) {

        #extract coordinates for island 1
        p1 <- points_transformed$geometry[f]
        #extract coordinates for island 2
        p2 <- points_transformed$geometry[t]

        message(paste0("Calculating distance between points ",p1_names[f]," and ",p2_names[t],"... "))

        if (f == t) { #write NA for self-comparisons
          iwidth <- c(iwidth,NA)
        } else if ((nrow(invector[p1, ,op=sf::st_contains]) || nrow(invector[p2, ,op=sf::st_contains])) == 0) { #check to see if island1 or island2 is submerged under water
          iwidth <- c(iwidth,NA)
          message("One or both islands are underwater during this interval, writing a value of NA")
        } else {
          combinedpoints <- sf::st_sfc(rbind(p1,p2),crs=epsg)
          islandpair <- invector[combinedpoints, ,op=sf::st_contains]
          if (nrow(islandpair) == 1) { #check to see if both points are on the same island, if TRUE write distance of 0 for island-based comparisons
            iwidth <- c(iwidth,NA)
            message("Both points are on the same island, writing a value of 0 for centroid-to-centroid distance")
          } else {
            island1 <- invector[p1, ,op=sf::st_contains]
            island2 <- invector[p2, ,op=sf::st_contains]
            island1_centroid <- sf::st_centroid(island1)
            island2_centroid <- sf::st_centroid(island2)
            centroidline <- sf::st_linestring(rbind(sf::st_coordinates(island1_centroid),sf::st_coordinates(island2_centroid)))

            m = -1/(stats::line(centroidline)$coefficients[2])
            c = sf::st_coordinates(island1_centroid)[,2]-(m*sf::st_coordinates(island1_centroid)[,1])

            #calculate 2 points on widthline at edges of bounding box
            x1 = sf::st_bbox(island1)$xmax
            y1 = m*x1+c

            x2 = sf::st_bbox(island1)$xmin
            y2 = m*x2+c

            coord1 <- sf::st_point(c(x1,y1))
            coord2 <- sf::st_point(c(x2,y2))
            widthline <- sf::st_sfc(sf::st_linestring(c(coord1,coord2)),crs=epsg)
            #because complex island shapes can slice widthlines into multiple fragments, convert widthline into coordinates,
            #generate a new line, and calculate length from this new line
            islandwidth<- sf::st_length(sf::st_linestring(sf::st_coordinates(sf::st_intersection(widthline,island1))))
            iwidth = c(iwidth,islandwidth)
          }
        }
      }
    }
    relativeislandwidth[intvlname] <- iwidth
  }
  relativeislandwidth$mean <- apply(relativeislandwidth[3:ncol(relativeislandwidth)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  utils::write.csv(relativeislandwidth,base::paste0(outdir,"/island_relativewidth.csv"))
}

#' Calculate area, perimeter and surface area of islands over time
#'
#' This function calculates the area, perimeter, and surface area of user-specified islands at each time/sea level interval, and
#' performs the same calculations as the pleistoshape_area(), pleistoshape_perimeter, and pleistoshape_surfacearea() functions combined.
#' This function is nonetheless provided since it allows users to calculate all three shape metrics simultaneously, which is faster than running
#' each constituent function individually. PleistoDist will return a value of NA if the island is submerged underwater for that particular interval. After calculating shape metrics for each interval, this function calculates
#' the time-weighted average area, perimeter, and surface area of each selected island.
#' @param points A user-generated multi-point shapefile (.shp) containing at least two points. If the shapefile attribute table contains a column
#' labelled 'Name', the output distance matrix will use the identifiers in this column to label each pairwise comparison. Otherwise, the
#' output distance matrix will use the attribute FID values instead.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param intervalfile This is the master control file generated using either the getintervals_time() or getintervals_sealvl()
#' function that defines the number of intervals, the sea level at each interval, and the duration of each interval. By default, this
#' function will use the "intervals.csv" file stored in the output folder, but users can also specify their own custom interval
#' file (with nice round mean sea level values, for example), although users need to ensure that the same column names are preserved, and
#' be aware that custom interval files may lead to inaccurate weighted mean calculations.
#' @param mapdir The directory containing map outputs from the makemaps() function
#' (i.e. the directory containing the raster_flat, raster_topo, and shapefile folders).
#' @param outdir Output directory for the island area, perimeter, and surface area matrix files.
#' If the specified output directory doesn't already exist, PleistoDist will create the output directory.
#' @return This function outputs three matrices of island shape estimates in long format, with one column per interval in each matrix.
#' @examples
#' \donttest{
#' #create temp directory
#' path <- file.path(tempdir())
#' #create points shapefile
#' points <- sf::st_multipoint(rbind(c(179.41256,-17.79426), c(179.27600,-17.97850)))
#' #convert points to feature geometry
#' points <- sf::st_sfc(points)
#' points <- sf::st_cast(points, "POINT")
#' #set default projection (WGS84)
#' sf::st_crs(points) <- 4326
#' #save shapefile
#' sf::write_sf(points,paste0(path,"/points.shp"))
#' #load bathymetry file
#' fiji <- system.file("extdata","FJ.asc",package="PleistoDist")
#' #generate interval file for 1 interval and 20 kya cutoff time, binning by time
#' getintervals_time(time=20,intervals=1,outdir=path)
#' #generate maps based on interval file, projecting map using EPSG:3141
#' makemaps(inputraster=fiji,epsg=3141,intervalfile=paste0(path,"/intervals.csv"),outdir=path)
#' #calculate island shape, perimeter, and surface area, projecting points using EPSG:3141
#' pleistoshape_all(points=paste0(path,"/points.shp"),epsg=3141,
#'     intervalfile=paste0(path,"/intervals.csv"),mapdir=path,outdir=path)
#' }
#'
#' @export
pleistoshape_all <- function (points,epsg,intervalfile,mapdir,outdir) {
  #check for existence of output directory, create output directory if it doesn't already exist.
  if(base::dir.exists(outdir)==FALSE) {
    base::dir.create(outdir)
  }
  intervalfile <- utils::read.csv(intervalfile)
  points = sf::st_read(points,fid_column_name="FID")

  #check projection on input points, set to WGS84 (EPSG:4326) if undefined
  if (is.na(sf::st_crs(points)$input)) {
    sf::st_crs(points) = 4326
  }

  #reproject input points to target projection
  points_transformed <- sf::st_transform(points,epsg)

  numintervals = max(intervalfile$Interval)

  #check to see if the points contains a column of unique names, otherwise use FID numbers as identifiers
  if ("Name" %in% colnames(points_transformed) && anyDuplicated(points_transformed$Name) == 0) {
    p1_names <- points_transformed$Name
  } else {
    p1_names <- points_transformed$FID
    message("No 'Name' column with unique identifiers detected in points file, defaulting to FID values instead")
  }

  #create dataframe for storing island area
  islandarea <- tibble::tibble(
    Island = p1_names
  )
  #create dataframe for storing island perimeter
  islandperimeter <- tibble::tibble(
    Island = p1_names
  )
  #create dataframe for storing island perimeter
  islandsurfacearea <- tibble::tibble(
    Island = p1_names
  )

  #for each interval...
  for (i in 0:numintervals) {
    #load interval shapefile
    invector <- sf::st_read(base::paste0(mapdir,"/shapefile/interval",i,".shp"))
    inraster <- terra::rast(base::paste0(mapdir,"/raster_topo/interval",i,".tif"))-intervalfile$MeanDepth[i+1]
    intvlname <- base::paste0("interval",i)
    ia <- c() #create blank list for temporarily storing area values (before offloading into the main dataframe)
    ip <- c() #create blank list for temporarily storing perimeter values (before offloading into the main dataframe)
    isa <- c() #create blank list for temporarily storing surface area values (before offloading into the main dataframe)

    #start of first order for-loop to select island 1
    for (f in 1:nrow(points_transformed)) {

      #extract coordinates for island
      p1 <- points_transformed$geometry[f]

      message("Calculating island area and perimeter for point ",p1_names[f],sep="")

      if (nrow(invector[p1, ,op=sf::st_contains]) == 0) { #check to see if island is submerged under water
        ia <- c(ia,0)
        ip <- c(ip,0)
        isa <- c(isa,0)
        message("Island is underwater during this interval, writing a value of 0")
      } else {
        iarea <- as.numeric(sf::st_area(invector[p1, ,op=sf::st_contains]))
        iperimeter <- as.numeric(lwgeom::st_perimeter(invector[p1, ,op=sf::st_contains]))
        islandmask <- terra::vect(invector[p1, ,op=sf::st_contains]) #use shapefile to create a mask layer
        islandtopo <- terra::mask(inraster,islandmask) #mask topo raster to isolate single island
        islandtopo_cropped <- terra::crop(islandtopo,islandmask) #crop raster extent to single island
        islandtopo_cropped_sp <- methods::as(raster::raster(islandtopo_cropped),"SpatialPixelsDataFrame") #convert spatraster to sp format
        isurfacearea <- sp::surfaceArea(islandtopo_cropped_sp) #calculate surface area (incorporating elevation)
        ia <- c(ia,iarea)
        ip <- c(ip,iperimeter)
        isa <- c(isa,isurfacearea)
        message(paste0("Island perimeter is ",iperimeter,"m"))
        message(paste0("Island area is ",iarea,"m^2"))
        message(paste0("Island surface area is ",isurfacearea,"m^2"))
      }
    }
    islandarea[intvlname] <- ia
    islandperimeter[intvlname] <- ip
    islandsurfacearea[intvlname] <- isa
  }
  islandarea$mean <- apply(islandarea[2:ncol(islandarea)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  islandperimeter$mean <- apply(islandperimeter[2:ncol(islandperimeter)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  islandsurfacearea$mean <- apply(islandsurfacearea[2:ncol(islandsurfacearea)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  utils::write.csv(islandarea,base::paste0(outdir,"/island_area.csv"))
  utils::write.csv(islandperimeter,base::paste0(outdir,"/island_perimeter.csv"))
  utils::write.csv(islandsurfacearea,base::paste0(outdir,"/island_surfacearea.csv"))
}

#' Calculate island perimeter over time
#'
#' This function calculates the perimeter of user-specified islands at each time/sea level interval, and calculates the weighted average perimter over the user-defined time period.
#' PleistoDist will return a value of NA if the island is submerged underwater for that particular interval.
#' @param points A user-generated multi-point shapefile (.shp) containing at least two points. If the shapefile attribute table contains a column
#' labelled 'Name', the output distance matrix will use the identifiers in this column to label each pairwise comparison. Otherwise, the
#' output distance matrix will use the attribute FID values instead.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param intervalfile This is the master control file generated using either the getintervals_time() or getintervals_sealvl()
#' function that defines the number of intervals, the sea level at each interval, and the duration of each interval. By default, this
#' function will use the "intervals.csv" file stored in the output folder, but users can also specify their own custom interval
#' file (with nice round mean sea level values, for example), although users need to ensure that the same column names are preserved, and
#' be aware that custom interval files may lead to inaccurate weighted mean calculations.
#' @param mapdir The directory containing map outputs from the makemaps() function
#' (i.e. the directory containing the raster_flat, raster_topo, and shapefile folders).
#' @param outdir Output directory for the island perimeter matrix file.
#' If the specified output directory doesn't already exist, PleistoDist will create the output directory.
#' @return This function outputs a matrix of island perimeter estimates in long format, with one column per interval.
#' @examples
#' #create temp directory
#' path <- file.path(tempdir())
#' #create points shapefile
#' points <- sf::st_multipoint(rbind(c(179.41256,-17.79426), c(179.27600,-17.97850)))
#' #convert points to feature geometry
#' points <- sf::st_sfc(points)
#' points <- sf::st_cast(points, "POINT")
#' #set default projection (WGS84)
#' sf::st_crs(points) <- 4326
#' #save shapefile
#' sf::write_sf(points,paste0(path,"/points.shp"))
#' #load bathymetry file
#' fiji <- system.file("extdata","FJ.asc",package="PleistoDist")
#' #generate interval file for 2 intervals and 20 kya cutoff time, binning by time
#' getintervals_time(time=20,intervals=2,outdir=path)
#' #generate maps based on interval file, projecting map using EPSG:3141
#' makemaps(inputraster=fiji,epsg=3141,intervalfile=paste0(path,"/intervals.csv"),outdir=path)
#' #calculate island perimeter, projecting points using EPSG:3141
#' pleistoshape_perimeter(points=paste0(path,"/points.shp"),epsg=3141,
#'     intervalfile=paste0(path,"/intervals.csv"),mapdir=path,outdir=path)
#' @export
pleistoshape_perimeter <- function (points,epsg,intervalfile,mapdir,outdir) {
  #check for existence of output directory, create output directory if it doesn't already exist.
  if(base::dir.exists(outdir)==FALSE) {
    base::dir.create(outdir)
  }
  intervalfile <- utils::read.csv(intervalfile)
  points = sf::st_read(points,fid_column_name="FID")

  #check projection on input points, set to WGS84 (EPSG:4326) if undefined
  if (is.na(sf::st_crs(points)$input)) {
    sf::st_crs(points) = 4326
  }

  #reproject input points to target projection
  points_transformed <- sf::st_transform(points,epsg)
  numintervals = max(intervalfile$Interval)

  #check to see if the points contains a column of unique names, otherwise use FID numbers as identifiers
  if ("Name" %in% colnames(points_transformed) && anyDuplicated(points_transformed$Name) == 0) {
    p1_names <- points_transformed$Name
  } else {
    p1_names <- points_transformed$FID
    message("No 'Name' column with unique identifiers detected in points file, defaulting to FID values instead")
  }
  #create dataframe for storing island perimeter
  islandperimeter <- tibble::tibble(
    Island = p1_names
  )
  #for each interval...
  for (i in 0:numintervals) {
    #load interval shapefile
    invector <- sf::st_read(base::paste0(mapdir,"/shapefile/interval",i,".shp"))
    intvlname <- base::paste0("interval",i)
    ip <- c() #create blank list for temporarily storing perimeter values (before offloading into the main dataframe)

    #start of first order for-loop to select island 1
    for (f in 1:nrow(points_transformed)) {

      #extract coordinates for island
      p1 <- points_transformed$geometry[f]

      message("Calculating island area and perimeter for point ",p1_names[f],sep="")

      if (nrow(invector[p1, ,op=sf::st_contains]) == 0) { #check to see if island is submerged under water
        ip <- c(ip,0)
        message("Island is underwater during this interval, writing a value of 0")
      } else {
        iperimeter <- as.numeric(lwgeom::st_perimeter(invector[p1, ,op=sf::st_contains]))
        ip <- c(ip,iperimeter)
        message(base::paste0("Island perimeter is ",iperimeter,"m"))
      }
    }
    islandperimeter[intvlname] <- ip
  }
  islandperimeter$mean <- apply(islandperimeter[2:ncol(islandperimeter)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  utils::write.csv(islandperimeter,base::paste0(outdir,"/island_perimeter.csv"))
}

#' Calculate island area over time
#'
#' This function calculates the 2-dimensional area of user-specified islands at each time/sea level interval, and calculates the weighted average area over the user-defined time period.
#' PleistoDist will return a value of NA if the island is submerged underwater for that particular interval.
#' @param points A user-generated multi-point shapefile (.shp) containing at least two points. If the shapefile attribute table contains a column
#' labelled 'Name', the output distance matrix will use the identifiers in this column to label each pairwise comparison. Otherwise, the
#' output distance matrix will use the attribute FID values instead.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param intervalfile This is the master control file generated using either the getintervals_time() or getintervals_sealvl()
#' function that defines the number of intervals, the sea level at each interval, and the duration of each interval. By default, this
#' function will use the "intervals.csv" file stored in the output folder, but users can also specify their own custom interval
#' file (with nice round mean sea level values, for example), although users need to ensure that the same column names are preserved, and
#' be aware that custom interval files may lead to inaccurate weighted mean calculations.
#' @param mapdir The directory containing map outputs from the makemaps() function
#' (i.e. the directory containing the raster_flat, raster_topo, and shapefile folders).
#' @param outdir Output directory for the 2D island area matrix file.
#' If the specified output directory doesn't already exist, PleistoDist will create the output directory.
#' @return This function outputs a matrix of 2D island area estimates in long format, with one column per interval.
#' @examples
#' #create temp directory
#' path <- file.path(tempdir())
#' #create points shapefile
#' points <- sf::st_multipoint(rbind(c(179.41256,-17.79426), c(179.27600,-17.97850)))
#' #convert points to feature geometry
#' points <- sf::st_sfc(points)
#' points <- sf::st_cast(points, "POINT")
#' #set default projection (WGS84)
#' sf::st_crs(points) <- 4326
#' #save shapefile
#' sf::write_sf(points,paste0(path,"/points.shp"))
#' #load bathymetry file
#' fiji <- system.file("extdata","FJ.asc",package="PleistoDist")
#' #generate interval file for 2 intervals and 20 kya cutoff time, binning by time
#' getintervals_time(time=20,intervals=2,outdir=path)
#' #generate maps based on interval file, projecting map using EPSG:3141
#' makemaps(inputraster=fiji,epsg=3141,intervalfile=paste0(path,"/intervals.csv"),outdir=path)
#' #calculate 2D island area, projecting points using EPSG:3141
#' pleistoshape_area(points=paste0(path,"/points.shp"),epsg=3141,
#'     intervalfile=paste0(path,"/intervals.csv"),mapdir=path,outdir=path)
#' @export
pleistoshape_area <- function (points,epsg,intervalfile,mapdir,outdir) {
  #check for existence of output directory, create output directory if it doesn't already exist.
  if(base::dir.exists(outdir)==FALSE) {
    base::dir.create(outdir)
  }
  intervalfile <- utils::read.csv(intervalfile)
  points = sf::st_read(points,fid_column_name="FID")

  #check projection on input points, set to WGS84 (EPSG:4326) if undefined
  if (is.na(sf::st_crs(points)$input)) {
    sf::st_crs(points) = 4326
  }

  #reproject input points to target projection
  points_transformed <- sf::st_transform(points,epsg)
  numintervals = max(intervalfile$Interval)

  #check to see if the points contains a column of unique names, otherwise use FID numbers as identifiers
  if ("Name" %in% colnames(points_transformed) && anyDuplicated(points_transformed$Name) == 0) {
    p1_names <- points_transformed$Name
  } else {
    p1_names <- points_transformed$FID
    message("No 'Name' column with unique identifiers detected in points file, defaulting to FID values instead")
  }

  #create dataframe for storing island area
  islandarea <- tibble::tibble(
    Island = p1_names
  )

  #for each interval...
  for (i in 0:numintervals) {
    #load interval shapefile
    invector <- sf::st_read(base::paste0(mapdir,"/shapefile/interval",i,".shp"))
    intvlname <- base::paste0("interval",i)
    ia <- c() #create blank list for temporarily storing area values (before offloading into the main dataframe)

    #start of first order for-loop to select island 1
    for (f in 1:nrow(points_transformed)) {
      #extract coordinates for island
      p1 <- points_transformed$geometry[f]
      message("Calculating island area of ",p1_names[f],sep="")

      if (nrow(invector[p1, ,op=sf::st_contains]) == 0) { #check to see if island is submerged under water
        ia <- c(ia,0)
        message("Island is underwater during this interval, writing a value of 0")
      } else {
        iarea <- as.numeric(sf::st_area(invector[p1, ,op=sf::st_contains]))
        ia <- c(ia,iarea)
        message(base::paste0("Island area is ",iarea,"m^2"))
      }
    }
    islandarea[intvlname] <- ia
  }
  islandarea$mean <- apply(islandarea[2:ncol(islandarea)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  utils::write.csv(islandarea,base::paste0(outdir,"/island_area.csv"))
}

#' Calculate island surface area over time
#'
#' This function calculates the surface area of user-specified islands at each time/sea level interval, and calculates the weighted average surface area over the user-defined time period.
#' PleistoDist will return a value of NA if the island is submerged underwater for that particular interval.
#' @param points A user-generated multi-point shapefile (.shp) containing at least two points. If the shapefile attribute table contains a column
#' labelled 'Name', the output distance matrix will use the identifiers in this column to label each pairwise comparison. Otherwise, the
#' output distance matrix will use the attribute FID values instead.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param intervalfile This is the master control file generated using either the getintervals_time() or getintervals_sealvl()
#' function that defines the number of intervals, the sea level at each interval, and the duration of each interval. By default, this
#' function will use the "intervals.csv" file stored in the output folder, but users can also specify their own custom interval
#' file (with nice round mean sea level values, for example), although users need to ensure that the same column names are preserved, and
#' be aware that custom interval files may lead to inaccurate weighted mean calculations.
#' @param mapdir The directory containing map outputs from the makemaps() function
#' (i.e. the directory containing the raster_flat, raster_topo, and shapefile folders).
#' @param outdir Output directory for the island surface area matrix file.
#' If the specified output directory doesn't already exist, PleistoDist will create the output directory.
#' @return This function outputs a matrix of island surface area estimates in long format, with one column per interval.
#' @examples
#' \donttest{
#' #create temp directory
#' path <- file.path(tempdir())
#' #create points shapefile
#' points <- sf::st_multipoint(rbind(c(179.41256,-17.79426), c(179.27600,-17.97850)))
#' #convert points to feature geometry
#' points <- sf::st_sfc(points)
#' points <- sf::st_cast(points, "POINT")
#' #set default projection (WGS84)
#' sf::st_crs(points) <- 4326
#' #save shapefile
#' sf::write_sf(points,paste0(path,"/points.shp"))
#' #load bathymetry file
#' fiji <- system.file("extdata","FJ.asc",package="PleistoDist")
#' #generate interval file for 2 intervals and 20 kya cutoff time, binning by time
#' getintervals_time(time=20,intervals=2,outdir=path)
#' #generate maps based on interval file, projecting map using EPSG:3141
#' makemaps(inputraster=fiji,epsg=3141,intervalfile=paste0(path,"/intervals.csv"),outdir=path)
#' #calculate island surface area, projecting points using EPSG:3141
#' pleistoshape_surfacearea(points=paste0(path,"/points.shp"),epsg=3141,
#'     intervalfile=paste0(path,"/intervals.csv"),mapdir=path,outdir=path)
#' }
#' @export
pleistoshape_surfacearea <- function (points,epsg,intervalfile,mapdir,outdir) {
  #check for existence of output directory, create output directory if it doesn't already exist.
  if(base::dir.exists(outdir)==FALSE) {
    base::dir.create(outdir)
  }
  intervalfile <- utils::read.csv(intervalfile)
  points = sf::st_read(points,fid_column_name="FID")

  #check projection on input points, set to WGS84 (EPSG:4326) if undefined
  if (is.na(sf::st_crs(points)$input)) {
    sf::st_crs(points) = 4326
  }

  #reproject input points to target projection
  points_transformed <- sf::st_transform(points,epsg)
  numintervals = max(intervalfile$Interval)

  #check to see if the points contains a column of unique names, otherwise use FID numbers as identifiers
  if ("Name" %in% colnames(points_transformed) && anyDuplicated(points_transformed$Name) == 0) {
    p1_names <- points_transformed$Name
  } else {
    p1_names <- points_transformed$FID
    message("No 'Name' column with unique identifiers detected in points file, defaulting to FID values instead")
  }

  #create dataframe for storing island perimeter
  islandsurfacearea <- tibble::tibble(
    Island = p1_names
  )

  #for each interval...
  for (i in 0:numintervals) {
    #load interval shapefile
    invector <- sf::st_read(base::paste0(mapdir,"/shapefile/interval",i,".shp"))
    inraster <- terra::rast(base::paste0(mapdir,"/raster_topo/interval",i,".tif"))-intervalfile$MeanDepth[i+1]
    intvlname <- base::paste0("interval",i)
    isa <- c() #create blank list for temporarily storing surface area values (before offloading into the main dataframe)

    #start of first order for-loop to select island 1
    for (f in 1:nrow(points_transformed)) {

      #extract coordinates for island
      p1 <- points_transformed$geometry[f]
      message("Calculating island surface area for point ",p1_names[f],sep="")

      if (nrow(invector[p1, ,op=sf::st_contains]) == 0) { #check to see if island is submerged under water
        isa <- c(isa,0)
        message("Island is underwater during this interval, writing a value of 0")
      } else {
        islandmask <- terra::vect(invector[p1, ,op=sf::st_contains]) #use shapefile to create a mask layer
        islandtopo <- terra::mask(inraster,islandmask) #mask topo raster to isolate single island
        islandtopo_cropped <- terra::crop(islandtopo,islandmask) #crop raster extent to single island
        islandtopo_cropped_sp <- methods::as(raster::raster(islandtopo_cropped),"SpatialPixelsDataFrame") #convert spatraster to sp format
        isurfacearea <- sp::surfaceArea(islandtopo_cropped_sp) #calculate surface area (incorporating elevation)
        isa <- c(isa,isurfacearea)
        message(paste0("Island surface area is ",isurfacearea,"m^2"))
      }
    }
    islandsurfacearea[intvlname] <- isa
  }
  islandsurfacearea$mean <- apply(islandsurfacearea[2:ncol(islandsurfacearea)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  utils::write.csv(islandsurfacearea,base::paste0(outdir,"/island_surfacearea.csv"))
}

#' Estimate net inter-island migration over time
#'
#' This function makes use of the equation from chapter 6 of MacArthur & Wilson (1967) that estimates the number of propagules from a
#' source island arriving at a receiving island as a function of inter-island distance, area, and relative width. Although estimating
#' the unidirectional rate of migration from one island to another requires several other coefficients (alpha: the coefficient of population density,
#' and lambda: the coefficient of overwater dispersal), calculating the ratio of migration rate for bidirectional migration causes
#' these coefficients to cancel out, allowing for net migration rate to be estimated using purely geomorphological and distance-based
#' metrics. The net migration rate can therefore be calculated using the equation (atan(t2/(2*d))*A1)/(atan(t1/(2*d))*A2), where t is
#' the island width relative to the other island, d is the inter-island distance, and A is the island area, for islands 1 and 2. If the
#' net migration rate is greater than 1, then the net movement of propagules occurs from island 1 to island 2. If the net migration rate
#' is less than 1, then the net movement of propagules occurs from island 2 to island 1.
#' @param points A user-generated multi-point shapefile (.shp) containing at least two points. If the shapefile attribute table contains a column
#' labelled 'Name', the output distance matrix will use the identifiers in this column to label each pairwise comparison. Otherwise, the
#' output distance matrix will use the attribute FID values instead.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param disttype The type of inter-island distance matrix to use for this calculation, either "centroid", "leastshore", or "meanshore".
#' @param intervalfile This is the master control file generated using either the getintervals_time() or getintervals_sealvl()
#' function that defines the number of intervals, the sea level at each interval, and the duration of each interval. By default, this
#' function will use the "intervals.csv" file stored in the output folder, but users can also specify their own custom interval
#' file (with nice round mean sea level values, for example), although users need to ensure that the same column names are preserved, and
#' be aware that custom interval files may lead to inaccurate weighted mean calculations.
#' @param mapdir The directory containing map outputs from the makemaps() function
#' (i.e. the directory containing the raster_flat, raster_topo, and shapefile folders).
#' @param outdir The directory containing inter-island distance, island area, and relative island width matrices,
#' either generated by other PleistoDist functions, or generated de novo by this function if these files don't already exist.
#' All pleistodist_netmig() outputs will also be written to this folder.
#' @return This function outputs a matrix of island surface area estimates in long format, with one column per interval.
#' @examples
#' \donttest{
#' #create temp directory
#' path <- file.path(tempdir())
#' #create points shapefile
#' points <- sf::st_multipoint(rbind(c(179.41256,-17.79426), c(179.27600,-17.97850)))
#' #convert points to feature geometry
#' points <- sf::st_sfc(points)
#' points <- sf::st_cast(points, "POINT")
#' #set default projection (WGS84)
#' sf::st_crs(points) <- 4326
#' #save shapefile
#' sf::write_sf(points,paste0(path,"/points.shp"))
#' #load bathymetry file
#' fiji <- system.file("extdata","FJ.asc",package="PleistoDist")
#' #generate interval file for 2 intervals and 20 kya cutoff time, binning by time
#' getintervals_time(time=20,intervals=2,outdir=path)
#' #generate maps based on interval file, projecting map using EPSG:3141
#' makemaps(inputraster=fiji,epsg=3141,intervalfile=paste0(path,"/intervals.csv"),outdir=path)
#' #calculate net migration, for centroid-to-centroid distances, projecting points using EPSG:3141
#' pleistodist_netmig(points=paste0(path,"/points.shp"),epsg=3141,disttype="centroid",
#'     intervalfile=paste0(path,"/intervals.csv"),mapdir=path,outdir=path)
#' }
#' @importFrom dplyr %>%
#' @export
pleistodist_netmig <- function(points,epsg,disttype,intervalfile,mapdir,outdir) {
  #check to make sure a valid disttype is supplied
  if (disttype %in% c("centroid","leastshore","meanshore")==FALSE) {
    stop("Please enter a valid distance matrix type (centroid, leastshore, or meanshore)")
  }
  #parse disttype input, reading the specified distance matrix file from the output folder
  if (file.exists(base::paste0(outdir,"/island_",disttype,"dist.csv"))==TRUE) {
    distfile = utils::read.csv(base::paste0(outdir,"/island_",disttype,"dist.csv"))
    message(base::paste0("Reading distance matrix: ",outdir,"/island_",disttype,"dist.csv"))
  } else {
    #if PleistoDist doesn't detect the correct distance matrix, then call the associated function to generate the distance matrix
    message("No distance matrix file found in output folder, generating distance matrix from scratch (this might take a while...)")
    do.call(base::paste0("pleistodist_",disttype),list(points,epsg,intervalfile,mapdir,outdir))
    distfile <- utils::read.csv(base::paste0(outdir,"/island_",disttype,"dist.csv"))
  }
  widthfile <- utils::read.csv(base::paste0(outdir,"/island_relativewidth.csv")) #read relative island widths
  if ((disttype != "euclidean") && ncol(widthfile)<=4) {
    message("Regenerating relative island width file...")
    pleistodist_relativewidth(points,epsg,intervalfile,mapdir,outdir)
    widthfile <- utils::read.csv(base::paste0(outdir,"/island_relativewidth.csv"))
  }
  if (file.exists(base::paste0(outdir,"/island_area.csv"))==TRUE) {
    areafile <- utils::read.csv(base::paste0(outdir,"/island_area.csv"))
  } else {
    message("No island area file found in output folder, generating island areas from scratch")
    pleistoshape_area(points,epsg,intervalfile,mapdir,outdir)
    areafile <- utils::read.csv(base::paste0(outdir,"/island_area.csv"))
  }

  intervalfile <- utils::read.csv(intervalfile)
  points = sf::st_read(points,fid_column_name="FID")

  #check projection on input points, set to WGS84 (EPSG:4326) if undefined
  if (is.na(sf::st_crs(points)$input)) {
    sf::st_crs(points) = 4326
  }

  #reproject input points to target projection
  points_transformed <- sf::st_transform(points,epsg)
  numintervals = max(intervalfile$Interval)

  #check to see if the points contains a column of unique names, otherwise use FID numbers as identifiers
  if ("Name" %in% colnames(points_transformed) && anyDuplicated(points_transformed$Name) == 0) {
    p1_names <- points_transformed$Name
    p2_names <- points_transformed$Name
  } else {
    p1_names <- points_transformed$FID
    p2_names <- points_transformed$FID
    message("No 'Name' column with unique identifiers detected in points file, defaulting to FID values instead")
  }

  #create dataframe for storing relative migration distance matrix
  relmig <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names,include.equals=T)[,2],
    Island2 = expand.grid(p1_names,p2_names,include.equals=T)[,1],
  )
  for (i in 0:numintervals) {
    intvlname <- base::paste0("interval",i)
    mr <- c()
    message(base::paste0("Calculating relative migration for ",intvlname))
    for (f in 1:nrow(points_transformed)) {
      for (t in 1:nrow(points_transformed)) {
        d12 <- distfile %>% dplyr::filter(Island1==p1_names[f],Island2==p2_names[t]) %>% dplyr::select(dplyr::all_of(intvlname)) %>% as.numeric()
        d21 <- distfile %>% dplyr::filter(Island1==p2_names[t],Island2==p1_names[f]) %>% dplyr::select(dplyr::all_of(intvlname)) %>% as.numeric()
        d = (d12+d21)/2
        if (is.na(d12)==TRUE) {
          mr <- c(mr,NA)
        } else if (d12 == 0) {
          mr <- c(mr,1)
        } else {
          A1 <- areafile %>% dplyr::filter(Island==p1_names[f]) %>% dplyr::select(dplyr::all_of(intvlname)) %>% as.numeric()
          A2 <- areafile %>% dplyr::filter(Island==p2_names[t]) %>% dplyr::select(dplyr::all_of(intvlname)) %>% as.numeric()
          t1 <- widthfile %>% dplyr::filter(Island1==p1_names[f],Island2==p2_names[t]) %>% dplyr::select(dplyr::all_of(intvlname)) %>% as.numeric()
          t2 <- widthfile %>% dplyr::filter(Island1==p2_names[t],Island2==p1_names[f]) %>% dplyr::select(dplyr::all_of(intvlname)) %>% as.numeric()
          migratio <- (atan(t2/(2*d))*A1)/(atan(t1/(2*d))*A2)
          mr <- c(mr,migratio)
        }
      }
    }
    relmig[intvlname] <- mr
  }
  relmig$mean <- apply(relmig[3:ncol(relmig)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  utils::write.csv(relmig,base::paste0(outdir,"/island_netmigration_",disttype,".csv"))
}

#' Estimate island visibility relative to horizon
#'
#' This function estimates the visibility of a destination island relative to a user-defined point located on or above on an origin island for each sea level/time interval.
#' This function works by first calculating the expected horizon distance relative to the origin point, accounting for the local curvature of the Earth's surface
#' as well as the height of both the ground level and the height of the observer above the ground level (users can set this relative height to 0 for non-flying organisms).
#' If the destination island is within the bounds of the horizon radius, this function then calculates the area of the viewshed on the destination island relative to the origin point, to account
#' for the effect of geographical features such as mountains blocking an observer's view. This viewshed area is thus the region visible to an observer standing at a particular
#' height above the origin point. This function also calculates the visible area of the island relative to the total area of the island. If the destination island is not
#' visible, then this function will return a value of 0. Do note that the horizon distance calculation is made relative to the centroid of the destination island, rather than
#' to the specific coordinates of the destination point. Also note that this calculation is a very rough estimate and the calculated visible area may not be identical across different runs (due to the
#' sampling method applied), and does not account for the effect of fog, mist, haze, etc. Lastly, this function calculates the weighted mean visibility over time/sea level of the destination
#' island relative to the origin point.
#' @param points A user-generated multi-point shapefile (.shp) containing at least two points. If the shapefile attribute table contains a column
#' labelled 'Name', the output distance matrix will use the identifiers in this column to label each pairwise comparison. Otherwise, the
#' output distance matrix will use the attribute FID values instead.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param intervalfile This is the master control file generated using either the getintervals_time() or getintervals_sealvl()
#' function that defines the number of intervals, the sea level at each interval, and the duration of each interval. By default, this
#' function will use the "intervals.csv" file stored in the output folder, but users can also specify their own custom interval
#' file (with nice round mean sea level values, for example), although users need to ensure that the same column names are preserved, and
#' be aware that custom interval files may lead to inaccurate weighted mean calculations.
#' @param mapdir The directory containing map outputs from the makemaps() function
#' (i.e. the directory containing the raster_flat, raster_topo, and shapefile folders).
#' @param outdir The output directory for the inter-island visibility matrices. This function generates two files, one estimating
#' the absolute visible area of the destination island, and another estimating the proportion of the destination island visible.
#' @param height The height of the observer relative to the ground level. This can be used to model the visibility of nearby islands to
#' flying animals (e.g. birds, bats, insects). As such, if an origin point is defined on a pixel that is 100m above sea level, a height value of 100 will result in
#' visibility estimates for an organism flying 200m above sea level (with sea level calibrated relative to the sea level of that particular interval).
#' @param plotfigs TRUE/FALSE option to generate island visibility maps for each pairwise island comparison. The function will generate a prompt if the number of
#' pairwise comparisons is very high (>5 points in the points file).
#' @return This function returns two matrices, both in long format, of the visible area of island 2 relative to and origin point on island 1, and the proportion of the
#' visible area on island 2 relative to the entire island area.
#' @examples
#' #create temp directory
#' path <- file.path(tempdir())
#' #create points shapefile
#' points <- sf::st_multipoint(rbind(c(178.68408,-16.82573), c(179.13585,-16.59487)))
#' #convert points to feature geometry
#' points <- sf::st_sfc(points)
#' points <- sf::st_cast(points, "POINT")
#' #set default projection (WGS84)
#' sf::st_crs(points) <- 4326
#' #save shapefile
#' sf::write_sf(points,paste0(path,"/points.shp"))
#' #load bathymetry file
#' fiji <- system.file("extdata","FJ.asc",package="PleistoDist")
#' #generate interval file for 1 interval and 20 kya cutoff time, binning by time
#' getintervals_time(time=20,intervals=1,outdir=path)
#' #generate maps based on interval file, projecting map using EPSG:3141
#' makemaps(inputraster=fiji,epsg=3141,intervalfile=paste0(path,"/intervals.csv"),outdir=path)
#' #calculate inter-island visibility, projecting points using EPSG:3141,
#' #with an observer height of 0 and no figure plotting
#' pleistodist_visibility(points=paste0(path,"/points.shp"),epsg=3141,
#'     intervalfile=paste0(path,"/intervals.csv"),mapdir=path,outdir=path,height=0,plotfigs=FALSE)
#' @export
pleistodist_visibility <- function(points,epsg,intervalfile,mapdir,outdir,height=0,plotfigs=FALSE) {
  intervalfile <- utils::read.csv(intervalfile)
  points <- sf::st_read(points,fid_column_name="FID")

  #check projection on input points, set to WGS84 (EPSG:4326) if undefined
  if (is.na(sf::st_crs(points)$input)) {
    sf::st_crs(points) = 4326
  }

  #extract geodetic latitude for each point and append to points file
  points$geolat <- sf::st_coordinates(sf::st_transform(points,4326))[,2]
  #reproject input points to target projection
  points_transformed <- sf::st_transform(points,epsg)
  numintervals = max(intervalfile$Interval)

  figtest <- TRUE

  if ((nrow(points_transformed) > 5) && (plotfigs == TRUE)) {
    figtest <- utils::askYesNo(
      msg=base::paste0("This function may generate up to a maximum of ",numintervals*as.numeric(nrow(points_transformed))^2," visibility maps. Are you sure you want to plot visibility maps? Yes/No/Cancel"))
  }

  if (is.na(figtest) == TRUE) {
    stop()
  } else if (figtest == FALSE) {
    plotfigs = FALSE
  }

  #check to see if the points contains a column of unique names, otherwise use FID numbers as identifiers
  if ("Name" %in% colnames(points_transformed) && anyDuplicated(points_transformed$Name) == 0) {
    p1_names <- points_transformed$Name
    p2_names <- points_transformed$Name

  } else {
    p1_names <- points_transformed$FID
    p2_names <- points_transformed$FID
    message("No 'Name' column with unique identifiers detected in points file, defaulting to FID values instead")
  }
  #create dataframe for storing visible island area
  visiblearea <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names,include.equals=T)[,2],
    Island2 = expand.grid(p1_names,p2_names,include.equals=T)[,1]
  )
  #create dataframe for storing visible island proportion
  visibleprop <- tibble::tibble(
    Island1 = expand.grid(p1_names,p2_names,include.equals=T)[,2],
    Island2 = expand.grid(p1_names,p2_names,include.equals=T)[,1]
  )
  #define basic earth shape parameters
  a = 6378137.0 #equatorial radius (in metres)
  b = 6356752.3 #polar radius (in metres)
  e2 = 1 - (b/a)^2 #squared eccentricity
  #start looping through each interval
  for (i in 0:numintervals) {
    invector <- sf::st_read(base::paste0(mapdir,"/shapefile/interval",i,".shp")) #load interval shapefile
    #load topo raster and correct for current sea level
    inraster <- terra::rast(base::paste0(mapdir,"/raster_topo/interval",i,".tif"))-intervalfile$MeanDepth[i+1]
    intvlname <- base::paste0("interval",i)
    #calculate resolution of inraster
    reso <- terra::res(inraster)
    #create empty vector for storing visible area and proportion
    visible <- c()
    prop <- c()
    #start 2nd order for-loop to parse source islands
    for (f in 1:nrow(points_transformed)) {
      #calculate parametric latitude
      theta <- as.numeric(atan((b/a)*tan(points$geolat[f]*(pi/180))))
      #start 3rd order for-loop to parse sink islands
      for (t in 1:nrow(points_transformed)) {
        #extract coordinates for island 1
        p1 <- points_transformed$geometry[f]
        #extract coordinates for island 2
        p2 <- points_transformed$geometry[t]
        message(base::paste0("Calculating visibility of ",p2_names[t]," from ",p1_names[f]," for ",intvlname))
        if (f==t) { #write NA for self-comparisons
          visible <- c(visible,NA)
          prop <- c(prop,NA)
        } else if ((nrow(invector[p1, ,op=sf::st_contains]) || nrow(invector[p2, ,op=sf::st_contains])) == 0) { #check to see if island1 or island2 is submerged under water
          visible <- c(visible,NA)
          prop <- c(prop,NA)
          message("One or both islands are underwater during this interval, writing a value of NA")
        } else {
          islandpair <- invector[sf::st_sfc(rbind(p1,p2),crs=epsg), ,op=sf::st_contains]
          if (nrow(islandpair) == 1) {
            visible <- c(visible,NA)
            prop <- c(prop,NA)
            message("Both points are on the same island, writing a value of NA")
          } else {
            island2 <- invector[p2, ,op=sf::st_contains]
            island2_centroid <- sf::st_centroid(island2)
            #create an object with source point from island 1 and centroid point from island 2
            combinedpoints <- sf::st_sfc(rbind(p1,island2_centroid$geometry[1]),crs=epsg)
            #transform points to WGS84 to facilitate calculation of azimuth
            combinedpoints_wgs84 <- sf::st_transform(combinedpoints,4326)
            #calculate angle relative to meridional plane (azimuth)
            phi <- as.numeric(lwgeom::st_geod_azimuth(combinedpoints_wgs84))
            #calculate meridional arc radius (earth's radius for that particular inter-island great circle arc), accounting for sea level change
            Re <- intervalfile$MeanDepth[i+1]+(a*(1-e2*cos(phi)*cos(phi)*cos(theta)*cos(theta))^(3/2))/((1-e2*(1-sin(phi)*sin(phi)*cos(theta)*cos(theta)))^(1/2))
            #calculate observer height (ground height + additional user-defined height)
            htotal <- as.numeric(terra::extract(inraster,sf::st_coordinates(p1)))+height
            #calculate horizon distance
            horizond <- Re*acos(1/(1+(htotal/Re)))
            message(base::paste0("The horizon distance as viewed from ",p1_names[f]," is ",horizond,"m. "))
            #check to see if horizon intersects island 2 and write 0 value if there is no intersection
            if (sf::st_intersects(island2,sf::st_buffer(p1,horizond),sparse=FALSE)==FALSE) {
              visible <- c(visible,0)
              prop <- c(prop,0)
              message(base::paste0(p2_names[t]," is not visible from ",p1_names[f],". Visible area = 0."))
            } else {
              #isolate fragments of island 2 that are within horizon distance
              island2_clipped <- sf::st_intersection(island2,sf::st_buffer(p1,horizond))
              #number of points to sample for line-of-sight calculation, based on inraster resolution
              nsample <- round(as.numeric((sf::st_area(island2_clipped))/(reso[1]*reso[2])))
              #generate sampling points on clipped fragment of island 2
              testpoints <- sf::st_sample(island2_clipped,nsample,type="regular")
              #check to see which points are visible from p1 (code written by Barry Rowlingson (spacedman))
              visibilitymatrix <- viewTo(r=inraster,xy=sf::st_coordinates(p1),xy2=sf::st_coordinates(testpoints),h1=height)
              #subset testpoints visible from p1
              testpoints_visible <- testpoints[which(visibilitymatrix==TRUE)]
              if (length(testpoints_visible)==0) {
                visible <- c(visible,0)
                prop <- c(prop,0)
                message(base::paste0(p2_names[t]," is not visible from ",p1_names[f],". Visible area = 0."))
              } else {
                #finally, estimate visible area, based on the resolution of input raster
                area <- as.numeric(reso[1]*reso[2]*length(testpoints_visible))
                proportion <- as.numeric(area/sf::st_area(island2))
                message(base::paste0(area,"m^2 of ",p2_names[t],", (",signif(proportion,3)," of the total island) is visible from ",p1_names[f]))
                visible <- c(visible,area)
                prop <- c(prop,proportion)
                #plot figures
                if (plotfigs == TRUE) {
                  islandmask <- terra::vect(islandpair)
                  islandcrop <- terra::crop(inraster,islandmask)
                  mapfig <- ggplot2::ggplot()+
                    ggspatial::layer_spatial(raster::raster(islandcrop))+
                    ggplot2::scale_fill_viridis_c(na.value = "transparent",name="Elevation")+
                    ggplot2::geom_sf(data=p1,size=4,shape=17,ggplot2::aes(colour="Observer"))+
                    ggplot2::geom_sf(data=testpoints_visible,size=0.7,ggplot2::aes(colour="Visible Area"))+
                    ggplot2::scale_colour_manual(name=NULL,values = c("Observer"="darkorange2","Visible Area"="darkorange2"),guide=ggplot2::guide_legend(override.aes = list(shape=c(17,16))))+
                    ggplot2::ggtitle(base::paste0("Visibility of ",p2_names[t]," from ",p1_names[f]))+
                    ggplot2::theme(legend.key = ggplot2::element_rect(fill = "darkgrey"))
                  ggplot2::ggsave(filename = base::paste0(outdir,"/visibilitymap_",p1_names[f],"_",p2_names[t],"_interval",i,".png"),device="png",plot=mapfig)
                }
              }
            }
          }
        }
      }
    }
    visiblearea[intvlname] <- visible
    visibleprop[intvlname] <- prop
  }
  visiblearea$mean <- apply(visiblearea[3:ncol(visiblearea)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  utils::write.csv(visiblearea,base::paste0(outdir,"/island_visiblearea.csv"))
  visibleprop$mean <- apply(visibleprop[3:ncol(visibleprop)],1,stats::weighted.mean,intervalfile$TimeInterval,na.rm=TRUE)
  utils::write.csv(visibleprop,base::paste0(outdir,"/island_visibleproportion.csv"))
}
g33k5p34k/PleistoDistR documentation built on Oct. 9, 2022, 5:27 a.m.
rdrr.io home R language documentation Run R code online
CRAN packages Bioconductor packages R-Forge packages GitHub packages
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
g33k5p34k/PleistoDistR
Calculating Island Geomorphologies Over Pleistocene Time

R/calcmetrics.R
In g33k5p34k/PleistoDistR: Calculating Island Geomorphologies Over Pleistocene Time

Defines functions pleistodist_visibility pleistodist_netmig pleistodist_meanshore pleistodist_leastcost pleistodist_leastshore pleistodist_centroid pleistodist_euclidean

Documented in pleistodist_centroid pleistodist_euclidean pleistodist_leastcost pleistodist_leastshore pleistodist_meanshore pleistodist_netmig pleistodist_visibility

R Package Documentation

Browse R Packages

We want your feedback!

g33k5p34k/PleistoDistR Calculating Island Geomorphologies Over Pleistocene Time

R/calcmetrics.R In g33k5p34k/PleistoDistR: Calculating Island Geomorphologies Over Pleistocene Time

Defines functions pleistodist_visibility pleistodist_netmig pleistodist_meanshore pleistodist_leastcost pleistodist_leastshore pleistodist_centroid pleistodist_euclidean

Documented in pleistodist_centroid pleistodist_euclidean pleistodist_leastcost pleistodist_leastshore pleistodist_meanshore pleistodist_netmig pleistodist_visibility

R Package Documentation

Browse R Packages

We want your feedback!

g33k5p34k/PleistoDistR
Calculating Island Geomorphologies Over Pleistocene Time

R/calcmetrics.R
In g33k5p34k/PleistoDistR: Calculating Island Geomorphologies Over Pleistocene Time
