google-earth-enginegoogle-earth-engine-python-apimatplotlibnumpy
I am using Earth Engine's Python API to convert ee.Image() to Python plots. I am currently using the ee.sampleRectangle() function to convert ee.Image() to a NumPy array and then using matplotlib's imshow() to generate the plot and save them as PNG.
What I have found is that even for the smallest of regions it takes atleast 1 sec for the conversion to NumPy arrays. I have lakhs of images I want to convert. Is there a better way to avoid the bottleneck I am facing?
Matplotlib knows nothing about georeferenced surfaces, it only knows x,y,z coordinates. You can also use Visvis or Mayavi.
you must first extract the x,y, z coordinates of the grid ( a raster is a grid of pixels and with a DEM, the value of the pixel is the elevation, z) with osgeo.gdal. No script here because it is possible to find the solutions on Gis StackExchange or on the web.
after, you can plot the points in 3D
from mpl_toolkits.mplot3d.axes3d import *
import matplotlib.pyplot as plt
from matplotlib import cm
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter3D(x,y,z,c=z,cmap=plt.cm.jet)
plt.show()
and you must reconstruct a 3D grid (surface) with the function griddata of matplotlib (Delaunay)
import numpy as np
from matplotlib.mlab import griddata
# craation of a 2D grid
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)
# interpolation
Z = griddata(x, y, z, xi, yi)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet,linewidth=1, antialiased=True)
plt.show()
The grid (with visvis):
!
The coloured grid (with matplotlib):
but you can also use others interpolation algorithms (Scipy thin-plate spline here) and drape colours
import scipy as sp
import scipy.interpolate
# 2D grid construction
spline = sp.interpolate.Rbf(x,y,z,function='thin-plate')
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)
# interpolation
Z = spline(X,Y)
With Visvis, you can make animations:
You can even plot the 3D contours:
You can improve your i/o time by considering the blocksize of the underlying raster. "Blocksize" refers to the size of the chunks that the raster is written to the hard disk. Typical blocksizes are 128 x 128, or 256 x 256, but it is not uncommon for the blocksize to be N x 1, where N is the number of columns in your raster (or some other non-square format). If this is your case, then considering blocksize will really speed up your reading to memory.
I had a very similar issue; here is some code that I developed that tries to optimize your tile size based on the tiff's blocksize. You want your tile size to be multiples of your blocksize (in x- and y-dimensions). Otherwise, I don't think there's much you can do unless you're willing/able to manipulate the original raster or buy an SSD.
Best Answer
sampleRectangle()is not the ideal way to retrieve image pixels from Earth Engine, because they'll end up converted to JSON and back. Try usingee.Image.getDownloadURL()— this can give you “NumPy binary format” data directly. (I'm not familiar with NumPy to tell you how exactly to use that, but I assume it'll suit your needs.)That said, if your problem is the speed independent of region size, you might find it useful to modify your program to request multiple images in parallel rather than one at a time. (Not too many or your requests will be rejected, but it's entirely normal to do that at all — for example, the Earth Engine Code Editor requests many individual image tiles to make up the interactive map.)