背景
我正在开发一个名为starBliss的自定义星图包。虽然我认为我成功地从我的question中找出了我需要的东西,但看起来有一些边缘情况,事情开始破裂。
有关更多背景信息,请点击此处查看GitHub问题
问题
我目前正在开发一个函数,该函数获取一个星座线shapefile,并基于经度和纬度值使用Lambert Azimuthal Equal Area投影将其重新投影,以创建给定位置在任何给定时间的夜空。
在当前starBliss
包中可以看到一个成功的实施:
#devtools::install_github("benyamindsmith/starBliss")
library(ggplot2)
library(starBliss)
p<- plot_starmap(location= "Toronto, ON, Canada",
date="2022-01-17",
style="black",
line1_text="Toronto",
line2_text ="January 17th, 2023",
line3_text="43.6532° N, 79.3832° W")
ggsave('toronto_black.png', plot = p, width = unit(10, 'in'),
height = unit(15, 'in'))
然而,这种方法确实会遇到一些问题,例如:
(图像被裁剪)
plot_starmap(
location= "Caracas, Venezuela",
date = as.Date("1991-03-17"),
style = "black")
上面的代码创建了一些非对角线(用红色圈出)
一个可重现的示例
我组合了一个函数,它将获取星座线数据并将其转换为当前状态。当我使用该函数并将其绘制在ggplot2
中时,问题仍然存在。
library(tidyverse)
library(sf)
library(tidygeocoder)
library(lubridate)
custom_starmap <- function(location,
date){
# Formatting Date properly
date<- as.Date(date)
# Formatted date
dt<- lubridate::ymd(date)
# Get Latitude and Longitude for ProjString
# For Latitude
suppressMessages(
capture.output(
lat <- tibble(singlelineaddress = location) %>%
geocode(address=singlelineaddress,method = 'arcgis') %>% .[["lat"]]
)
)
# Reference date used for calculating longitude
ref_date <- paste0(year(dt),"01","01",sep="-") %>% ymd()
# Resulting longitude
lon <- (-as.numeric(difftime(ref_date,dt, units="days"))/365)*360
# The CRS
projString <- paste0("+proj=laea +x_0=0 +y_0=0 +lon_0=",round(lon,4), " +lat_0=", round(lat,4))
# Data Transformation
flip <- matrix(c(-1, 0, 0, 1), 2, 2)
hemisphere <- st_sfc(st_point(c(lon, lat)), crs = 4326) %>%
st_buffer(dist = 1e7) %>%
st_transform(crs = projString)
# Data source for constellation lines
url1 <- "https://raw.githubusercontent.com/benyamindsmith/starBliss/main/data/constellations.lines.json"
# Reading Data
invisible(
capture.output(
constellation_lines_sf <- invisible(st_read(url1, stringsAsFactors = FALSE)) %>%
st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=360")) %>%
st_transform(crs = projString) %>%
st_intersection(hemisphere) %>%
filter(!is.na(st_is_valid(.))) %>%
mutate(geometry = geometry * flip)
)
)
st_crs(constellation_lines_sf) <- projString
return(constellation_lines_sf)
}
# The data
df<-custom_starmap(location= "Caracas, Venezuela",
date = as.Date("1991-03-17"))
df
> Simple feature collection with 49 features and 2 fields
Geometry type: GEOMETRY
Dimension: XY
Bounding box: xmin: 8967611 ymin: -8898251 xmax: -8714977 ymax: 9004400
CRS: +proj=laea +x_0=0 +y_0=0 +lon_0=73.9726 +lat_0=10.488
First 10 features:
id rank geometry
1 And 1 MULTILINESTRING ((3542468 3...
2 Ant 3 LINESTRING (-6234955 -52010...
3 Aqr 2 MULTILINESTRING ((8967611 -...
4 Ari 1 LINESTRING (3098546 2071855...
5 Aur 1 MULTILINESTRING ((-1307725 ...
6 Cae 3 LINESTRING (557848.4 -59059...
7 Cam 2 MULTILINESTRING ((-24783.5 ...
8 Cnc 2 MULTILINESTRING ((-6264812 ...
9 CMa 1 MULTILINESTRING ((-2356827 ...
10 CMi 2 LINESTRING (-4432439 -32157..
当我绘制这些数据时,可以看到这些线:
df %>% ggplot()+geom_sf()
(为清楚起见,用红色圈出)
我该如何解决这个问题?我正在使用的CRS格式是否有问题?或者我需要裁剪行吗?
1条答案
按热度按时间46scxncf1#
我认为使用
s2
进行此类练习更安全:创建于2023年1月23日,使用reprex v2.0.2