Script converts total station data (Northwest Instruments NTS02) to lat/long or state plane (or other) coordinate systems using three control points.

NorthwestInstruments_conversion.R

# Script converts total station data (Northwest Instruments NTS02) to 
# lat/long or state plane (or other) coordinate systems using three control points.

##### 
# User inputs
#
# set input data: raw output from the NorthWest Instruments software
tsData <- read.csv("C:/RDATA/FelixNeck/data/data_elevation_20170504.csv", header=TRUE, na.strings=".")

# set the location of output file containing northing/easting and lat/long data 
exportedData <- "C:/RDATA/ESRIDATA/FelixNeck_data/exportedData.csv"

# create a vector of point labels for replicated points (they'll be averaged)
# for a dataset where replicated points all contain 'SET' or 'BM':
repList <- paste0(c("SET", "BM", "WLL"), collapse = "|")

# enter your control points - insert names as they appear in the dataset, and their 
# known coordinates (i.e., lat/long, state plane) here:
knowns <- data.frame( # note: these coordinates are very approximate (handheld GPS)
  code = as.character(c("BM1", "BM2", "BM3", "BM4", "BM5", "WLL_post")),
  y = c(4586547.788, 4586514.305, 4586427.574, 4586690.426, 4586732.158, 4586614.593),
  x = c(369867.517, 369792.297, 369755.3878, 369759.5583, 369733.7123, 369869.449),
  z = c(1.36300, 1.02750, 0.59650, 0.917, 0.62550, -0.5535), # elevations rel. to MHW 1983-2001
  y2 = c(NA, NA, 41.4188623, 41.42122976, 41.42160133, 41.42056474),
  x2 = c(NA, NA, -70.55854339, -70.5585501, -70.55886827, 	-70.55721922)
)

# or, load data: 
# knowns <- read.csv("C:/RDATA/surveying/data_RTKswdSETs.csv", header=TRUE, na.strings=".")
# names(knowns) <- c("code", "y", "x", "z", "y2", "x2")
#####



# install.packages("plyr")
library(plyr)



#####
# this function generates a transform matrix:
genTransform <- function(knownData = knowns, y.col = "y", x.col = "x", # columns with known x and y data 
                         unknownData = SETpoints, y.unk = "N", x.unk = "E", 
                         n = 3) {
  # If "unknownData" object doesn't have columns 
  # named "code", "N" and "E",
  # insert those column names for 'y.unk' and 'x.unk' arguments
  # 'n' sets the number of control points to use. this will be developed in future.
  # if user provides more than three points with known coordinates, 
  # use the first three found in both datasets
  
  tiePoints <- c(as.character(knownData$code[knownData$code %in% unknownData$code])[1:n])
  
  # re-format to get the datasets in the same order
  data <- knownData[knownData$code %in% tiePoints, ]
  data$fill.col <- 1
  data$x.un <- data$y.un <- NA
  for (i in 1:nrow(data)) {
    data$y.un[i] <- unknownData[unknownData$code %in% tiePoints[i], y.unk]
    data$x.un[i] <- unknownData[unknownData$code %in% tiePoints[i], x.unk]
  }
  print(noquote(c("Used the following control points: ", paste0(data$code[1:n], collapse = ", "))))
  # get known coordinates
  k <- as.matrix(data[c(x.col, y.col, "fill.col")])
  u <- as.matrix(data[c("x.un", "y.un", "fill.col")])
  
  A <- t(solve(u, k)) # this matches matlab output 'A'
}


# function to convert between coordinate systems
# be sure that transform matrix and NEunknowns have the same column order for northings/eastings
# inputData: dataset to append to
# transformMatrix: matrix used for conversion, generated by genTransform()
# xyPoints: unknown x/y points - it's very important that these appear in the order x, y
# colNames: a vector of strings containing names of columns to append
coordConv <- function (inputData, transformMatrix, xyPoints, colNames) {
  # add columns to receive data
  inputData[, colNames[2]] <- inputData[, colNames[1]] <- NA
  
  xyPoints$fill.col <- 1
  if (!nrow(xyPoints) == 1) {
    for (i in 1:nrow(xyPoints)) {
      outputTemp <- transformMatrix %*% as.numeric(xyPoints[i, ])
      inputData[i, colNames[1]] <- outputTemp[1]
      inputData[i, colNames[2]] <- outputTemp[2]
    }    
  } else { 
    outputTemp <- transformMatrix %*% as.numeric(xyPoints)
    inputData[colNames[1]] <- outputTemp[1]
    inputData[colNames[2]] <- outputTemp[2]
  }
  
  inputData
}
#####

if ("Code" %in% names(tsData)) {
  names(tsData) <- gsub("Code", "code", names(tsData))
  }
  
# make sure eastings and northings are numeric
tsData$N <- as.numeric(as.character(tsData$N))
tsData$E <- as.numeric(as.character(tsData$E))
tsData$Z <- as.numeric(as.character(tsData$Z))
tsData$code <- as.character(tsData$code)


# set total station coordinates. for NWI this always seems to be zero
tsEasting  <- 0
tsNorthing <- 0

# if user provides more than three points with known coordinates, 
# build a transform matrix based on the the first three found in both datasets
pointsToUse <- c(as.character(knowns$code[knowns$code %in% tsData$code])[1:3])


# find average values for replicated measurements
ctrlPoints <- ddply(tsData[grep(repList, tsData$code), ], .(code), summarise,
                    N = mean(N, na.rm = T),
                    E = mean(E, na.rm = T),
                    Z = mean(Z, na.rm = T)
)

# now, isolate just the control points you'll be using for the transform matrix
SETpoints <- ctrlPoints[grep(paste0(pointsToUse, collapse = "|"), ctrlPoints$code), ]

# remove replicated points from the dataset
tsData2 <- tsData[-c(grep(c(repList), tsData$code)), ]

# columns sought from tsData
seeker.cols <- which(names(tsData) %in% names(SETpoints))
tsData2 <- rbind(ctrlPoints, tsData2[, seeker.cols])





# convert total station coordinates to northing/eastings and lat/long
NETransMat <- genTransform(unknownData = ctrlPoints[grep(paste0(pointsToUse, collapse = "|"), ctrlPoints$code), ])
LatTransMat <- genTransform(unknownData = ctrlPoints[grep(paste0(pointsToUse, collapse = "|"), ctrlPoints$code), ], 
                            x.col = "x2", y.col = "y2") # I put the long/lat coordinates in columns labeled x2 and y2

rawInput <- tsData2[, c("E", "N")] # column order is v. important
newPoints <- coordConv(inputData = tsData2, transformMatrix = NETransMat, 
                       xyPoints = rawInput, 
                       colNames = c("Easting.SP", "Northing.SP"))
newPoints2 <- coordConv(inputData = newPoints, transformMatrix = LatTransMat, 
                       xyPoints = rawInput, 
                       colNames = c("Long", "Lat"))
newPoints3 <- cbind(newPoints, newPoints2[, c("Long", "Lat")])

# export data
write.csv(newPoints3, file = exportedData, row.names = FALSE)