I’m trying to work out how to extract the data along a transect in Google Earth Engine - ideally at defined intervals. Using the Python API in a Jupyter notebook I've mapped some data, buffered a point (to define the region of interest), and mapped a transect. I don't know whether I should use an ee method to extract data constrained along the lineString (which I assume isn't a shape?), or whether I'm heading in the wrong direction and should be exporting the buffered area as a GeoTIFF to process the transect in QGIS/ArcGIS.
import ee
ee.Initialize()
def maskS2clouds(image):
qa = image.select('QA60')
# Bits 10 and 11 are clouds and cirrus, respectively.
cloudBitMask = 1 << 10
cirrusBitMask = 1 << 11
# Both flags should be set to zero, indicating clear conditions.
mask = qa.bitwiseAnd(cloudBitMask).eq(0) \
.And(qa.bitwiseAnd(cirrusBitMask).eq(0))
return image.updateMask(mask).divide(10000)
centre = ee.Geometry.Point(-115.435272, 35.584001,)
centreBuffer = centre.buffer(5000)
dataset = ee.ImageCollection('COPERNICUS/S2_SR') \
.filterDate('2020-07-01', '2020-07-31') \
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20)) \
.map(maskS2clouds)
visualization = {
'min': 0.0,
'max': 0.3,
'bands': ['B4', 'B3', 'B2'], # Red, Green, BLue. 10m resolution.
}
Map.setCenter(-115.435272, 35.584001, 12)
# Map = geemap.Map(center=[35.584001, -115.435272], zoom=14)
Map.addLayer(dataset.mean(), visualization, 'RGB')
Map.addLayer(centre,
{'color': 'black'},
'Geometry [black]: point');
Map.addLayer(centreBuffer,
{'color': 'red'},
'Result [red]: point.buffer')
transect = ee.Geometry.LineString([[-115.4, 35.584001], [-115.45, 35.584001]]);
Map.addLayer(transect, {'color': 'Green'}, 'transect');
Map
You may find these resources by Qiusheng Wu helpful.
His Common module API has the function extract_transect(image, line, reducer='mean', n_segments=100, dist_interval=None, scale=None, crs=None, crsTransform=None, tileScale=1.0, to_pandas=False, **kwargs)
This extracts a transect from an image.
From the third link above, the Jupyter notebook is:
import ee
import geemap
from bqplot import pyplot as plt
Map = geemap.Map()
Map.add_basemap("TERRAIN")
Map
image = ee.Image('USGS/SRTMGL1_003')
vis_params = {
'min': 0,
'max': 4000,
'palette': ['006633', 'E5FFCC', '662A00', 'D8D8D8', 'F5F5F5']}
Map.addLayer(image, vis_params, 'SRTM DEM', True, 0.5)
# Use the drawing tool to draw any line on the map.
line = Map.user_roi
if line is None:
line = ee.Geometry.LineString([[-120.223279, 36.314849], [-118.926969, 36.712192], [-117.202217, 36.756215]])
Map.addLayer(line, {}, "ROI")
Map.centerObject(line)
line.getInfo()
reducer='mean' # Any ee.Reducer, e.g., mean, median, min, max, stdDev
transect = geemap.extract_transect(image, line, n_segments=100, reducer=reducer, to_pandas=True)
transect
fig = plt.figure()
plt.plot(transect['distance'], transect[reducer])
plt.xlabel('Distance')
plt.ylabel("Elevation")
plt.show()
A JavaScript version is:
/***
* Reduces image values along the given line string geometry using given reducer.
*
* Samples image values using image native scale, or opt_scale
*/
function reduceImageProfile(image, line, reducer, scale, crs) {
var length = line.length();
var distances = ee.List.sequence(0, length, scale)
var lines = line.cutLines(distances, ee.Number(scale).divide(5)).geometries();
lines = lines.zip(distances).map(function(l) {
l = ee.List(l)
var geom = ee.Geometry(l.get(0))
var distance = ee.Number(l.get(1))
geom = ee.Geometry.LineString(geom.coordinates())
return ee.Feature(geom, {distance: distance})
})
lines = ee.FeatureCollection(lines)
// reduce image for every segment
var values = image.reduceRegions( {
collection: ee.FeatureCollection(lines),
reducer: reducer,
scale: scale,
crs: crs
})
return values
}
// Define a line across the Olympic Peninsula, USA.
// Import a digital surface model and add latitude and longitude bands.
var elevImg = ee.Image('JAXA/ALOS/AW3D30/V2_2').select('AVE_DSM');
var profile = reduceImageProfile(elevImg, transect, ee.Reducer.mean(), 100)
print(ui.Chart.feature.byFeature(profile, 'distance', ['mean']))