convert Topcon total station data to known coordinates (e.g., lat/long; northing/easting)

TopCon_conversion.R

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

##### 
# User inputs
#
# data should be the raw TopConLink output of the points tab from the "measurement" file
tsData <- read.delim("C:/RDATA/surveying/data_150526SWD.txt", header=TRUE, sep = "\t", na.strings=".")

# the TSObs dataset is needed to convert between total station orientations
tsobs <- read.delim("C:/RDATA/surveying/data_150526SWD_tsobs.txt", header=TRUE, sep = "\t", na.strings=".")

# the "coordinate" dataset is really just for the total station coordinates
tsCoord <- read.delim("C:/RDATA/surveying/data_150526SWD_coord.txt", header=TRUE, sep = "\t", na.strings=".")

# create a vector to capture replicated points
# in my data these labels all contain 'SET' or 'BM'
repList <- paste0(c("SET", "BM"), 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(
  Code = as.character(c("SET13", "SET14", "SET15")),
  y = c(183992.8752847, 183989.2218590, 183997.5237405),
  x = c(256468.7365637, 256480.5637925, 256501.1030425),
  z = c(1.104, 1.190, 1.210),
  y2 = c(41.1163561, 41.1163239, 41.1163999),
  x2 = c(-73.3254669, -73.3253258, -73.3250819)
)

# or, load as a datafile: 
# knowns <- read.delim("C:/RDATA/surveying/data_RTKswdSETs.txt", header=TRUE, sep = "\t", na.strings=".")
# names(knowns) <- c("Code", "y", "x", "z", "y2", "x2")
# knowns$Code <- as.character(knowns$Code)
#####

#####
# custom functions
# This function detects the total station orientation and adds a column labeling the side
getSide <- function(TSobs, ZenithAngle = "Zenith.Angle") {
  # TSObs = "TSObs" dataset. A column will be added to this, and 
  #  ZenithAngle must be a column in this dataset
  TSobs$side <- as.character("NA")
  startData <- sapply(strsplit(as.character(TSobs[, ZenithAngle]), "'"), head, 1)
  midData <- substr(startData, 1, nchar(startData) - 3) # get degrees: 180/360 is straight up/down
  
  for (i in 1:nrow(TSobs)) {
  if (as.numeric(midData[i]) >= 180) {
    TSobs$side[i] <- "a"
  } else { # 
    TSobs$side[i] <- "b"
  }
  }
  TSobs
}

# this function generates a transform matrix:
# If "unknownData" object doesn't have columns 
# named "Code", "Ground.Northing..m." and "Ground.Easting..m.",
# insert those column names as 'y.unk' and 'x.unk'
genTransform <- function(knownData = knowns, y.col = "y", x.col = "x", 
                         unknownData = SETpoints, y.unk = "Ground.Northing..m.", x.unk = "Ground.Easting..m.", 
                         n = 3) {
  
  # 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 careful that transform matrix and NEunknowns have the same order of 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
}
#####


# make sure eastings and northings are numeric
tsData$Ground.Northing..m. <- as.numeric(as.character(tsData$Ground.Northing..m.))
tsData$Ground.Easting..m. <- as.numeric(as.character(tsData$Ground.Easting..m.))
tsData$Elevation..m. <- as.numeric(as.character(tsData$Elevation..m.))

# label total station orientations
tsobs <- getSide(tsobs)

# get total station coordinate range
tsEasting <- as.numeric(as.character(tsCoord$Ground.Easting[1]))
tsNorthing <- as.numeric(as.character(tsCoord$Ground.Northing[1]))

# 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])

# join the total station datasets
tsData <- cbind(tsData[-1, ], tsobs) # sometimes this is 'HAngle.Residual', other times it's 'HAngle Residual'
tsData$Code <- as.character(tsData$Code)

# find side with the most datapoints
# this requires a user-entered "side" column in the raw data 
# indicating the total station orientation - a or b
if (sum(tsData$side %in% "a") > sum(tsData$side %in% "b")) {
  dominantSide <- "a"
} else {
  dominantSide <- "b"
}
  


# the short way to convert between total station orientations:
for (i in 1:nrow(tsData)) {
  if(!tsData$side[i] %in% dominantSide) { 
    HD.h <- tsData$Ground.Northing..m.[i] - tsNorthing
    VD.h <- tsData$Ground.Easting..m.[i] - tsEasting

    tsData$Ground.Northing..m.[i] <- tsNorthing - HD.h
    tsData$Ground.Easting..m.[i] <- tsEasting - VD.h
  }
}


# find average values for replicated measurements
ctrlPoints <- ddply(tsData[grep(repList, tsData$Code), ], .(Code), summarise,
                    Ground.Northing..m. = mean(Ground.Northing..m., na.rm = T),
                    Ground.Easting..m. = mean(Ground.Easting..m., na.rm = T),
                    Elevation..m. = mean(Elevation..m., 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 <- grep(paste(names(SETpoints), collapse = "|"), names(tsData))[1:ncol(SETpoints)]
tsData2 <- rbind(ctrlPoints, tsData2[, seeker.cols])



# try with single observation
# marker <- "SET13"
# test.pts <- ctrlPoints[ctrlPoints$Code %in% marker, c("Ground.Easting..m.", "Ground.Northing..m.")]
# coordConv(inputData = test.pts, transformMatrix = genTransform(), xyPoints = test.pts,
#           colNames = c("Easting.SP", "Northing.SP"))
# # try with dataframe
# test <- coordConv(inputData = SETpoints, transformMatrix = genTransform(), xyPoints = SETpoints[, c(3, 2)], # column order is v. important
#           colNames = c("Easting.SP", "Northing.SP"))
# check result:
# test
# knownData

# 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("Ground.Easting..m.", "Ground.Northing..m.")] # column order is v. important
newPoints <- coordConv(inputData = tsData2, transformMatrix = NETransMat, 
                      xyPoints = rawInput, 
                      colNames = c("Easting.SP", "Northing.SP"))
newPoints <- coordConv(inputData = newPoints, transformMatrix = LatTransMat, 
                       xyPoints = rawInput, 
                       colNames = c("Long", "Lat"))


# export data
write.csv(newPoints, file = "exportedData.csv")