I’m working on a large dataset of woodland plots with in the SEOSAW network . As part of this work I want to compare the spatial distribution of tree stems within the plots. Some of the plots have their stem locations recorded as decimal latitude/longitude co-ordinates, presumably done with a handheld GPS, while other plots have their stems recorded as metre coordinates from the plot edges, presumably these plots were done with a set of tape measures.
I can’t easily compare plots with these different methods and it’s good to be consistent, so I set about writing a function in R that can take stem locations and plot corner loctions in latitude/longitude coordinates and convert them to rought X Y metre coordinates. The function can be found below, and also here . This is what I came up with:
# Function allowing lat long to x y coordinate conversion
# John Godlee (johngodlee@gmail.com)
# 2018_04_20
# Packages ----
library(dplyr) # Data manip.
library(rgdal) # CRS stuff
library(raster) # spLines()
library(rgeos) # gDistance()
# Function ----
# stem_id = a unique ID string for each stem
# stem_lon, stem_lat = latitude longitude coordinates for each stem
# corner_id = unique ID string for each plot corner
# corner_lon, corner_lat = latitude longitude coordinates for each plot corner
latlong_xy <- function(stem_id, stem_lon, stem_lat, corner_id, corner_lon, corner_lat){
long_2_utm <- function(x,y) {
paste("UTM zone ",
(floor((x + 180)/6) %% 60) + 1,
ifelse(y < 0, "S", "N"),
sep = "")
}
epsg <- make_EPSG() # Create list of EPSG dataset to search for CRS
wgs84 <- epsg[grep("WGS 84", epsg$note, ignore.case=TRUE),] # Search for wgs84
wgs84[grep("longlat", wgs84$prj4, ignore.case=TRUE),] # grep proj4string to check
wgs84_crs <- CRS(wgs84[grep("longlat", wgs84$prj4, ignore.case=TRUE),]$prj4[2]) # Store string as vector
# Must change the UTM zone to match location of your plot.
utm_id <- long_2_utm(mean(corner_lon),
mean(corner_lat))
utm_zone_crs <- CRS(wgs84[grep(utm_id, wgs84$note, ignore.case=TRUE),]$prj4[1]) # grep for UTM zone and store
# Convert stem data to utm
stems_points <- SpatialPointsDataFrame(as.matrix(cbind(stem_lon, stem_lat)), # extract only long lat coords
proj4string = wgs84_crs,
data = data.frame(stem_id))
# Transform SPDF to utm
stems_points_sp_utm <- spTransform(stems_points, utm_zone_crs)
# Convert back to dataframe
stems_points_df_utm <- as.data.frame(stems_points_sp_utm)
# Give column names
colnames(stems_points_df_utm) <- c("id", "x_utm", "y_utm")
# Convert plot corners to utm ----
plot_corners_clean <- data.frame(corner_id, "x" = corner_lon, "y" = corner_lat)
plot_corners_points <- SpatialPointsDataFrame(plot_corners_clean[,2:3],
proj4string = wgs84_crs,
data = data.frame(plot_corners_clean[,1]))
# Transform SPDF to utm
plot_corners_points_utm <- spTransform(plot_corners_points, utm_zone_crs)
# Convert back to dataframe
plot_corners_df_utm <- as.data.frame(plot_corners_points_utm)
# Give column names
colnames(plot_corners_df_utm) <- c("id", "x_coord", "y_coord")
# Make spatial lines from corners ----
# Get corner locations
corner_nw <- plot_corners_df_utm[1,]
corner_ne <- plot_corners_df_utm[2,]
corner_sw <- plot_corners_df_utm[4,]
# Create x axis line
x_line <- rbind(corner_nw, corner_ne) %>%
dplyr::select(x_coord, y_coord) %>%
as.matrix(.) %>%
spLines(., crs = utm_zone_crs)
# Create y axis line
y_line <- rbind(corner_nw, corner_sw) %>%
dplyr::select(x_coord, y_coord) %>%
as.matrix(.) %>%
spLines(., crs = utm_zone_crs)
# Calculate x y distances from line to point for each stem and append to data frame ----
stems_points_df_utm$x_coord <- as.vector(gDistance(stems_points_sp_utm, y_line, byid = T))
stems_points_df_utm$y_coord <- as.vector(gDistance(stems_points_sp_utm, x_line, byid = T))
# Clean up old and unnecessary columns
stems_loc_df <- stems_points_df_utm %>%
dplyr::select("id", "x_coord", "y_coord")
stems_loc_df
}