Issue when transforming an EE.Image into numpy array

Question

I use Google Earth Engine with Python 3 with Colab's notebooks.

I would like to transform the Google Images into Numpy arrays to be used for further processing.

I am using the toList Reducer to have a list of values corresponding to B8 band of a Sentinel2 tiles. This is a one dimension list so I need to reshape it properly. However, I have a problem when doing a reshape on a large tile. Indeed, the output image contains a visual shift when the scale is too low (100m). When I do my tests with a much larger scale (1000m), this offset no longer exists. Obviously I would like to have this result with much smaller scales.

There is the ouput for the 100m processing (KO) :

And there is for the 1000m (OK) :

You can find an example of the code here to test it for yourself (with scale of 1000 and scale of 100) : https://colab.research.google.com/drive/1YFWxy3JOWiJ9UAaWsCVRvOwxx_0mobsd

import ee
import numpy as np

ee.Authenticate()
ee.Initialize()
 ############### Function ######################
 #---- Compute the transformation Earth Engine to numpy
def prepareData(collection, testingScale):
    CRS=collection.first().projection().getInfo()['crs'] # Get CRS
    img=collection.median().float()  # Compute an image and cast it into float
    countRes=img.reduceRegion(reducer=ee.Reducer.count(),\ # Get the size of the image
                 crs=CRS,\
                 geometry=collection.geometry().bounds(),\
                 maxPixels=1e13,\
                 scale=testingScale) 
    img = img.reduceRegion(reducer=ee.Reducer.toList(),\ # Transform earth engine Image into a one dimension array
                 crs=CRS,\
                 geometry=collection.geometry().bounds(),\
                 maxPixels=1e13,\
                 scale=testingScale)
    dataB8 = np.array((ee.Array(img.get("B8")).getInfo()))
    dim=np.sqrt(countRes.get('B8').getInfo()).astype(np.int16) # Get the dimension for the reshape into two dimension array
    nrows, ncols = dim, dim 
    image = dataB8.reshape((ncols,nrows))
    return image

############### Main body ######################
#---- Get an example tile of Sentinel2
point=ee.Geometry.Point([-70.066310,-21.341900]) 
collection = ee.ImageCollection("COPERNICUS/S2").filterBounds(point)\
                                  .filterDate("2018-01-01","2019-01-10")\
                                      .filterMetadata("CLOUDY_PIXEL_PERCENTAGE","less_than",10)\
                                  .select(['B8'])
#---- Prepare a numpy image
image = prepareData(collection, 1000) # Reprojection with scale 1000 OK !     Scale 100 non ok !
display(image)

#----- Plot the numpy image
import matplotlib.pyplot as plt
plt.imshow(image)
plt.show()

Has anyone ever had this kind of problem or had any ideas for explanation or improvement?

Answer 1

There is either a mismatch between the dimensions you are trying to reshape the array to and the number of pixels being transferred from Earth Engine servers to the Colab client or that the pixel order is not correct. There are two problems with this method:

1. getInfo() is not going to be reliable for transferring an 1,205,604 element list server to client. This is anecdotally evident in the long time is takes to complete the transfer (timeout concerns) and that I get a different looking result than you do at the 100m scale.

2. You are assuming that the region is a square (np.sqrt()) to get the dimensions. It works in this case, but in general, that is huge assumption.

A better method is to use ee.Image.sampleRectangle(), which returns a 2D array of pixels - no need to reshape. Please see this post for an example.