Skip to content

Instantly share code, notes, and snippets.

@ernestguevarra
Last active December 14, 2018 20:24
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ernestguevarra/e8fcb497729cb435bae1c1ef41a4f306 to your computer and use it in GitHub Desktop.
Save ernestguevarra/e8fcb497729cb435bae1c1ef41a4f306 to your computer and use it in GitHub Desktop.
library(dplyr)
library(sp)
library(rgdal)
library(tibble)
find_UTM_zone <- function(longitude, latitude) {
# Special zones for Svalbard and Norway
if (latitude >= 72.0 && latitude < 84.0 )
if (longitude >= 0.0 && longitude < 9.0)
return(31);
if (longitude >= 9.0 && longitude < 21.0)
return(33)
if (longitude >= 21.0 && longitude < 33.0)
return(35)
if (longitude >= 33.0 && longitude < 42.0)
return(37)
(floor((longitude + 180) / 6) %% 60) + 1
}
find_UTM_hemisphere <- function(latitude) {
ifelse(latitude > 0, "north", "south")
}
# returns a DF containing the UTM values, the zone and the hemisphere
longlat_to_UTM <- function(long, lat, units = 'm') {
df <- data.frame(
id = seq_along(long),
x = long,
y = lat
)
sp::coordinates(df) <- c("x", "y")
hemisphere <- find_UTM_hemisphere(lat)
zone <- find_UTM_zone(long, lat)
sp::proj4string(df) <- sp::CRS("+init=epsg:4326")
CRSstring <- paste0(
"+proj=utm +zone=", zone,
" +ellps=WGS84",
" +", hemisphere,
" +units=", units)
if (dplyr::n_distinct(CRSstring) > 1L)
stop("multiple zone/hemisphere detected")
res <- sp::spTransform(df, sp::CRS(CRSstring[1L])) %>%
tibble::as_data_frame() %>%
dplyr::mutate(
zone = zone,
hemisphere = hemisphere
)
res
}
UTM_to_longlat <- function(utm_df, zone, hemisphere) {
CRSstring <- paste0("+proj=utm +zone=", zone, " +", hemisphere)
utmcoor <- sp::SpatialPoints(utm_df, proj4string = sp::CRS(CRSstring))
longlatcoor <- sp::spTransform(utmcoor, sp::CRS("+init=epsg:4326"))
tibble::as_data_frame(longlatcoor)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment