ReducetoVectors in Google Earth Engine not working with the Mollweide transformation. Need to calculate zonal statistics for global scale analysis

Question

My study is a global-scale analysis of the factors driving forest gain, with grid cells as the unit of analysis. I need to calculate the amount of forest gain per ~65km2 hexagonal grid cell at the global scale but I've run into problems with the projection system and the image to vector conversion.

I used the code developed by Noel Gorelick to construct a hexagonal grid over the whole world (https://gorelick.medium.com/more-buffered-samples-with-hex-cells-b9a9bd36120d).The grid cells need to be equal area and I want to use the Mollweide equal area projected co-ordinate system to construct the grid.

To calculate the area of forest gain per grid cell, I have tried using reduceToVectors to convert the hexagonal grid image to vector so that then I can use ReduceRegions. However, the Image to Vector conversion is failing (it works when I tried it at the smaller scale, for the island of Madagascar). It gives the error shown in the following image. It seems it is not correctly incorporating the crs: proj command... I have also tried to do the zonal statistics keeping the hexagonal grid as an image and using reduceConnectedComponents but I'm not sure this is right. And when I try to export the resulting grid, it will not do it in the desired co-ordinate system. I get errors saying it cannot parse the co-ordinate system.

Any advice? It seems it is fundamentally a problem with the Mollweide projection system I am using. I'm also open to other suggestions for doing these zonal statistics. I am also very new to GEE.

var Forest_dynamics = ee.Image('projects/glad/GLCLU2020/Forest_type');
// Extract only the forest gain pixels

var Forest_gain = Forest_dynamics.updateMask(Forest_dynamics.eq(3));

// Reclassify to 1

var Forest_gainRC = Forest_gain.remap([3], [1]);

/** Create a hexagonal grid within which to extract values of the amount of forest gain.
* 
// Load map of all the countries in the world. 

var region = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
    .union();
          //Union merges them all into one feature layer representing Earth's land masses.

//Map.addLayer(region, {}, 'Region');

var hexGrid = function(proj, diameter, strict) {
  var size = ee.Number(diameter).divide(Math.sqrt(3)) // Distance from center to vertex
  
  var coords = ee.Image.pixelCoordinates(proj)
  var vals = {
    // Switch x and y here to get flat top instead of pointy top hexagons.
    x: coords.select("x"),
    u: coords.select("x").divide(diameter),  // term 1
    v: coords.select("y").divide(size),      // term 2
    r: ee.Number(diameter).divide(2),
  }
  var i = ee.Image().expression("floor((floor(u - v) + floor(x / r))/3)", vals)
  var j = ee.Image().expression("floor((floor(u + v) + floor(v - u))/3)", vals)
  
  // Turn the hex coordinates into a single "ID" number.
  var cells = i.long().leftShift(32).add(j.long()).rename("hexgrid")
  return cells
}

// Projection for World Mollweide
var wkt = ' \
  PROJCS["World_Mollweide", \
    GEOGCS["GCS_WGS_1984", \
      DATUM["WGS_1984", \
        SPHEROID["WGS_1984",6378137,298.257223563]], \
      PRIMEM["Greenwich",0], \
      UNIT["Degree",0.017453292519943295]], \
    PROJECTION["Mollweide"], \
    PARAMETER["False_Easting",0], \
    PARAMETER["False_Northing",0], \
    PARAMETER["Central_Meridian",0], \
    UNIT["Meter",1], \
    AUTHORITY["EPSG","54009"]]';

var proj = ee.Projection(wkt);

// Generate a hex grid and mask off cells that don't touch the region. This removes grid cells in the sea.

var grid = hexGrid(proj, 65657, true)
var regionImg = ee.Image(0).byte().paint(region, 1)
var mask = grid.addBands(regionImg)
    .reduceConnectedComponents(ee.Reducer.max(), "hexgrid", 128)
grid = grid.updateMask(mask)

Map.addLayer(grid.randomVisualizer(), {}, 'grid')
 
// To convert to polygon need to set the geometry to global. To do this create a bounding box variable:

var globalBoundingBox = ee.Geometry.Rectangle([-180, -90, 180, 90]);  // I have also tried this using the bounds for the Mollweide projection system from https://epsg.io/54009 and ee.Geometry.Rectange([-18040095.7, -9020047.85, 18040095.7, 9020047.85], proj, false). Also drew errors.     

Map.addLayer(globalBoundingBox, {}, 'Global bounding box')


var vector = grid.reduceToVectors({
   geometry: globalBoundingBox,
   crs: proj,
   scale: 1000,
   maxPixels: 1e16,
   eightConnected: true
 })


*** These is the alternative attempt using reduceConnectedComponents and keeping both layers as images ***

// Add the grid cell IDs to the Forest Gain layer as another band

var Forest_gain_grid = Forest_gainRC.addBands(grid);

var Gain_per_hex = Forest_gain_grid.reduceConnectedComponents(ee.Reducer.sum(), "hexgrid", 256)

Map.addLayer(Gain_per_hex.randomVisualizer(), {}, 'Gain per hex')

Export.image.toDrive({
  image: Gain_per_hex,
  description: 'Forest_Gain_per_hex',
  region: globalBoundingBox,
  maxPixels:1e13,
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true
  }
}); 

Answer 1

When you specify a polygon the way you did, you get a geodesic polygon, which has curved top and bottom, that essentially wraps up over the pole and back down, resulting in essentially a 0-size box.

Instead, you want a non-geodesic polygon, and since you're doing forest and don't care about the north or south pole, you should only use up to something like 88 degrees, so you don't have to deal with the singularities at the poles.

If you're going -180 to +180, tt's also always a good idea to use a 6-point polygon instead of a 2 -corner rectangle, to guarantee that it doesn't collapse as being degenerate.

var globalBoundingBox = ee.Geometry.Polygon([-180, -88, 0, -88, 180, -88, 180, 88, 0, 88, -180, 88], null, false)