R/internal_clean_coordinate.R

In azizka/CoordinateCleaner: Automated Cleaning of Occurrence Records from Biological Collections

# Check projection of custom reference and reproject to wgs84 if necessary
reproj <- function(ref) {
  wgs84 <- "+proj=longlat +datum=WGS84 +no_defs"
  ref_crs <- terra::crs(ref, proj = TRUE)
  # if no projection information is given assume wgs84
  if (ref_crs == "") {
    warning(
      "no projection information for reference found,
            assuming '+proj=longlat +datum=WGS84 +no_defs'"
    )
    terra::crs(ref) <- wgs84
  } else {
    if (is.na(ref_crs)) {
      warning(
        "no projection information for reference found,
            assuming '+proj=longlat +datum=WGS84 +no_defs'"
      )
      ref <- terra::project(ref, wgs84)
    } else if (ref_crs != wgs84) {
      #otherwise reproject
      ref <- terra::project(ref, wgs84)
      warning("reprojecting reference to '+proj=longlat +datum=WGS84 +no_defs'")
    }
  }
  return(ref)
}

# A function to create a raster from an input dataset, used in cc_outl for spatial thinning, 
# based on a point dataset and a raster resolution

ras_create <- function(x, lat, lon,  thinning_res){
  # get data extend
  ex <- terra::ext(terra::vect(x[, c(lon, lat)], 
                               geom = c(lon, lat))) + thinning_res * 2
  
  # check for boundary conditions
  if (ex[1] < -180 | ex[2] > 180 | ex[3] < -90 | ex[4] > 90) {
    warning("fixing raster boundaries, assuming lat/lon projection")
    
    if (ex[1] < -180) {
      ex[1] <- -180
    }
    
    if (ex[2] > 180) {
      ex[2] <- 180
    }
    
    if (ex[3] < -90) {
      ex[3] <- -90
    }
    
    if (ex[4] > 90) {
      ex[4] <- 90
    }
  }
  
  # create raster
  ras <- terra::rast(x = ex, resolution = thinning_res)
  
  # set cell ids
  vals <- seq_len(terra::ncell(ras))
  ras <- terra::setValues(ras, vals)
  
  return(ras)
}


# A function to get the distance between raster midpoints and 
#output a data.frame with the distances and the cell IDs as row and column names for cc_outl

ras_dist <-  function(x, lat, lon, ras, weights) {
  #x = a data.frame of point coordinates, ras = a raster with cell IDs as layer,
  #weight = logical, shall the distance matrix be weightened by the number of
  #points per cell? assign each point to a raster cell
  pts <- terra::extract(x = ras, 
                        y = terra::vect(x[, c(lon, lat)],
                                        geom = c(lon, lat),
                                        crs = ras))
  
  # convert to data.frame
  midp <- data.frame(terra::as.points(ras), 
                     terra::xyFromCell(ras, 1:terra::ncell(ras)))
  
  # retain only cells that contain points
  midp <- midp[midp$lyr.1 %in% unique(pts$lyr.1), , drop = FALSE]
  
  # order
  midp <- midp[match(unique(pts$lyr.1), midp$lyr.1), , drop = FALSE]
  
  
  # calculate geospheric distance between raster cells with points
  dist <- geosphere::distm(midp[, c("x", "y")], 
                           fun = geosphere::distHaversine) / 1000
  
  # set rownames and colnames to cell IDs
  dist <- as.data.frame(dist, row.names = as.integer(midp$lyr.1))
  names(dist) <- midp$lyr.1
  
  if (weights) {
    # approximate within cell distance as half 
    # the cell size, assumin 1 deg = 100km
    # this is crude, but doesn't really matter
    dist[dist == 0] <- 100 * mean(terra::res(ras)) / 2
    
    # weight matrix to account for the number of points per cell
    ## the number of points in each cell
    cou <- table(pts$lyr.1)
    
    ## order
    cou <- cou[match(unique(pts$lyr.1), names(cou))]
    
    # weight matrix, representing the number of distances between or within the cellse (points cell 1 * points cell 2)
    wm <- outer(cou, cou)
    
    # multiply matrix elements to get weightend sum
    dist <- round(dist * wm, 0)
    
    dist <- list(pts = pts, dist = dist, wm = wm)
  } else {
    # set diagonale to NA, so it does not influence the mean
    dist[dist == 0] <- NA
    
    dist <- list(pts = pts, dist = dist)
  }
  
  return(dist)
}