2 Extract NICFI Planet data at points (GEDI)

In [1]:

# import earth engine and intialise high volume end-point
import ee
# ee.Authenticate() # you may need to authenticate
ee.Initialize(opt_url='https://earthengine-highvolume.googleapis.com')

# import earth engine and intialise high volume end-point import ee # ee.Authenticate() # you may need to authenticate ee.Initialize(opt_url='https://earthengine-highvolume.googleapis.com')

In [2]:

import os

import csv
import geedim
import geemap

from geeml.utils import eeprint, getCountry, createGrid
from geeml.extract import extractor

import os import csv import geedim import geemap from geeml.utils import eeprint, getCountry, createGrid from geeml.extract import extractor

In [3]:

%load_ext watermark
%watermark
%watermark --iversions

%load_ext watermark %watermark %watermark --iversions

Last updated: 2022-06-29T16:30:55.194517+02:00

Python implementation: CPython
Python version       : 3.9.12
IPython version      : 8.4.0

Compiler    : MSC v.1929 64 bit (AMD64)
OS          : Windows
Release     : 10
Machine     : AMD64
Processor   : Intel64 Family 6 Model 158 Stepping 10, GenuineIntel
CPU cores   : 12
Architecture: 64bit

ee    : 0.2
geemap: 0.13.4
geedim: 1.2.0
csv   : 1.0

Extract data for Sparse imagery (GEDI) and NICFI planet at intersecting points¶

Step 1: Import datasets

Step 2: Prepare data (filter, reduce)

Step 3: Extract data

We convert the GEDI image to points and extract the Planet data at these corresponding points. We use the ESA Worldcover data to remove water, ice and urban areas.¶

Step 1: Import datasets

In [4]:

# Import data
GEDI = ee.ImageCollection("LARSE/GEDI/GEDI02_A_002_MONTHLY")
planet = ee.ImageCollection("projects/planet-nicfi/assets/basemaps/africa")
landcover = ee.ImageCollection("ESA/WorldCover/v100")

# A point in Kenya
poi = ee.Geometry.Point([37.857884,-0.002197])
aoi = getCountry(poi)#kenya

# Import data GEDI = ee.ImageCollection("LARSE/GEDI/GEDI02_A_002_MONTHLY") planet = ee.ImageCollection("projects/planet-nicfi/assets/basemaps/africa") landcover = ee.ImageCollection("ESA/WorldCover/v100") # A point in Kenya poi = ee.Geometry.Point([37.857884,-0.002197]) aoi = getCountry(poi)#kenya

Step 2: Prepare data

In [5]:

# Mask to remove buildings, snow/ice/ and open water
lcmask = landcover.filterBounds(aoi).mosaic().eq([50, 70, 80]).reduce(ee.Reducer.max()).eq(0)

# Filter GEDI data (remove low quality data)
def qualityMask(img):
    return img.updateMask(img.select('quality_flag').eq(1))\
          .updateMask(img.select('degrade_flag').eq(0))\
          .updateMask(lcmask)

# Filter data (to aoi and apply qualityMask and select rh98 band)
dataset = GEDI.filterBounds(aoi).map(qualityMask)\
                .select(['rh98']);

# Set projection and scale
projection = dataset.first().projection()
scale = projection.nominalScale()

mosaic = dataset.mosaic().setDefaultProjection(**{'crs':projection, 'scale':5}).clip(aoi);

# Planet data- get percentiles across all monthly composite planet data
monthlyPlanet = planet.filterBounds(aoi).filter(ee.Filter.eq('cadence','monthly'))\
.reduce(ee.Reducer.percentile([5,25,50,75,95])).clip(aoi)

# Mask to remove buildings, snow/ice/ and open water lcmask = landcover.filterBounds(aoi).mosaic().eq([50, 70, 80]).reduce(ee.Reducer.max()).eq(0) # Filter GEDI data (remove low quality data) def qualityMask(img): return img.updateMask(img.select('quality_flag').eq(1))\ .updateMask(img.select('degrade_flag').eq(0))\ .updateMask(lcmask) # Filter data (to aoi and apply qualityMask and select rh98 band) dataset = GEDI.filterBounds(aoi).map(qualityMask)\ .select(['rh98']); # Set projection and scale projection = dataset.first().projection() scale = projection.nominalScale() mosaic = dataset.mosaic().setDefaultProjection(**{'crs':projection, 'scale':5}).clip(aoi); # Planet data- get percentiles across all monthly composite planet data monthlyPlanet = planet.filterBounds(aoi).filter(ee.Filter.eq('cadence','monthly'))\ .reduce(ee.Reducer.percentile([5,25,50,75,95])).clip(aoi)

spcvGridSize defines a grid, in this case a 30 km by 30 km grid. This is useful for performing Spatial Cross validation¶

In [ ]:

# Download directory
dd = r'D:\Scratch'

# Initialise extractor
gedi = extractor(covariates = monthlyPlanet, target = mosaic, aoi = aoi, scale = 5, dd=dd, spcvGridSize= 30000)

# Extract data in batches of 30 000 points
gedi.extractPoints(gridSize = 50000, batchSize = 30000, filename = 'height.csv')

# Download directory dd = r'D:\Scratch' # Initialise extractor gedi = extractor(covariates = monthlyPlanet, target = mosaic, aoi = aoi, scale = 5, dd=dd, spcvGridSize= 30000) # Extract data in batches of 30 000 points gedi.extractPoints(gridSize = 50000, batchSize = 30000, filename = 'height.csv')

In [ ]: