Eclipse map

R-bloggers 2024-04-12

[This article was first published on r.iresmi.net, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

2017 Solar Eclipse with Totality – Composite – CC-BY-NC by Jeff Geerling

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…

Note

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…

Warning

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)

To leave a comment for the author, please follow the link and comment on their blog: r.iresmi.net.

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.

Continue reading: Eclipse map