#' Coordinate transform #' #' WGS84 to GCJ02, GCJ02 to BD09 and vice versa. #' WG84 is used by Google maps and hardware, #' GCJ02 is used by Gaode and Tencent, #' BD09 is used by Baidu. #' #' @param wgs_lon WGS84 longitude #' @param wgs_lat WGS84 latitude #' @param gcj_lon GCJ02 longitude, gaode/tencent #' @param gcj_lat GCJ02 latitude, gaode/tencent #' @param bd_lon BD09 longitude, baidu map #' @param bd_lat BD09 latitude, baidu map #' @return Transformed vector of longitude and latitude. #' @name coordinate_transform #' @examples #' test <- #' structure( #' list( #' startLon = c( #' 113.58850000000001, #' 113.54700000000003, #' 113.54050000000001, #' 113.52866660000001, #' 113.53083329999998, #' 113.54083329999997 #' ), #' startLat = c( #' 33.43299999999999, #' 33.466666599999996, #' 33.49250000000001, #' 33.4581666, #' 33.45766660000001, #' 33.49250000000001 #' ), #' endLon = c( #' 113.5468333, #' 113.54050000000001, #' 113.52866660000001, #' 113.53083329999998, #' 113.54083329999997, #' 113.52249999999998 #' ), #' endLat = c( #' 33.465, #' 33.49250000000001, #' 33.4581666, #' 33.457166599999994, #' 33.49233330000001, #' 33.49250000000001 #' ) #' ), #' row.names = c(NA, -6L), #' class = c("tbl_df", "tbl", "data.frame") #' ) #' #' wgs_lon <- test$startLon[[1]] #' wgs_lat <- test$startLat[[1]] #' wgs2gcj(wgs_lon, wgs_lat)# gaode #' wgs2bd(wgs_lon, wgs_lat)#baidu #' #' t(apply(test, 1, function(x) { #' start_gps <- wgs2bd(x[1], x[2]) #' end_gps <- wgs2bd(x[3], x[4]) #' c(start_gps, end_gps) #' })) #' NULL # long semidiameter a <- 6378245.0 # oblateness ee <- 0.00669342162296594323 # Check whether the coordinate in China. out_of_china <- function(lon, lat) { if(lon < 72.004 || lon > 137.8347) { return(TRUE) } else if(lat < 0.8293 || lat > 55.8271) { return(TRUE) } else { FALSE } } transform_lat <- function(lon, lat) { ret <- -100.0 + 2.0 * lon + 3.0 * lat + 0.2 * lat * lat + 0.1 * lon * lat + 0.2 * sqrt(abs(lon)) ret <- ret + (20.0 * sin(6.0 * lon * pi) + 20.0 * sin(2.0 * lon * pi)) * 2.0 / 3.0 ret <- ret + (20.0 * sin(lat * pi) + 40.0 * sin(lat / 3.0 * pi)) * 2.0 / 3.0 ret <- ret + (160.0 * sin(lat / 12.0 * pi) + 320 * sin(lat * pi / 30.0)) * 2.0 / 3.0 ret } transform_lon <- function(lon, lat) { ret <- 300.0 + lon + 2.0 * lat + 0.1 * lon * lon + 0.1 * lon * lat + 0.1 * sqrt(abs(lon)) ret <- ret + (20.0 * sin(6.0 * lon * pi) + 20.0 * sin(2.0 * lon * pi)) * 2.0 / 3.0 ret <- ret + (20.0 * sin(lon * pi) + 40.0 * sin(lon / 3.0 * pi)) * 2.0 / 3.0 ret <- ret + (150.0 * sin(lon / 12.0 * pi) + 300.0 * sin(lon / 30.0 * pi)) * 2.0 / 3.0 ret } #' @export #' @rdname coordinate_transform wgs2gcj <- function(wgs_lon, wgs_lat) { if(out_of_china(wgs_lon, wgs_lat)){ warning("gps out of China, set to 0.",call. = FALSE) return(c(0.0, 0.0)) } dlat <- transform_lat(wgs_lon - 105.0, wgs_lat - 35.0) dlon <- transform_lon(wgs_lon - 105.0, wgs_lat - 35.0) rad_lat <- wgs_lat / 180.0 * pi magic <- 1 - ee * (sin(rad_lat))^2 dlat <- (dlat * 180.0) / ((a * (1-ee)) / (magic * sqrt(magic)) * pi) dlon <- (dlon * 180.0) / (a / sqrt(magic) * cos(rad_lat) * pi) gd_lat <- wgs_lat + dlat gd_lon <- wgs_lon + dlon c(gd_lon, gd_lat) } #' @export #' @rdname coordinate_transform gcj2wgs <- function(gcj_lon, gcj_lat) { if(out_of_china(gcj_lon, gcj_lat)){ warning("gps out of China, set to 0.",call. = FALSE) return(c(0.0, 0.0)) } dlat <- transform_lat(gcj_lon - 105.0, gcj_lat - 35.0) dlon <- transform_lon(gcj_lon - 105.0, gcj_lat - 35.0) rad_lat <- gcj_lat/ 180.0 * pi magic <- 1 - ee * (sin(rad_lat))^2 dlat <- (dlat * 180.0) / ((a * (1- ee)) / (magic * sqrt(magic)) * pi) dlon <- (dlon * 180.0) / (a / sqrt(magic) * cos(rad_lat) * pi) wgs_lat <- gcj_lat * 2 - (gcj_lat + dlat) wgs_lon <- gcj_lon * 2 - (gcj_lon + dlon) c(wgs_lon, wgs_lat) } #' @export #' @rdname coordinate_transform gcj2bd <- function(gcj_lon, gcj_lat) { if(out_of_china(gcj_lon, gcj_lat)){ warning("gps out of China, set to 0.",call. = FALSE) return(c(0.0, 0.0)) } x_pi <- pi * 3000.0 / 180.0 z <- sqrt(gcj_lon^2 + gcj_lat^2) + 0.00002 * sin(gcj_lat * x_pi) theta <- atan2(gcj_lat, gcj_lon) + 0.000003 * cos(gcj_lon * x_pi) bd_lon <- z * cos(theta) + 0.0065 bd_lat <- z * sin(theta) + 0.006 c(bd_lon, bd_lat) } #' @export #' @rdname coordinate_transform bd2gcj <- function(bd_lon, bd_lat) { if(out_of_china(bd_lon, bd_lat)){ warning("gps out of China, set to 0.",call. = FALSE) return(c(0.0, 0.0)) } x_pi <- pi * 3000.0 / 180.0 x <- bd_lon - 0.0065 y <- bd_lat - 0.006 z <- sqrt(x * x + y * y) - 0.00002 * sin(y * x_pi) theta <- atan2(y, x) - 0.000003 * cos(x * x_pi) gd_lon <- z * cos(theta) gd_lat <- z * sin(theta) c(gd_lon, gd_lat) } #' @export #' @rdname coordinate_transform wgs2bd <- function(wgs_lon, wgs_lat){ gcj <- wgs2gcj(wgs_lon,wgs_lat) bd <- gcj2bd(gcj[[1]],gcj[[2]]) bd_lon <- bd[[1]] bd_lat <- bd[[2]] c(bd_lon, bd_lat) } #' @export #' @rdname coordinate_transform bd2wgs <- function(bd_lon, bd_lat){ gcj <- bd2gcj(bd_lon, bd_lat) wgs <- gcj2wgs(gcj[[1]], gcj[[2]]) wgs_lon <- wgs[[1]] wgs_lat <- wgs[[2]] c(wgs_lon, wgs_lat) } #' Generate GPS lines data for echarts #' #' Your can see the package/misc files. #' #' @param df data.frame with starLon, starLat, endLon, endLat. #' #' @return json txt #' @export #' #' @examples #' \dontrun{ #' gps_trans <- function(df) { #' test <- select(df, startLon, startLat, endLon, endLat) #' as_tibble(t(apply(test, 1, function(x) { #' x <- as.numeric(x) #' start_gps <- wgs2bd(x[1], x[2]) #' end_gps <- wgs2bd(x[3], x[4]) #' c(start_gps[1], start_gps[2], end_gps[1], end_gps[2]) #' }))) #' } #' usr_data <- purrr::map(usr_ids, ~ filter(pam_data, userId == .x)) #' all <- map(lapply(usr_data, gps_trans), ~ fmt_gps_json(.x)) #' #'} fmt_gps_json <- function(df) { list_mat <- df %>% purrr::pmap( ~ matrix(c(..1, ..2, ..3, ..4), byrow = TRUE, ncol = 2)) tibble::tibble(coords = list_mat) %>% jsonlite::toJSON() } # to rad to_rad <- function(gps){ gps / 180 * pi } #' Calculate earth distance between two points gps longitude and latitude. #' #' @param lon_s Start gps longitude #' @param lat_s Start gps latitude #' @param lon_e End gps longitude #' @param lat_e End gps latitude #' #' @return Distance of double type in meters. #' @export #' #' @examples #' cal_gps_dist(121.434174,31.158213,121.444366,31.169403) #' # meters cal_gps_dist <- function(lon_s, lat_s, lon_e, lat_e){ lon_s <- to_rad(lon_s) lat_s <- to_rad(lat_s) lon_e <- to_rad(lon_e) lat_e <- to_rad(lat_e) r <- 6378160 # meter tmp <- (sin((lat_e - lat_s) / 2))^2 + cos(lat_s)*cos(lat_e) * (sin((lon_e - lon_s)/2))^2 2 * r * asin(sqrt(tmp)) }