Adequate CRS when creating a grid over a whole continent in R

Question

I'm trying to generate a grid of 6.72 km by 6.72 km points over the whole continent of Africa, but I'm having a lot of trouble in doing so.

I'm mainly running into 2 problems:

My first issue arises when I try to generate a grid on a single country. I downloaded the country level map of the world from GADM, kept only countries in Africa, and then I used the "Simplify Tool" In QGIS to make the polygons more simple so the code runs in a reasonable time.

Here is some example code for Angola based heavily on this answer on Stack Overflow, where I work with the original non-simplified GADM map, so I also changed the grid size to 200km so the code runs in a reasonable time. I had to change the CRS to EPSG:26910 for the code to run.

library(sf)
library(tidyverse)

#importing map of Angola
angola <- readRDS(gzcon(url("https://biogeo.ucdavis.edu/data/gadm3.6/Rsf/gadm36_AGO_0_sf.rds")))
#checking shape
plot(angola$geometry)
#checking CRS
st_crs(angola)
#changing CRS
angola <- st_transform(angola, 26910)
st_crs(angola)

#building grid
grid <- angola %>%
  st_make_grid(cellsize =200000) %>%
  st_intersection(angola) %>%
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>%
  mutate(id = row_number())

#checking grid
plot(grid$geometry)

The problem here with the CRS is that the map with the grid is flipped. I tried changing back the CRS of the grid to the original after doing the calculations, but the dimensions of the pixels are not the desired when going back to the original CRS.

What CRS is the most adequate for this task? I'm very new to GIS, but I need a CRS that works properly in all Africa and that allows me to define a grid in meters and not in degrees.

My second issue is that I can't generate a grid over the whole continent of Africa. I downloaded the country level map of the world from GADM, kept only countries in Africa, and then I used the "Dissolve Tool" in QGIS to combine the shapes of all the countries and get the shape of the continent.

I can't provide code for people to replicate what I did in QGIS, so I attach this Dropbox link with the continent shapefile, and an example again with a grid with 200km, and that obviously suffer from the CRS problems mentioned in the previous point.

africa_comb <-st_read("africa_comb.shp")

#changing CRS
africa_comb <- st_transform(africa_comb, 26910)

#making the grid
grid_comb <- africa_comb %>%
  st_make_grid(cellsize =200000) %>%
  st_intersection(africa_comb) %>%
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>%
  mutate(id = row_number())

Here I get the following error:

Error in CPL_geos_op2(op, x, y) : 
  Evaluation error: TopologyException: Input geom 1 is invalid: Self-intersection at or near point...

I have tried before doing the same without dissolving the countries, but I run into the same error. Both maps look like this after changing the CRS, so I'm suspecting the CRS might be causing the issue here:

Answer 1

find here a full solution.

First, let's check if EPSG:26910 is suitable for Africa:

library(giscoR)
library(sf)
library(tidyverse)


# Wrong projection - full world
gisco_coastallines %>%
  st_transform(26910) %>%
  ggplot() +
  geom_sf()

It doesn't look very nice. If we check the crs we see the wkt string:

st_crs(26910)
#> Coordinate Reference System:
#>   User input: EPSG:26910 
#>   wkt:
#> PROJCRS["NAD83 / UTM zone 10N",
#>     BASEGEOGCRS["NAD83",
#>         DATUM["North American Datum 1983",
#>             ELLIPSOID["GRS 1980",6378137,298.257222101,
#>                 LENGTHUNIT["metre",1]]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         ID["EPSG",4269]],
#>     CONVERSION["UTM zone 10N",
#>         METHOD["Transverse Mercator",
#>             ID["EPSG",9807]],
#>         PARAMETER["Latitude of natural origin",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8801]],
#>         PARAMETER["Longitude of natural origin",-123,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8802]],
#>         PARAMETER["Scale factor at natural origin",0.9996,
#>             SCALEUNIT["unity",1],
#>             ID["EPSG",8805]],
#>         PARAMETER["False easting",500000,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8806]],
#>         PARAMETER["False northing",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8807]]],
#>     CS[Cartesian,2],
#>         AXIS["(E)",east,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1]],
#>         AXIS["(N)",north,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1]],
#>     USAGE[
#>         SCOPE["Engineering survey, topographic mapping."],
#>         AREA["North America - between 126°W and 120°W - onshore and offshore. Canada - British Columbia; Northwest Territories; Yukon. United States (USA) - California; Oregon; Washington."],
#>         BBOX[30.54,-126,81.8,-119.99]],
#>     ID["EPSG",26910]]

Well, AREA is pretty descriptive.

Using 3857

Here we select all the countries of Africa and create a single shapefile. I am using a 60m resolution, that is the lowest level provided on giscoR and combining all the shapes:

# Full solution
# Get Africa (low res = 60M)

africa <- gisco_get_countries(region = "Africa", resolution = 60) %>%
  st_union() %>%
  st_transform(3857)
#> although coordinates are longitude/latitude, st_union assumes that they are planar

ggplot(africa) +
  geom_sf()

That's much better.

Using crsuggest

Now we can use crsuggest to check a better crs

# Using crsuggest
library(crsuggest)
#> Using the EPSG Dataset v10.019, a product of the International Association of Oil & Gas Producers. 
#> Please view the terms of use at https://epsg.org/terms-of-use.html.

suggest_crs(africa)
#> # A tibble: 10 x 6
#>    crs_code crs_name      crs_type crs_gcs crs_units crs_proj4                  
#>    <chr>    <chr>         <chr>      <dbl> <chr>     <chr>                      
#>  1 32732    WGS 84 / UTM~ project~    4326 m         +proj=utm +zone=32 +south ~
#>  2 32532    WGS 72BE / U~ project~    4324 m         +proj=utm +zone=32 +south ~
#>  3 32332    WGS 72 / UTM~ project~    4322 m         +proj=utm +zone=32 +south ~
#>  4 32731    WGS 84 / UTM~ project~    4326 m         +proj=utm +zone=31 +south ~
#>  5 32531    WGS 72BE / U~ project~    4324 m         +proj=utm +zone=31 +south ~
#>  6 32331    WGS 72 / UTM~ project~    4322 m         +proj=utm +zone=31 +south ~
#>  7 5523     WGS 84 / Gab~ project~    4326 m         +proj=tmerc +lat_0=0 +lon_~
#>  8 5223     WGS 84 / Gab~ project~    4326 m         +proj=tmerc +lat_0=0 +lon_~
#>  9 26692    M'poraloko /~ project~    4266 m         +proj=utm +zone=32 +south ~
#> 10 32633    WGS 84 / UTM~ project~    4326 m         +proj=utm +zone=33 +datum=~

# We try 32732
africa_new <- africa %>% st_transform(32732)

ggplot(africa_new) +
  geom_sf()

Ok, it works too. Let's stick to it and make the grid

Grid (200 km)

# Make grid 200 km

grid <- st_make_grid(africa_new, cellsize = 200000) %>%
  st_intersection(africa_new) %>%
  st_cast("MULTIPOLYGON") %>%
  st_sf() %>%
  mutate(id = row_number())

ggplot(africa_new) +
  geom_sf() +
  geom_sf(data = grid)

^{Created on 2021-05-15 by the reprex package (v2.0.0)}