0

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:

Satellite image around South Africa

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.

2
  • Hi Nicolas, it looks like you're preparing to export a lot of small chunks of imagery. I wonder if there are ways to do more processing in EE and then export larger chunks. What it is you are trying to accomplish in the end? Maybe I can suggest ways you might do more of the processing in EE. I don't mean to be presumptuous, just looking for opportunities to help :) May 28 '21 at 18:36
  • Hi Justin, thanks for the suggestion. I will generate a lot of satellite imagery in GEE to then use it to train a Convolutional Neural Network. I will also do more stuff in R with the grid, like computing some distances and intersecting the grid with other rasters, so I don't think I can do much more in GEE. May 28 '21 at 20:17
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));

enter image description here

1
  • Thank you very much for your answer! I had no idea how to pass those projection and error margin parameters to my functions. I'm now sure that the grid has the intended size, and I can add these parameters to my function that exports satellite imagery. May 28 '21 at 20:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.