LatLong to UTM, NAD83

R_utm.R

to_utm <- function(x, y){
  E <- 0.00669438
  R <- 6378137
  
  K0 <- 0.9996
  
  E2 <- E * E
  E3 <- E2 * E
  E_P2 <- E / (1.0 - E)
  
  SQRT_E <- (1.0 - E) ^ (1/2)
  E.1 <- (1 -SQRT_E) / (1 + SQRT_E)
  E.2 <- E.1 * E.1
  E.3 <- E.2 * E.1
  E.4 <- E.3 * E.1
  E.5 <- E.4 * E.1
  
  M1 <- (1 - E / 4 - 3 * E2 / 64 - 5 * E3 / 256)
  M2 <- (3 * E / 8 + 3 * E2 / 32 + 45 * E3 / 1024)
  M3 <- (15 * E2 / 256 + 45 * E3 / 1024)
  M4 <- (35 * E3 / 3072)
  
  P2 <- (3. / 2 * E.1 - 27. / 32 * E.3 + 269. / 512 * E.5)
  P3 <- (21. / 16 * E.2 - 55. / 32 * E.4)
  P4 <- (151. / 96 * E.3 - 417. / 128 * E.5)
  P5 <- (1097. / 512 * E.4)
  

   zone_number <- as.integer(as.matrix(((x + 180) / 6) + 1))
   zone_number <- data.frame(zone_number)
  
  lat_rad <- (pi / 180) * y
  lat_sin <- sin(lat_rad)
  lat_cos <- cos(lat_rad)
  lat_tan <- lat_sin / lat_cos
  lat_tan2 <- lat_tan * lat_tan
  lat_tan4 <- lat_tan2 * lat_tan2
  
  lon_rad <- (pi / 180) * x
  central_lon = ((zone_number - 1) * 6) - 180 + 3
  central_lon_rad = (pi / 180) * central_lon
  
  n <- R / (1 - E * lat_sin ^ 2) ^ (1/2)
  c <- E_P2 * lat_cos ^ 2
  
  a <- lat_cos * (lon_rad - central_lon_rad)
  a2 <- a * a
  a3 <- a2 * a
  a4 <- a3 * a
  a5 <- a4 * a
  a6 <- a5 * a
  
  m <- R * (M1 * lat_rad - 
              M2 * sin(2 * lat_rad) +
              M3 * sin(4 * lat_rad) - 
              M4 * sin(6 * lat_rad))
  
  easting <-  K0 * n * (a +
                          a3 / 6 * (1 - lat_tan2 + c) +
                          a5 / 120 * (5 - 18 * lat_tan2 * lat_tan4 + 72 * c - 58 * E_P2)) + 500000
  northing <- K0 * (m + n * lat_tan * (a2 / 2 +
                                         a4 / 24 * (5 - lat_tan2 + 9 * c + 4 * c**2) +
                                         a6 / 720 * (61 - 58 * lat_tan2 + lat_tan4 + 600 * c - 330 * E_P2)))

  output <- c(easting, northing, zone_number)
  return(output)
}