Context
I am in the process of developing a custom star map package called
A reproducible example
I put together a function that will get the constellation lines data and transform it as it presently is. When I use that and plot in ggplot2
it with geom_sf
the problem still exists.
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..
When I plot this data the lines in question can be seen:
df %>% ggplot() geom_sf()
(circled in red for clarity)
How do I fix this? Is there an issue with the format of the CRS that I am using? Or do I need to crop the lines?
CodePudding user response:
I think it is safer to use s2
for this type of exercise:
library(tidyverse)
library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(tidygeocoder)
library(lubridate)
#> Loading required package: timechange
#>
#> Attaching package: 'lubridate'
#> The following objects are masked from 'package:base':
#>
#> date, intersect, setdiff, union
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 with s2
hemisphere <- s2::s2_buffer_cells(
s2::as_s2_geography(paste0("POINT(", lon, " ", lat, ")")),
distance = 1e7,
max_cells = 5000)
# 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")) %>%
# Use s2 for the cut
st_as_s2() %>%
s2::s2_intersection(hemisphere) %>%
# Back to sf
st_as_sf() %>%
st_transform(crs = projString) %>%
filter(!is.na(st_is_valid(.))) %>%
mutate(geometry = geometry * flip) %>%
# Filter if empty, since the cut can produce empty geometries
filter(!st_is_empty(.))
)
)
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 48 features and 0 fields
#> Geometry type: GEOMETRY
#> Dimension: XY
#> Bounding box: xmin: -8700015 ymin: -8913303 xmax: 8922028 ymax: 8998639
#> CRS: proj=laea x_0=0 y_0=0 lon_0=73.9726 lat_0=10.488
#> First 10 features:
#> geometry
#> 1 MULTILINESTRING ((3542468 3...
#> 2 LINESTRING (-6234955 -52010...
#> 3 MULTILINESTRING ((8922028 -...
#> 4 LINESTRING (3098546 2071855...
#> 5 MULTILINESTRING ((-1307725 ...
#> 6 LINESTRING (557848.4 -59059...
#> 7 MULTILINESTRING ((-24783.5 ...
#> 8 MULTILINESTRING ((-6264812 ...
#> 9 MULTILINESTRING ((-2356827 ...
#> 10 LINESTRING (-4432439 -32157...
ggplot(df)
geom_sf()
Created on 2023-01-23 with reprex v2.0.2