Getting DEM Values as NumPY Array in Earth Engine

Question

I'm trying to convert Earth Engine data to a local GeoTIFF manually, using more or less the method found here. One crucial step is converting the DEM raster values to a NumPy array. I've been trying to do so as follows:

ee.Initialize()

# define ROI
area = ee.Geometry.Rectangle(bounds)

# set image
img = ee.Image('USGS/SRTMGL1_003')

# get the lat lon and add the elevation data back in
latlng = ee.Image.pixelLonLat().addBands(img)

# reduce to list
latlng = latlng.reduceRegion(reducer=ee.Reducer.toList(), geometry=area, maxPixels=1e8, scale=20)

# put the values into lists
elev_values = np.array((ee.Array(latlng.get("elevation")).getInfo()))
lats = np.array((ee.Array(latlng.get("latitude")).getInfo()))
lngs = np.array((ee.Array(latlng.get("longitude")).getInfo()))

While I have had success converting unique latitude and longitude values to a NumPy array, I keep ending up with an empty array of elevation values. The current code throws the following error:

Traceback (most recent call last):
  File "test_one.py", line 392, in <module>
    gDem = GoogleDEM([40.01, -21.01, 40.02, -21.00])
  File "test_one.py", line 41, in __init__
    elev_values = np.array((ee.Array(latlng.get("elevation")).getInfo()))
  File "C:\Users\lmonn\AppData\Local\Programs\Python\Python37-32\lib\site-packages\ee\computedobject.py", line 95, in getInfo
    return data.computeValue(self)
  File "C:\Users\lmonn\AppData\Local\Programs\Python\Python37-32\lib\site-packages\ee\data.py", line 490, in computeValue
    return send_('/value', ({'json': obj.serialize(), 'json_format': 'v2'}))
  File "C:\Users\lmonn\AppData\Local\Programs\Python\Python37-32\lib\site-packages\ee\data.py", line 1186, in send_
    raise ee_exception.EEException(json_content['error']['message'])
ee.ee_exception.EEException: Array: No numbers in 'values', must provide a type.

My guess would be that there is some nesting of objects taking place that I don't quite have a grasp on.

How can I fix this issue? Is there a simpler way to get a NumPy Array of the elevation values for a DEM on Earth Engine?

I can't tell from your question how GoogleDEM() is defined, but have you checked that you are correctly ordering the x/y values? — Tyler Erickson
– Tyler Erickson, Commented May 3, 2019 at 23:13
@TylerErickson Thanks for the response! The bounds are in the order provided by the ee.Geometry.Rectangle documentation. GoogleDEM() is actually a class I'm working on. The code above is part of the implementation of its constructor. — LIAM M
– LIAM M, Commented May 4, 2019 at 0:11

Guido Lemoine · Accepted Answer · 2019-05-05 10:45:30Z

This can be done directly, without the need to go through numpy array construction via getInfo(). getInfo() is limited to 5000 records, so is not a good solution for larger selections. Check out ee.Image.getDownloadURL().

import ee
import requests
import zipfile

import numpy as np

ee.Initialize()
bounds = [-97.94, 26.81, -96.52, 26.84] ## sample land / sea bounds
area = ee.Geometry.Rectangle(-97.94, 26.81, -96.52, 26.84)
img = ee.Image('USGS/SRTMGL1_003').clip(area)

url = img.getDownloadURL({'name': 'MyPieceOfDEM', 'crs': 'EPSG:4326', 'scale': 90})

print(url)

filename = 'MyDEM.zip'

# Download the subset
r = requests.get(url, stream=True)
with open(filename, 'wb') as fd:
    for chunk in r.iter_content(chunk_size=1024):
        fd.write(chunk)

# Extract the GeoTIFF for the zipped download
z = zipfile.ZipFile(filename)
z.extractall()

This will produce the following file, which you can inspect with gdalinfo

gdalinfo MyPieceOfDEM.elevation.tif
Driver: GTiff/GeoTIFF
Files: MyPieceOfDEM.elevation.tif
Size is 1758, 41
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-97.940530650170658,26.842469173247011)
Pixel Size = (0.000808483755708,-0.000808483755708)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -97.9405307,  26.8424692) ( 97d56'25.91"W, 26d50'32.89"N)
Lower Left  ( -97.9405307,  26.8093213) ( 97d56'25.91"W, 26d48'33.56"N)
Upper Right ( -96.5192162,  26.8424692) ( 96d31' 9.18"W, 26d50'32.89"N)
Lower Right ( -96.5192162,  26.8093213) ( 96d31' 9.18"W, 26d48'33.56"N)
Center      ( -97.2298734,  26.8258953) ( 97d13'47.54"W, 26d49'33.22"N)
Band 1 Block=1758x2 Type=Int16, ColorInterp=Gray

Andrew Cutts · Accepted Answer · 2019-05-05 08:55:37Z

As you correctly point out the issue is to do with values over the ocean return no values. I've just added a try and except statement to your code to catch if this is the case to retun a Numpy array containing NaN's of the same dimensions as the lats array (though could be lngs if you prefer). I have also tried it on a bounds containing both land and sea and it worked for me.

import ee
import numpy as np
ee.Initialize()
bounds = [-97.94, 26.81, -96.52, 26.84] ## sample land / sea bounds
area = ee.Geometry.Rectangle(bounds)
img = ee.Image('USGS/SRTMGL1_003')
latlng = ee.Image.pixelLonLat().addBands(img)
latlng = latlng.reduceRegion(reducer=ee.Reducer.toList(), geometry=area, maxPixels=1e8, scale=20)
lats = np.array((ee.Array(latlng.get("latitude")).getInfo()))
lngs = np.array((ee.Array(latlng.get("longitude")).getInfo()))
try:
    elev_values = np.array((ee.Array(latlng.get("elevation")).getInfo()))
except:
    elev_values = np.full_like(lats, np.nan,dtype=np.float64)
print(list(elev_values)) ## print as list to check

This is a really nice lightweight solution. The only thing I would note is that if you make larger requests you'll run into issues with getInfo(). — LIAM M
– LIAM M, Commented May 5, 2019 at 19:08

LIAM M · Accepted Answer · 2019-05-04 18:22:04Z

0

The issue has not to do with negative latitude and longitude values as I previously suggested. But, rather, is an issue of selected values over the ocean. I am going to use this error to help check appropriate bounds. However, if someone else encounters this problem and does not want to deal with throwing an error, the void-filled 3 arc-sec dataset might be a solution.

answered May 4, 2019 at 18:22

LIAM M

634 bronze badges

Add a comment |

Stack Exchange Network

Getting DEM Values as NumPY Array in Earth Engine

3 Answers 3

Your Answer

Hot Network Questions

Getting DEM Values as NumPY Array in Earth Engine

3 Answers 3

Your Answer

Sign up or log in

Post as a guest

Related

Hot Network Questions