How to get back a single pixel in the MODIS sinusoidal projection?

This seems like it should give me just one pixel because I'm asking for a 1x1 meter rectangle from a dataset whose scale is 463m:

# Not using SR-ORG:6974, as this is an EE-specific code and not all operations understand it
SIN = ("""PROJCS["MODIS Sinusoidal",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Sinusoidal"],
    PARAMETER["false_easting",0.0],
    PARAMETER["false_northing",0.0],
    PARAMETER["central_meridian",0.0],
    PARAMETER["semi_major",6371007.181],
    PARAMETER["semi_minor",6371007.181],
    UNIT["m",1.0],
    AUTHORITY["SR-ORG","6974"]]""")

ee.Initialize()

kwargs = { 
  'engine': 'ee',
  'scale': 463.3127165279165,
  'geometry': ee.Geometry.Rectangle(
      [4447800,-2223901,4448801,-2222900],
      ee.Projection(SIN), False, True),
  'crs': SIN,
}
image = ee.ImageCollection(ee.ImageCollection('MODIS/061/MCD43A4').first())
ds = xarray.open_dataset(image, **kwargs)
print(ds)

Problem 1: I get ee_exception.EEException: Geometry.bounds: Unable to perform this geometry operation. Please specify a non-zero error margin. because this call doesn't have an error margin rpcs.append(('bounds', self.geometry.bounds()))

If I change it to rpcs.append(('bounds', self.geometry.bounds(1))), I get

Problem 2: The code now works, but returns a 3x2 array, not 1x1.

<xarray.Dataset>
Dimensions:                                   (time: 1, X: 3, Y: 2)
Coordinates:
  * time                                      (time) datetime64[ns] 2000-02-24
  * X                                         (X) float32 42.6 42.61 42.61
  * Y                                         (Y) float32 -19.89 -19.89

I think this happens because of a double conversion from projection units (meters) to degrees and then back in

    # Parse the dataset bounds from the native projection (either from the CRS
    # or the image geometry) and translate it to the representation that will be
    # used for all internal `computePixels()` calls.
    try:
      if isinstance(geometry, ee.Geometry):
        x_min_0, y_min_0, x_max_0, y_max_0 = _ee_bounds_to_bounds(
            self.get_info['bounds']
        )
      else:
        x_min_0, y_min_0, x_max_0, y_max_0 = self.crs.area_of_use.bounds
    except AttributeError:
      # `area_of_use` is probable `None`. Parse the geometry from the first
      # image instead (calculated in self.get_info())
      x_min_0, y_min_0, x_max_0, y_max_0 = _ee_bounds_to_bounds(
          self.get_info['bounds']
      )
    # TODO(#40): Investigate data discrepancy (off-by-one) issue.
    x_min, y_min = self.transform(x_min_0, y_min_0)
    x_max, y_max = self.transform(x_max_0, y_max_0)
    self.bounds = x_min, y_min, x_max, y_max

At the end self.bounds is (4454267.556539465, -2212366.2541716346, 4455777.874164572, -2211369.618782846) while we started from (4447800,-2223901,4448801,-2222900)

If I simply hardcode self. bounds to be exactly what I pass, I get

Problem 3. Even so I get back a 2x2 array, not 1x1

Dimensions:                                   (time: 1, X: 2, Y: 2)
Coordinates:
  * time                                      (time) datetime64[ns] 2000-02-24
  * X                                         (X) float32 42.57 42.57
  * Y                                         (Y) float32 -19.99 -20.0

Hello @simonff I thought 2x2 is true as I tried the below native approach and in this approach the result is same. Can you tell us why are you expecting 1x1 if scale is 463m.

import numpy as np

image = ee.Image.pixelCoordinates(ee.Projection(SIN))
scale = 463.3127165279165
x_min, y_min, x_max, y_max = 4447800, -2223901, 4448801, -2222900

x_size = int(np.round((x_max - x_min) / np.abs(scale)))
y_size = int(np.round((y_max - y_min) / np.abs(scale)))
kwargs = {
    "grid": {
        "dimensions": {"width": x_size, "height": y_size},
        "affineTransform": {
            "translateX": x_min,
            "translateY": y_max,
            "scaleX": scale,
            "scaleY": scale * (-1),
        },
        "crsCode": SIN,            
    },
}

params = {
    "expression": image,
    "fileFormat": "NUMPY_NDARRAY",
    **kwargs,
}
result = ee.data.computePixels(params)
result

result:

array([[(4448031.65635826, -2223131.65635826),
        (4448494.96907479, -2223131.65635826)],
       [(4448031.65635826, -2223594.96907479),
        (4448494.96907479, -2223594.96907479)]],
      dtype=[('x', '<f8'), ('y', '<f8')])

P.S. I updated my code with this PR that's why I got the result in meter.

google / Xee

How to get back a single pixel in the MODIS sinusoidal projection? #112