Generating grid that matches Google Earth Engine polygons

Question

I'm trying to create a grid in R over the whole continent of Africa (I received some help in this GIS StackExchange post), but I'm having trouble generating a grid that matches the CRS used in the Google Earth Engine, which I'm going to use later to download satellite imagery.

The discrepancies are higher when I depart from the Equator line. For instance, here is some code to generate a grid of 10 km by 10 km patches over South Africa, using the stars library:

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

# country shapefile
south_africa <- gisco_get_countries(country="ZAF", resolution = 10)  %>% st_transform(3857)

# method 1: stars
grid = st_as_stars(st_bbox(south_africa), dx = 10000, dy = 10000)
grid = st_as_sf(grid)
grid = grid[south_africa, ]

Here I used CRS 3857, so I have my units in meters to define the grid. I can compute the area of the pixels, and get the geometry of a pixel so we can see that the area is correct (1, because I normalized it) under this coordinate reference system:

grid$area = st_area(grid$geometry)/(10000)^2
print(paste('area:',grid[[3]][[10000]]))
grid[[2]][[10000]]

[1] "area: 1"
POLYGON ((2331832 -3517624, 2341832 -3517624, 2341832 -3527624, 2331832 -3527624, 2331832 -3517624))

I also checked the grid using st_make_grid, and I obtain exactly the same results when checking the same patch in the grid:

grid_2=st_make_grid(south_africa, cellsize = 10000)
grid_2 = st_as_sf(grid_2)
grid_2 = grid[south_africa, ]
grid_2$area = st_area(grid_2$geometry)/(10000)^2  
print(paste('area:',grid_2[[3]][[10000]]))
grid_2[[2]][[10000]]

[1] "area: 1"
POLYGON ((2331832 -3517624, 2341832 -3517624, 2341832 -3527624, 2331832 -3527624, 2331832 -3517624))

The problem is that if I take this polygon and import it to Google Earth Engine, I see that the dimensions of the polygon are not right. In this case, I get that the square root of the area is 8,637 meters, and the average of the length of the sides of the polygon is 8,632 meters, but both measures should be 10,000 meters if the geometries where properly defined.

I also take a Landsat-8 Satellite image around the centroid of the polygon, covering 10km by 10km, and I obtain an image that is bigger than my polygon. Here you can see the satellite image, and the polygon in black:

Here is the code to visualize this in the code editor of Google Earth Engine.

Is there any way in which I could fix this? I probably need to define a CRS that matches the one in Google Earth Engine. I tried working with CRS 4326, but I couldn't define a grid in meters with that CRS.

Answer 1

Whenever you are calculating areas and lengths, the projection is important to consider. Most ee.Geometry methods have parameters for projection and maxError; set them to specify what projection and tolerance to use for calculations.

In your script, I've added these arguments to the geometry methods and it computes the correct area and perimeter for the defined polygon and the buffered centroid matches the polygon (when buffer distance is defined in the desired 'EPSG:3857' projection).

Full Code Editor script

For example:

var MAXERROR = 0.001;
var PROJ = 'EPSG:3857';

var stars= ee.Geometry.Polygon([
  2331832, -3517624, 2341832, -3517624, 2341832,
  -3527624, 2331832, -3527624, 2331832, -3517624], PROJ, false);

print('square root of area', stars.area(MAXERROR, PROJ).sqrt());
print('average side size', stars.perimeter(MAXERROR, PROJ).divide(4));