Difference in clip results between Earth Engine and Python Gdal using the same geometry

Question

I made scripts that take geometries as input, and outputs NDVI time-series that I can use for classification later. I use ESA's Sentinel-2 Copernicus database via the Python API sentinelsat to get satellite images, then clip them according to the geometries, and calculate the NDVI values for each geometry. Now I am trying to add Google Earth Engine cloud detection layers (for reason's I won't go into here, I need to use ESA's database for the satellite images themselves) so I can filter cloudy pixels. This because Google’s cloud layers are better than ESA's. I manage to download a clip of my geometry from the cloud layer that Google has (currently using their EE package). The issue is however that the result of clipping using the geometry is different from EE than from what I get when I clip ESA's image (the shape and size are different; also if I download the RGB layers from Google and clip them the resulting shape and size are different). I have no clue why exactly that is the case; and I would really love to find out!

When I clip using EE, the resulting image is always is two pixels larger in each dimension (x and y are different but not the channels ofcourse, more on that below) than when I use GDAL to clip ESA's image. It is not just some difference in geometry border handling or buffering it seems, since the shape also slightly differs. I think it is me not understanding how the projections work in EE; but I’ve been trying everything to fix it with no result.

So let’s take the following WKT geometry (in epsg: 32631) as an example: wkt_geom = {'type': 'MultiPolygon', 'coordinates': [[[[614995.9723789722, 5712503.2703056475], [614977.9866070098, 5712516.114248688], [614979.4594190861, 5712518.141167168], [614984.882983979, 5712525.605807577], [614914.0006207882, 5712577.090088011], [614914.1450753526, 5712577.517933527], [614930.0140337138, 5712609.852333707], [614938.7555573103, 5712605.131935651], [614974.024492798, 5712584.449255214], [614977.8123583449, 5712580.957946738], [614980.4038983571, 5712577.924560902], [615003.3580864383, 5712552.329715296], [615019.1074205302, 5712535.583137482], [614996.3756030932, 5712503.881753471], [614995.9723789722, 5712503.2703056475]]]]}

I am going to clip it from ESA’s S2A_MSIL2A_20210316T105031_N0214_R051_T31UFT_20210316T134446 satellite image which is in the same epsg (32631). I’ve requested it from the LTA so it should be available. Clipping this geometry with GDAL warp as such:

wkt_geom = ogr.CreateGeometryFromWkt(geometry.wkt)
drv = ogr.GetDriverByName('ESRI Shapefile')
feature_ds = drv.CreateDataSource("/vsimem/memory_name.shp")
feature_layer = feature_ds.CreateLayer("layer", srs ,geom_type=ogr.wkbPolygon)
featureDefnHeaders = feature_layer.GetLayerDefn()
outFeature = ogr.Feature(featureDefnHeaders)
outFeature.SetGeometryDirectly(wkt_geom)
feature_layer.CreateFeature(outFeature)
feature_ds.FlushCache()
clipped_sat = gdal.Warp('/vsimem/temp.tif', 
          gdal_ds, 
          cutlineDSName="/vsimem/memory_name.shp", 
          cropToCutline=True,
          copyMetadata=True,
          dstNodata = 0)

will result in an image of shape (9, 9, channels) with the following shape (nodata in blue vs data in yellow):

When I however use this geometry to clip in Python with EE as such:

# then I create an EE geom
ee_geom = ee.Geometry.MultiPolygon(
    coords = geom_wkt['coordinates'],
    proj = 'EPSG:32631'
    )
# then I get the relevant image from EE
im = ee.Image('COPERNICUS/S2_SR/20210316T105031_20210316T105442_T31UFT')
# then I clip the image
im = im.clip(ee_geom)
# then I get it in numpy using geemap, though I also tried using EE’s download option which results in the same
ee_im = geemap.ee_to_numpy(ime)

I get an image of shape (11, 11, channels) where again yellow is data and blue is nodata:

The clips are not only always two pixels larger in each dimension, but also are slightly different in shape (as you can see especially well in the corners of the relevant pixels). It might be that its due to differences in projection calculations, or a difference in how GDAL handles the border of a geometry vs EE. Though the documentations do not give me enough help to fix the issue. And I am completely at a loss what causes the differences, or even how to find out why they are different so let alone how I can fix it. Though I am using Python, helping me using Earth Engine itself would be enough, as I can translate it to Python myself. Moreover, for those trying to help me it might be easiest to start in EE and see if you can get it in the shape of my GDAL result since I think its me mishandling EE.

Answer 1

This is just a difference in how pixels that are partially clipped are handled. Earth Engine sets the mask of the partially clipped pixels to a value < 1. If you don't want those pixels, you can mask off any pixels where the mask < 1. That wont change the limits of your region though. If you really want the smaller footprint you can negative buffer by 1/2 a pixel and rely on the partial mask to make up the difference.