feat: add grid management

geetools/ee_feature.py

@@ -14,6 +14,57 @@ def __init__(self, obj: ee.Feature):
         """Initialize the class."""
         self._obj = obj
 
+    @classmethod
+    def cell(
+        cls,
+        size: int | ee.Number,
+        ix: int | ee.Number = 0,
+        iy: int | ee.Number = 0,
+        crs: str | ee.String = "EPSG:4326",
+    ) -> ee.Feature:
+        """Create a feature representing one cell in a regularly spaced grid.
+
+        Args:
+            size: The cell side length in the units of the CRS.
+            ix: The cell index in the x direction, starting at the CRS area-of-use minimum.
+            iy: The cell index in the y direction, starting at the CRS area-of-use minimum.
+            crs: The CRS of the grid.
+
+        Returns:
+            A feature representing the requested grid cell.
+        """
+        from pyproj import CRS, Transformer
+
+        size, crs = ee.Number(size), ee.String(crs)
+        crs_value = crs.getInfo()
+        # Earth Engine cannot currently compute projection bounds server-side; pyproj supplies the area of use.
+        # https://gis.stackexchange.com/questions/471100/how-to-get-the-bounds-of-a-projection-in-gee
+        # https://issuetracker.google.com/issues/313937471
+        crs_x_min, crs_y_min, crs_x_max, crs_y_max = CRS(crs_value).area_of_use.bounds
+        transformer = Transformer.from_crs(4326, crs_value, always_xy=True)
+        x_min, y_min = transformer.transform(crs_x_min, crs_y_min)
+        x_max, y_max = transformer.transform(crs_x_max, crs_y_max)
+
+        x_min, y_min = ee.Number(x_min).toFloat(), ee.Number(y_min).toFloat()
+        x_max, y_max = ee.Number(x_max).toFloat(), ee.Number(y_max).toFloat()
+        longitudes = ee.List.sequence(x_min, x_max, size).add(x_max).distinct()
+        latitudes = ee.List.sequence(y_min, y_max, size).add(y_max).distinct()
+        ix, iy = ee.Number(ix).toInt(), ee.Number(iy).toInt()
+
+        cell_id = ee.String("{proj}_{size}_E{ix}-S{iy}").geetools.format(
+            {
+                "proj": ee.Projection(crs).crs(),
+                "size": size.format(),
+                "ix": ix.format("%02d"),
+                "iy": iy.format("%02d"),
+            }
+        )
+        rectangle = ee.Geometry.Rectangle(
+            [longitudes.get(ix), latitudes.get(iy), longitudes.get(ix.add(1)), latitudes.get(iy.add(1))],
+            crs,
+        )
+        return ee.Feature(rectangle, {"id": cell_id})
+
     def toFeatureCollection(self) -> ee.FeatureCollection:
         """Convert a :py:class:`ee.Feature` composed of a multiGeometry geometry into a :py:class:`ee.FeatureCollection`.
 

geetools/ee_feature_collection.py

@@ -8,6 +8,7 @@
 import shapely
 from matplotlib import pyplot as plt
 from matplotlib.axes import Axes
+from pyproj import CRS, Transformer
 
 from .accessors import register_class_accessor
 from .utils import plot_data
@@ -30,6 +31,69 @@ def __init__(self, obj: ee.FeatureCollection):
         """Initialize the :py:class:`ee.FeatureCollection` class."""
         self._obj = obj
 
+    @classmethod
+    def grid(cls, size: int | float | ee.Number, crs: str | ee.String = "EPSG:4326") -> ee.FeatureCollection:
+        """Create a global feature grid with square cells in the requested CRS.
+
+        Warning:
+            Generating this grid is a very costly operation. It is strongly recommended to export the
+            resulting feature collection to an Earth Engine asset before using it in other computations.
+
+        Args:
+            size: The cell side length in the units of the CRS.
+            crs: The CRS of the grid. Defaults to ``"EPSG:4326"``.
+
+        Returns:
+            A feature collection whose cells have an ``id`` property.
+        """
+        size, crs = ee.Number(size), ee.String(crs)
+        crs_value = crs.getInfo()
+        # Earth Engine cannot currently compute projection bounds server-side; pyproj supplies the area of use.
+        # https://gis.stackexchange.com/questions/471100/how-to-get-the-bounds-of-a-projection-in-gee
+        # https://issuetracker.google.com/issues/313937471
+        crs_x_min, crs_y_min, crs_x_max, crs_y_max = CRS(crs_value).area_of_use.bounds
+        transformer = Transformer.from_crs(4326, crs_value, always_xy=True)
+        x_min, y_min = transformer.transform(crs_x_min, crs_y_min)
+        x_max, y_max = transformer.transform(crs_x_max, crs_y_max)
+
+        distance_transformer = Transformer.from_crs(crs_value, 3857, always_xy=True)
+        x1 = distance_transformer.transform(x_min, 0)[0]
+        x2 = distance_transformer.transform(x_min + size.getInfo(), 0)[0]
+        equator_cell_size = int(abs(x2 - x1))
+
+        x_min, y_min = ee.Number(x_min).toFloat(), ee.Number(y_min).toFloat()
+        x_max, y_max = ee.Number(x_max).toFloat(), ee.Number(y_max).toFloat()
+        longitudes = ee.List.sequence(x_min, x_max, size).add(x_max).distinct()
+        latitudes = ee.List.sequence(y_min, y_max, size).add(y_max).distinct()
+        longitude_indexes = ee.List.sequence(0, longitudes.size().subtract(2))
+        latitude_indexes = ee.List.sequence(0, latitudes.size().subtract(2))
+
+        def rectangle(i: ee.Number, j: ee.Number) -> ee.Feature:
+            i, j = ee.Number(i).toInt(), ee.Number(j).toInt()
+            coordinates = [
+                longitudes.get(i),
+                latitudes.get(j),
+                longitudes.get(i.add(1)),
+                latitudes.get(j.add(1)),
+            ]
+            cell_id = ee.String("{proj}_{size}_E{ix}-S{iy}").geetools.format(
+                {
+                    "proj": ee.Projection(crs).crs(),
+                    "size": size.format(),
+                    "ix": i.format("%02d"),
+                    "iy": j.format("%02d"),
+                }
+            )
+            return ee.Feature(
+                ee.Geometry.Rectangle(coordinates, crs),
+                {"id": cell_id, "equator_cell_size_m": equator_cell_size},
+            )
+
+        def column(i: ee.Number) -> ee.List:
+            return latitude_indexes.map(lambda j: rectangle(i, j))
+
+        return ee.FeatureCollection(ee.List(ee.List(longitude_indexes).map(column)).flatten())
+
     def toImage(
         self,
         color: str | ee.String | int | ee.Number = 0,

pyproject.toml

@@ -29,6 +29,7 @@ dependencies = [
     "requests",
     "pandas",
     "geopandas",
+    "pyproj",
     "deprecated",
     "xee>=0.1.0", # xee change API prototype for geometry setting
     "yamlable",

tests/test_Feature.py

@@ -4,6 +4,29 @@
 import geetools  # noqa: F401
 
 
+class TestCell:
+    """Test the ``cell`` method."""
+
+    def test_cell(self):
+        cell = ee.Feature.geetools.cell(20)
+        assert cell.get("id").getInfo() == "EPSG:4326_20_E00-S00"
+        assert cell.geometry().coordinates().getInfo() == [
+            [[-180, -90], [-160, -90], [-160, -70], [-180, -70], [-180, -90]]
+        ]
+
+    def test_cell_with_indexes(self):
+        cell = ee.Feature.geetools.cell(20, ix=5, iy=5)
+        assert cell.get("id").getInfo() == "EPSG:4326_20_E05-S05"
+
+    def test_cell_with_size_10_and_indexes(self):
+        cell = ee.Feature.geetools.cell(10, ix=10, iy=10)
+        assert cell.get("id").getInfo() == "EPSG:4326_10_E10-S10"
+
+    def test_cell_with_projected_crs(self):
+        cell = ee.Feature.geetools.cell(1000, crs="EPSG:3857")
+        assert cell.get("id").getInfo() == "EPSG:3857_1000_E00-S00"
+
+
 class TestToFeatureCollection:
     """Test the ``toFeatureCollection`` method."""
 

tests/test_FeatureCollection.py

@@ -33,6 +33,25 @@ def vatican(self):
         return ee.Geometry.Point([12.453386, 41.903282]).buffer(1)
 
 
+class TestGrid:
+    """Test the ``grid`` method."""
+
+    def test_grid_size_and_ids(self):
+        grid = ee.FeatureCollection.geetools.grid(20)
+        assert grid.size().getInfo() == 162
+        assert grid.first().get("id").getInfo() == "EPSG:4326_20_E00-S00"
+
+    def test_grid_with_float_size(self):
+        grid = ee.FeatureCollection.geetools.grid(60.0)
+        assert grid.size().getInfo() == 18
+        assert grid.first().get("id").getInfo() == "EPSG:4326_60_E00-S00"
+
+    def test_grid_projected_crs(self):
+        grid = ee.FeatureCollection.geetools.grid(10_000_000, "EPSG:3857")
+        assert grid.size().getInfo() > 0
+        assert grid.first().get("id").getInfo() == "EPSG:3857_10000000_E00-S00"
+
+
 class TestReduceColumnsGroup:
     """Test the ``reduceColumnsGroup`` method."""
 

tests/test_FeatureCollection/serialized_test_split.yml

@@ -174,15 +174,15 @@ values:
         count:
           functionInvocationValue:
             arguments:
-              left:
+              input:
                 functionInvocationValue:
                   arguments:
-                    input:
+                    left:
                       constantValue: 3
-                  functionName: Number.toInt
-              right:
-                constantValue: 1
-            functionName: Number.add
+                    right:
+                      constantValue: 1
+                  functionName: Number.add
+            functionName: Number.toInt
         end:
           constantValue: 1
         start: