0

I try to extract for a collection of 1 ha square polygons the bands coefficients from a preprocessed Landsat CCDC Google Earth Engine Asset ee.ImageCollection('projects/CCDC/v3').

My plan is the following:

  1. Try for a given ee.Geometry.Point() location: OK.
  2. Try for a 1 ha square polygon using ee.Image.sampleRectangle(): FAIL, returns the same results as for the point Geometry while the aoi contains several pixels (asset resolution = 30 m).
  3. Loop through the ee.FeatureCollection(): need input on an efficient way to iterate. I was thinking of using map() and then concatenate/merge the results into an xarray.

Notes:

  • I am aware of the gee-ccdc-tools package, but I am looking for a Python solution.
  • I wonder if this could be a projection issue as hinted here.
import ee
import geemap
import numpy as np

ee.Initialize()

poi = ee.Geometry.Point(-53.423000,-13.623000) # Point of interest
aoi = poi.buffer(50).bounds() # 1 ha rectangle centered on poi

ccdc = ee.ImageCollection('projects/CCDC/v3') # Load asset
ccdc_mosaic = ccdc.mosaic() # Create mosaic from 990 non-overlapping images covering the world

Map = geemap.Map()
Map.addLayer(ccdc_mosaic, {'bands':'GREEN_magnitude', 'min': 0, 'max': 0.01}, 'ccdc_mosaic')
Map.addLayer(poi, {'color': 'red'}, 'poi')
Map.addLayer(aoi, {'color': 'blue'}, 'aoi')
Map.centerObject(aoi)
Map # Display for inspection with 'Inspector'

tEnd_poi = np.array(ccdc_mosaic.sampleRectangle(poi).get('tEnd').getInfo())
tEnd_poi

# This pixel contains 4 time segments: OK
array([[[1999.78662184, 2007.64972982, 2013.694921  , 2019.67435735]]])


GREEN_coefs_poi = np.array(ccdc_mosaic.sampleRectangle(poi).get('GREEN_coefs').getInfo())
GREEN_coefs_poi

# The 8 coefficients for the GREEN band for the 4 time segments: OK
array([[[[ 2.31620175e+01, -1.15526059e-02,  2.27645180e-02,
           0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
           0.00000000e+00,  0.00000000e+00],
         [ 3.99429518e-02,  0.00000000e+00,  1.79411590e-03,
          -1.11117008e-03,  0.00000000e+00,  0.00000000e+00,
           0.00000000e+00,  0.00000000e+00],
         [ 1.17203301e-01, -3.64512731e-05,  4.46052439e-03,
          -1.55708945e-03,  0.00000000e+00,  0.00000000e+00,
           0.00000000e+00,  0.00000000e+00],
         [ 4.13283961e-02,  0.00000000e+00,  2.13285985e-03,
           0.00000000e+00, -1.16010972e-03,  0.00000000e+00,
           0.00000000e+00,  0.00000000e+00]]]])

# Get coefficients for the area of interest
GREEN_coefs_aoi = np.array(ccdc_mosaic.sampleRectangle(aoi).get('GREEN_coefs').getInfo())

GREEN_coefs_poi == GREEN_coefs_aoi
# True...
1
  • you should not use sampleRectangle but export the image either as an asset or a .Tif. be aware that if you export to tif you will have an issue as exported bands cannot be arrays. you can have a look here: gis.stackexchange.com/questions/393671/… Commented Dec 9, 2022 at 4:23

1 Answer 1

1

You can turn the CCDC segments array image into an image with a band per array element, then use getDownloadURL() to extract the data. You will end up with a lot of bands, and you cannot download more than 32MB. The below script of 1000x1000m and 6 segments gives you 250 bands and a 17MB download.

Have a look at Noel's article for some more tips on downloads.

import ee

ee.Authenticate() # Skip if you're running this in SEPAL
ee.Initialize()


def flattenSegments(segments):
    maxSegments = (segments
        .select('tStart')
        .arrayLength(0)
        .reduceRegion(
            reducer=ee.Reducer.max(),
            geometry=geometry,
            scale=30,
            maxPixels=1e9
        ).getNumber('tStart'))

    return ee.Image(
        ee.List.sequence(0, maxSegments.subtract(1))
        .iterate(
            lambda segmentIndex, acc:
                ee.Image(acc).addBands(toSegment(segmentIndex, segments)),
            ee.Image().select([])
        )
    )


def toSegment(segmentIndex, segments):
    segmentIndex = ee.Number(segmentIndex).int8()
    segment = segments.arraySlice(0, segmentIndex, segmentIndex.add(1))
    coefs = ee.Image(
        ee.List.sequence(0, 7)
        .iterate(
            lambda coefIndex, acc:
                ee.Image(acc).addBands(toCoef(coefIndex, segment)),
            ee.Image().select([])
        )
    )
    otherBands = (segment
        .select(
            segment
            .bandNames()
            .filter(
                ee.Filter.stringEndsWith('item', '_coefs').Not()
            )
        )
        .updateMask(coefs.select(0).mask())
        .arrayGet([0]))
    return (otherBands
        .addBands(coefs)
        .regexpRename('(.*)', segmentIndex.format('%d').cat('_$1'), False))


def toCoef(coefIndex, segment):
    coefIndex = ee.Number(coefIndex).int8()
    coef = (segment
        .select('.*_coefs')
        .arraySlice(1, coefIndex, coefIndex.add(1))
        .regexpRename(
            '(.*)_coefs',
            ee.String('$1_coef_').cat(coefIndex.format('%d'))
        ))
    return (coef
        .updateMask(coef.select(0).arrayLength(0))
        .arrayGet([0, 0]))


geometry = (ee.Geometry.Point([-64.44248584974648, -7.473629193079254])
    .buffer(1000).bounds())

segments = (ee.ImageCollection('projects/CCDC/v3')
    .filterBounds(geometry)
    .mosaic())

flattened = flattenSegments(segments)
url = flattened.getDownloadURL({
    'format': "NPY",
    'region':  geometry
})

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.