Reading Assets from Google Earth Engine with GDAL and R {terra}

A guide on how to directly read assets from Google Earth Engine

Intro and Aim

Google Earth Engine (GEE) stores petabytes of Earth Observation data and more. The great aspect is that you can perform analyses directly in the cloud. This means that, in principle, you only need to extract your results (summary statistics, plots, final maps, etc.) from GEE. However, there are times when you want to conduct part of your analysis using your own system. For these instances, accessing (or downloading) raw or preprocessed data from GEE becomes necessary.

Recently, XEE emerged as a potential game-changer. But there’s another, perhaps overlooked, method for accessing GEE data directly from your machine: using GDAL and any software built on it, like the R {terra} package.

I was inspired by Andrew Brown’s blog post and decided to try it out myself.

Stack Requirements

Ensure you have the following:

• GEE account (and a connected Google Cloud project, created with your GEE account)
• GDAL
• R and the {terra} package

Setup

Follow these steps to access GEE data from your machine:

1. Go to https://console.cloud.google.com/
2. Select the cloud project you use with Google Earth Engine.
3. Navigate to APIs and Services and enable the Google Earth Engine API.
4. Create a service account.
5. Generate a credentials .json file and save it on your machine.

GDAL (on Linux)

With everything set up, we can access the data. GDAL, as a command-line tool, is both convenient and powerful.

I chose the DEM Merit dataset for this example. It’s a single, global image with one band (“dem”), simplifying the process for demonstration purposes.

First, I create an environment variable with the location of my credentials file.

Then, I use gdalinfo on the dataset location, adding the “EEDAI:” prefix. This tells GDAL that the dataset is stored on Google Earth Engine and to use the appropriate protocol.

export GOOGLE_APPLICATION_CREDENTIALS="$HOME/gee/ee-credentials.json"
gdalinfo EEDAI:MERIT/DEM/v1_0_3

Driver: EEDAI/Earth Engine Data API Image
Files: none associated
Size is 432000, 174000
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-180.000416666666666,84.999583333333334)
Pixel Size = (0.000833333333333,-0.000833333333333)
Corner Coordinates:
Upper Left  (-180.0004167,  84.9995833) (180d 0' 1.50"W, 84d59'58.50"N)
Lower Left  (-180.0004167, -60.0004167) (180d 0' 1.50"W, 60d 0' 1.50"S)
Upper Right ( 179.9995833,  84.9995833) (179d59'58.50"E, 84d59'58.50"N)
Lower Right ( 179.9995833, -60.0004167) (179d59'58.50"E, 60d 0' 1.50"S)
Center      (  -0.0004167,  12.4995833) (  0d 0' 1.50"W, 12d29'58.50"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = dem
  Overviews: 216000x87000, 108000x43500, 54000x21750, 27000x10875, 13500x5437, 6750x2718, 3375x1359, 1687x679, 843x339, 421x169, 210x84

With this, you’ll get information about the dataset.

Once confirmed it’s working, you can use other GDAL tools like gdalwarp, gdal_translate, etc.

R Terra

{terra} uses GDAL to read raster data, meaning the same dataset can be read into R just as we did with GDAL.

Similar to the GDAL example, we set an environment variable for the credentials file and then read the dataset (with “EEDAI:”) using the rast() function.

library(terra)

terra 1.7.37

Sys.setenv(GOOGLE_APPLICATION_CREDENTIALS=file.path(path.expand('~'),"gee/ee-credentials.json"))

rst <- rast("EEDAI:MERIT/DEM/v1_0_3")

Now, the asset can be accessed with all the functions available in {terra}. We can print the object and plot it. Plotting is fast, even for a global layer at 90 meters resolution, because terra only retrieves a sample of it, just enough to provide a clear image at the needed resolution.

print(rst)

class       : SpatRaster 
dimensions  : 174000, 432000, 1  (nrow, ncol, nlyr)
resolution  : 0.0008333333, 0.0008333333  (x, y)
extent      : -180.0004, 179.9996, -60.00042, 84.99958  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : v1_0_3 
name        : dem 

plot(rst, range=c(0,5000)) # wide range for global DEM

Next, we select a spatial subset using the window() function by defining a bounding box.

xmin = 28 # bounding box xmin
xmax = 29 # bounding box xmax
ymin = -3 # bounding box ymin
ymax = -2 # bounding box ymax

bb = ext(c(xmin, xmax, ymin, ymax))
window(rst) <- NULL
window(rst) = bb
print(rst)

class       : SpatRaster 
dimensions  : 1200, 1200, 1  (nrow, ncol, nlyr)
resolution  : 0.0008333333, 0.0008333333  (x, y)
window      : 28.00042, 29.00042, -2.999583, -1.999583  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : v1_0_3 
name        : dem 

We can now obtain summary statistics from the area with the summary() function.

summary(rst)

Warning: [summary] used a sample

      dem        
 Min.   : 641.1  
 1st Qu.:1175.5  
 Median :1448.5  
 Mean   :1508.9  
 3rd Qu.:1750.8  
 Max.   :3431.0  

Finally, plot the area with an appropriate value range.

plot(rst, range=c(650,3500))