Skip to content

Instantly share code, notes, and snippets.

@mutolisp
Created October 8, 2020 09:03
Show Gist options
  • Save mutolisp/cd49355323330d91cb7cfb3c91058ecc to your computer and use it in GitHub Desktop.
Save mutolisp/cd49355323330d91cb7cfb3c91058ecc to your computer and use it in GitHub Desktop.
twTransfrom
# install sp package before use
# read coordinates csv
twTransform <- function(input = file.choose(), x='x', y='y' ,from = 'wgs84', to= 'twd97'){
require(sp)
coords <- read.csv(input, sep = ',')
coords <- cbind(coords[x], coords[y])
colnames(coords) <- c('x', 'y')
pts <- coordinates(coords)
# projection
wgs84 <- CRS('+proj=longlat +ellps=WGS84')
twd97 <- CRS('+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ')
twd67 <- CRS('+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +units=m +no_defs ')
# Set projection
if ( from == 'wgs84'){
pts.prj <- SpatialPoints(pts, wgs84)
} else if ( from == 'twd97' ){
pts.prj <- SpatialPoints(pts, twd97)
} else if ( from == 'twd67' ){
pts.prj <- SpatialPoints(pts, twd67)
} else {
print('Please use wgs84, twd97 or twd67')
}
# transform to
if ( to == 'wgs84'){
pts.reproj <- spTransform(pts.prj, wgs84)
} else if ( from == 'twd97' ){
pts.reproj <- spTransform(pts.prj, twd97)
} else if ( from == 'twd67' ){
pts.reproj <- spTransform(pts.prj, twd67)
} else {
print('Please use wgs84, twd97 or twd67')
}
pts.reproj.coords <- as.data.table(pts.reproj)
return(pts.reproj.coords)
}
# call function
twTransform('/tmp/observations-24196.csv', x = 'longitude', y = 'latitude')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment