Eclipse map
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
I saw this post showing a map of the year of the most recent total eclipse, and people mentioning that we can find the data on the Five Millennium Canon of Solar Eclipses Database (the data also mentioned in the thread on ArcGIS don’t go before the seventeenth century).
There is a bit of involved scraping because manually we can’t get more than 4 eclipses in a go and part of the generation is driven by javascript, but we can eventually get the whole 2538 eclipses between -1999 and 2024!
Setup
library(tidyverse) library(rvest) library(glue) library(httr) library(sf) library(mapview)
Data
First we make the query on the site: all total eclipse. It populates a select HTML control where we can manually scrape the eclipses dates with the browser tools. I saved them in a CSV file in the data directory.
We will next chose an area of interest (AOI). I’ll open a spatial data of Metropolitan France (built by merging all regions) ; pick your own, in EPSG:4326…
You can use data from Natural Earth with {rnaturalearth} for example. Add a character field eclipse_date to hold the dates of eclipse.
eclipses_dates <- read_csv("data/eclipses_dates.csv",
col_types = "c") |>
pull(astro_dates)
aoi <- read_sf("~/data/adminexpress/adminexpress_cog_simpl_000_2022.gpkg",
layer = "region") |>
filter(insee_reg > "06") |>
st_union() |>
tibble(geom = _) |>
st_sf() |>
mutate(eclipse_date = NA_character_, .before = 1)
Download files
Based on a half-hidden URL, we will, slowly, ask for the KMZ to be generated, find its URL, get the file and save it in the results directory. It should take an hour…
get_kmz <- function(eclipse_date) {
message(glue("Downloading: {eclipse_date}"))
tryCatch({
read_html(glue("http://xjubier.free.fr/en/site_pages/solar_eclipses/xSE_GoogleEarth.php?Ecl={eclipse_date}&Acc=2&Umb=1&Lmt=1&Mag=0")) |>
html_elements("fieldset a") |>
html_attr("href") |>
GET(write_disk(glue("results/eclipse_{eclipse_date}.kmz")))
},
error = function(e) {
message(" x Can't download/save")
print(e)
return(NULL)
})
}
eclipses_dates |>
walk(slowly(get_kmz), .progress = TRUE)
If you are eager to start, get a bundle of all of them from here (110 MB). Uncompress in results.
Prepare the processing
We need to open the right layer in each KMZ, as there are many of them with varying names, then we intersect it with our AOI, recursively. So we make one function for each task…
Many layers have geometry errors for the S2 engine: we will skip them ; so beware the map may be inaccurate!
get_polygon <- function(eclipse_date) {
tryCatch({
layer <- st_layers(glue("results/eclipse_{eclipse_date}.kmz")) |>
pluck("name") |>
keep(\(x) str_detect(x, "Umbral Path"))
read_sf(glue("results/eclipse_{eclipse_date}.kmz"),
layer = layer) |>
st_zm() |>
st_make_valid() |>
st_collection_extract("POLYGON") |>
transmute(eclipse_date = eclipse_date)
},
error = function(e) {
message(" x Error in opening")
print(e)
return(NULL)
})
}
crop_map <- function(current_map, eclipse_date) {
message(glue("Date: {eclipse_date}"))
p <- get_polygon(eclipse_date)
tryCatch({
if (any(apply(st_intersects(current_map, p), 1, any))) {
message(" -> Eclipse matches area of interest")
bind_rows(st_difference(current_map, p) |>
select(eclipse_date),
st_intersection(current_map, p) |>
mutate(eclipse_date = eclipse_date.1) |>
select(-eclipse_date.1))
} else {
current_map
}
},
error = function(e) {
message(" x Error in intersection")
print(e)
return(current_map)
})
}
And we run them:
aoi_eclipse <- eclipses_dates |>
reduce(crop_map, .init = aoi) |>
write_sf("aoi_eclipse.gpkg")
You can use the Geopackage in QGIS, or just display it here:
eclipse_year <- aoi_eclipse |> st_make_valid() |> group_by(year = str_sub(eclipse_date, 1, 5)) |> summarise() mapview(eclipse_year)
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.