# Creating polygon with one third of world surface using R

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?

• Take a line from the N pole down 0 longitude to the south pole, then up back to the N pole along the line of 120 longitude. This will be 1/3 of the spherical earth. How well this area approximates to 1/3 the spherical area when projected is tricky though... Jun 19 at 9:51
• yes that sounds about right! I confess that I have a hard time to image how I would draw such a line, when not on a globe 😬
– Lenn
Jun 19 at 10:03

## 1 Answer

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.

• thank you so much! why would I need a high number for `n`? For me it works with only ten points too
– Lenn
Jun 19 at 11:46
• It depends on what you are going to do with this polygon. If you reproject it then all the points get joined by straight lines, for example if you only had one point on each meridian the reprojected version in Mercator would be like a square at 45 degrees. Jun 19 at 12:13
• ah yes:) that is true. Thanks so much !
– Lenn
Jun 19 at 12:42