Oblique Geo Coords with GEE in R

Function to create Oblique geographic coordinates (OGC) in R from a Feature/FeatureCollection (or shapefile) with Google Earth Engine. This function is basically all taken from Anders Bjørn Møller and was developed in Møller et al 2020. This code is just simply converted for use in the GEEs framework and the original authors should be cited.

Load rgee

library(rgee)
ee_Initialize()

## ── rgee 1.1.5 ─────────────────────────────────────── earthengine-api 0.1.323 ── 
##  ✔ user: not_defined
##  ✔ Initializing Google Earth Engine:
 ✔ Initializing Google Earth Engine:  DONE!
## 
 ✔ Earth Engine account: users/trevanflynn 
## ────────────────────────────────────────────────────────────────────────────────

Only need sf library if you want to use your own shapefile. This example will be done with a FeatureCollection straight from GEE. If you use your own shapefile, import as shown below and skip to next section.

For your own use:

region = read_sf(“your_area.shp”)%>%st_agr(‘constant’)%>%sf_as_ee()

#province from GEE dataset
#Get province
worldProvince = ee$FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level1"); 
filterProvince = ee$Filter$eq('ADM1_NAME', 'Limpopo'); #filter province
region = worldProvince$filter(filterProvince); #get 

#Check data
Map$centerObject(region, 7)
Map$addLayer(region)

This code only needs a region and no actual raster image. The area defines the bounds which can be manipulated into coordinate systems and resolution. However, it is recommended that a projected coordinate system in m be used.

ogcCalc() function: Calculates OGC from a Geometry, Feature or FeatureCollection.

#function to get an image of coordinates from FeatureCollection
ogcCalc = function(area, crs = "EPSG: 3857", scale = 250, nr = 6){
  
  #create x and y raster
  img = ee$Image$pixelCoordinates(crs)$ #x and y from coordinates
    clip(area)$ #clip to region
    reproject(crs =crs, scale = scale)
  
  #select x ane y
  x = img$select("x")
  y = img$select("y")
  
  #angles and position
  pis <- pi*seq(0, 1, 1/nr)[2:(nr)] #sequence of angles
  pos <- seq(2, nr, 1)
  
  #need to edit locations for even numbers 
  is.even <-function(x, tol = .Machine$double.eps^0.5){
    abs(x - round(x)) < tol
    }
  
  if(is.even(nr/2)) {
    pis <- pis[ -(nr/2)]
    pos <- pos[ -(nr/2)]
  }
  
  #Make a list for each image calculation
  calc1 = list()
  calc2 = list()
  ogc = list() #list of ogc
  
  #calculation 1
  for(i in 1:length(pis)){
    
    #get constant image of angles
    calc1[[pos[i]]] = ee$Image$constant(ee$Number(pis[i]))$
      clip(region)
    
    #calc 1
    calc1[[pos[i]]] = calc1[[pos[i]]]$subtract((y$divide(x))$atan())$cos()
    
    #calc 2
    calc2 = (x$pow(2)$add(y$pow(2)))$sqrt()
    
    #calc 3
    ogc[[pos[i]]] = calc1[[pos[i]]]$multiply(calc2)
  }
  
  #make sure to account for even numbers (will just be null otherwise)
  ogc[[1]] = x
  ogc[[(nr/2)+1]] = y
  
  #create labels (cant have "." in GEE names so substitute out for "_" or just make sure there are no decimals)
  labs=  paste0("pi", format(seq(0, 1, 1/nr)[1:nr], digits = 1, nsmall = 2))%>%
    gsub("[.]", "_",.)
  
  #add xy back in (otherwise names will be wrong)
  ogc = ee$ImageCollection(ogc)$
    toBands()$
    rename(labs)
  
  #return OGC
  return(ogc)
}

Now run the function on the xyImage and they are ready.

#apply function
ogcImg = ogcCalc(region, crs = "EPSG:4326")

#check to see if worked (first and last are x and y)
print(ogcImg$bandNames()$getInfo())

## [1] "pi0_00" "pi0_17" "pi0_33" "pi0_50" "pi0_67" "pi0_83"

#print to see x axis
#get min and mx for visual
mn = ogcImg$reduceRegion(reducer = 'min', geometry = region, maxPixels = 1e13)
mx = ogcImg$reduceRegion(reducer = 'max', geometry = region, maxPixels = 1e13)

#visParam
vis = list(min = mn$get("pi0_00"), max = mx$get("pi0_00"), palette = c("blue", "red", "yellow"))

#Plot to make sure it worked
Map$addLayer(ogcImg$select("pi0_00"), vis)

#print to see transformation
#visParam
vis = list(min = mn$get("pi0_17"), max = mx$get("pi0_17"), palette = c("blue", "red", "yellow"))

#Plot to make sure it worked
Map$addLayer(ogcImg$select("pi0_17"), vis)

Very good covariates for DSM and almost always comes out as some of the most important in most of my models, especially for large areas.

REFERENCE

Møller, A. B., Beucher, A. M., and Pouladi, N., Greve, M. H. (2020) Oblique geographic coordinates as covariates for digital soil mapping, SOIL 6(2) 269 - 289 https:/10.5194/soil-6-269-2020