Google Earth Engine – Fix Invalid Geometry Constructed by EE Python API

geometrygoogle-earth-enginepython

I have a problem with Python API of Google Earth Engine.

I'm trying to create a vector mask of a small area based on NDVI. I take Sentinel-2 data as two bands and pass a rectangle (defined as series of hand-picked coordinates) to constrict the ROI. When I try to run the attached code I get an error message:

  File "C:/files/test.py", line 24, in ndROI
    ndviVectors = ndviMask.reduceToVectors(ROIgeometry, 10, 'label') # At this point the "EEException" is thrown
...

  File "C:\Users\glowa\Miniconda2\lib\site-packages\ee\geometry.py", line 83, in __init__
    raise ee_exception.EEException('Invalid GeoJSON geometry.')

EEException: Invalid GeoJSON geometry.

Apparently the constructor returns an invalid polygon, although when printed it looks fine. How can I make sure, that geometries are correct? What I tried and had no success with were:

different CRS and different parameter setup
different coordinate declarations (3 levels of brackets and two)
different places
even replacing Polygon function with Rectangle function (and changing coordinate format accordingly)

all ending up with the same exception

I'm using GEE Python package installed in Miniconda on Windows 10. Below the problematic code.

import ee

ee.Initialize()

# Define function, which produces a NDVI-based vector ROI
def ndROI(imgName,rectROI,diffBand1, diffBand2):

    # Production of rectangle limiting ROI - it's a ee.geometry defined from coordinates and staying in the GEE.
    ROIgeometry = ee.Geometry.Polygon(rectROI, False, 0, True)
    print ROIgeometry
    # Land looks similar to ice - do NDVI from Sentinel 2 to mask away. Get S2 from the latest time, compute ndvi, export. Will do once.
    s1forNDVI = ee.Image(imgName) # Choose nicest scene by hand
    ndvi = s1forNDVI.normalizedDifference([diffBand1, diffBand2]) # select proper bands of sentinel
    ndviMask = ndvi.gte(0.2).updateMask(ndvi.lte(0.2)) # Mask greater-than-0.2, that is plants have 1, water and ice have 0.

    # Vectorise NDVI
    ndviVectors = ndviMask.reduceToVectors(ROIgeometry, 10, 'label') # At this point the "EEException" is thrown
    ndviResultGeometry = ee.FeatureCollection(ndviVectors).geometry(); # convert NDVI product to feature collection and take its geometry
    return ndviResultGeometry

# Define inputs
imgName = 'COPERNICUS/S2/20150817T095016_20160923T153347_T34UDE'
rectROI = [[20.1149, 53.8063], [20.1369, 53.8063], [20.1369, 53.7913], [20.1149, 53.7913], [20.1149, 53.8063]]
diffBand1 = 'B8'
diffBand2 = 'B4'

# Run function
eeNDVIpoly = ndROI(imgName,rectROI,diffBand1, diffBand2)

Best Answer

You have to be careful with *args and **kwargs. Your code runs fine after you fix that. Except when constructing the Polygon, in which case I think there is a bug in the API (I reported it to developers).

import ee

ee.Initialize()

# Define function, which produces a NDVI-based vector ROI
def ndROI(imgName,rectROI,diffBand1, diffBand2):

    # Production of rectangle limiting ROI - it's a ee.geometry defined from coordinates and staying in the GEE.
    ROIgeometry = ee.Geometry.Polygon(rectROI, 
                                      # geodesic=False, 
                                      # maxError=0, 
                                      # evenOdd=True
                                     )
    print(ROIgeometry.getInfo())
    # Land looks similar to ice - do NDVI from Sentinel 2 to mask away. Get S2 from the latest time, compute ndvi, export. Will do once.
    s1forNDVI = ee.Image(imgName) # Choose nicest scene by hand
    ndvi = s1forNDVI.normalizedDifference([diffBand1, diffBand2]) # select proper bands of sentinel
    ndviMask = ndvi.gte(0.2).updateMask(ndvi.lte(0.2)) # Mask greater-than-0.2, that is plants have 1, water and ice have 0.

    # Vectorise NDVI
    ndviVectors = ndviMask.reduceToVectors(geometry=ROIgeometry, 
                                           scale=10, 
                                           labelProperty='label')
    ndviResultGeometry = ee.FeatureCollection(ndviVectors).geometry(); # convert NDVI product to feature collection and take its geometry
    return ndviResultGeometry

# Define inputs
imgName = 'COPERNICUS/S2/20150817T095016_20160923T153347_T34UDE'
rectROI = [[20.1149, 53.8063], [20.1369, 53.8063], [20.1369, 53.7913], [20.1149, 53.7913], [20.1149, 53.8063]]
diffBand1 = 'B8'
diffBand2 = 'B4'

# Run function
eeNDVIpoly = ndROI(imgName,rectROI,diffBand1, diffBand2)

If you switch to python3 you could even visualize it.

UPDATE For those interested in visualizing an interactive map using python as I show in the image, geetools package doesn't hold that ability anymore. Now it is in a new module called ipygee

import ipygee as ui
Map = ui.Map()
Map.show()

Related Solutions

[GIS] Invalid Geometry QGIS

A few options for you:

PostGIS (if available, easy to install locally if not and well worth it)

UPDATE TABLE schema.tablename set geom = ST_MakeValid(geom);

There is also the geometry checker plugin that ships with QGIS but is usually turned off (installed plugins), that has the option to fix all errors.

GRASS v.clean is also an option, available under processing tools in QGIS.

Google Earth Engine – Batch Processing Data of Multiple ROIs in Google Earth Engine

Here's an example that hopefully does something like what you want, but you'll need to modify it to meet your requirements for "manually selected marker points" and Sentinel 2.

In general, do NOT use for-loops or getInfo(). To be clear, DO NOT USE FOR LOOPS unless you have a very good reason to be doing so. If you don't know, use map() for all the reasons described here, here, here and here. (And check the same guides for why to NOT use getInfo() or convert to a list as I've done here.) The reason it's OK to do this here (necessary, in fact) is because Export is a client side function and you can't use a client function in map().

// Get some imagery to play with.
var landsat = ee.ImageCollection("LANDSAT/LC08/C01/T1")
    .filterDate('2016-01-01', '2017-01-01');

var composite = ee.Algorithms.Landsat.simpleComposite({
  collection: landsat,
  asFloat: true
});

var rgbVis = {bands: ["B4", "B3", "B2"], min:0, max: 0.3};
Map.addLayer(composite, rgbVis, "RGB");

// This is a hacky way to get a pixel grid at arbitrary resolution.
var pixels = ee.Image.random().multiply(10000000).toInt32()
    .reduceToVectors({
      reducer: ee.Reducer.countEvery(),
      geometry: Map.getBounds(true),
      geometryType: 'bb' ,
      eightConnected: false,
      scale: 20000,
      crs: 'EPSG:4326'
    });
Map.addLayer(pixels);

// Only do this is you have a few regions.  Not suitable for 
// large feature collections.
var pixelsList = pixels.toList(pixels.size());

// This is one of the few places in the EE API where you need 
// a for-loop and a getInfo() call.  Export is a client function.
for (var i=0; i<pixels.size().getInfo(); i++) {
  Export.image.toDrive({
    image: composite, 
    description: 'foo_' + i, 
    fileNamePrefix: 'foo_' + i, 
    region: ee.Feature(pixelsList.get(i)).geometry(), 
    scale: 30, 
  });
}

If you really need precise control over the export regions, you can make lists of coordinates and turn those into a collection of the ROIs you want to export. All that aside, you still need to click 'Run' on the exports. If you want it completely automated, use the Python API and ee.batch.Export followed by task.start().

Best Answer

Related Solutions

[GIS] Invalid Geometry QGIS

Google Earth Engine – Batch Processing Data of Multiple ROIs in Google Earth Engine

Related Question