Creating polygon with one third of world surface using R

Question

I would like to create a polygon with (roughly) one third (or any other fraction) of the worlds surface.

I assume I should use an equal area projeciton like Lambert cylindrical equal-area, robinson or maybe mollweide.

What I did (in R) was the following.

I went to epsg.io, grabbed the extents of the robinson-projection (54030) and did this:

world_bb = c(-17005833.33, -8625154.67, 17005833.33, 8625154.67) %>% setNames(c("xmin", "ymin", "xmax", "ymax")) %>% st_bbox() %>% 
  st_as_sfc() %>% st_as_sf(crs="+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs")
world_area_sqkm = st_area(world_bb) %>% units::set_units("km^2")
world_area_sqkm
#586711771 [km^2]

Google tells me that the world has a surface area of 510 million square kilometers.

Also, when I reproject that bounding box to 4326 to plot it on a web map I get this:

world_bb_4326 = st_transform(world_bb, 4326)
mapview(world_bb_4326)

Does anyone have any idea on how to do that better using R?

Answer 1

Write a function to return a segment of earth surface between two longitudes. Straightforward except a few complications - need to make sure the coordinates match at the poles and need to fill in lots of points along the meridians (n) to get accurate lines:

eslice <- function(n, l1, l2){
    latdown = seq(90,-90, length=n)
    londown = rep(l1, length(latdown))
    latup = seq(-90, 90, length=n)[-1]
    lonup = rep(l2, length(latup))

    lonlat = cbind(c(londown, lonup), c(latdown, latup))
    lonlat[nrow(lonlat),] = lonlat[1,]
    sf::st_sfc(sf::st_polygon(list(lonlat)), crs=4326)
}

Then you can:

s1 = eslice(100, 0, 120)
mapview(s1)

Its area will be computed by spherical geometry, so should be 1/3 the spherical earth surface area:

> st_area(s1) * 3
5.100661e+14 [m^2]

Yup.