Google Earth Engine integration

Author: jonas-eberle

Readme.md

@@ -1,7 +1,13 @@
 pyEOM
 =====
 
-Earth Observation Monitor
+pyEOM is the Python package for the Earth Observation Monitor (www.earth-observation-monitor.net), a project developed by the Friedrich-Schiller-University Jena, Department for Earth Observation (www.eo.uni-jena.de).
+
+You can use pyEOM to extract time-series data for different MODIS products and Landsat data and a given geometry (individual Pixel or Polygon).
+
+*Please note:* The development of this package is in an early state. Documentation, testing and completion of the work will follow.
+
+If you have any questions, please do not hesitate to contact me: Jonas Eberle <jonas.eberle@uni-jena.de>
 
 Installation
 -----------
@@ -18,6 +24,8 @@ Installation
 * gdalbuildvrt
 * gdal_merge.py
 
+Please make sure that the command-line tools are within the System path. Otherwise you have to specify the path to the tools within the GDAL class (class variable "path") in processing.py.
+
 *Further requirements:*
 For MODIS HDF data processing, HDF4 image format has to be available for gdal command-line tools.
 
@@ -38,3 +46,32 @@ Usage
         'source': 'LPDAAC'
     })
     output = ingest.start()
+
+Installation
+------------
+
+No setup.py is available at the moment. Please download the pyEOM folder and set the Python path to the parent directory.
+    import sys
+    sys.path.append('/directory/to/parent/folder/of/pyEOM')
+
+Available datasets
+------------------
+* MODIS Surface Reflectance: MOD09A1, MOD09Q1
+* MODIS Snow: MOD10A1, MYD10A1, MOD10A2, MYD10A2, MOD10C1, MYD10C1, MOD10C2, MYD10C2, MOD10CM, MYD10CM
+* MODIS Land Surface Temperature: MOD11A1, MYD11A1, MOD11A2, MYD11A2, MOD11B1, MYD11B1, MOD11C1, MYD11C1, MOD11C2, MYD11C2, MOD11C3, MYD11C3
+* MODIS Vegetation Indices: MOD13A1, MYD13A1, MOD13A2, MYD13A2, MOD13A3, MYD13A3, MOD13C1, MYD13C1, MOD13C2, MYD13C2, MOD13Q1, MYD13Q1
+* MODIS MCD15A2, MCD15A3
+* MODIS MCD43A3
+* MODIS MOD44B
+* MODIS MCD45A1
+
+Further information to the datasets can be found at the following sites:
+* http://lpdaac.usgs.gov/products/modis_products_table
+* NSIDC
+
+Google Earth Engine
+-------------------
+
+If you have an account for the Google Earth Engine (https://earthengine.google.org), you can also use the Python bindings to have access to specific MODIS products and Landsat datasets.
+
+Please set the variables MY_SERVICE_ACCOUNT and MY_PRIVATE_KEY_FILE in gee_init.py to your service account username and the private key file.

pyEOM/datasets/Landsat.py

@@ -1 +1,44 @@
 __author__ = 'Jonas Eberle <jonas.eberle@eberle-mail.de>'
+from pyEOM.sources import *
+
+class GEE(GEESource):
+    dataset = None
+    start = '1900-01-01'
+    end = '2100-01-01'
+    epsg = 'EPSG:4326'
+    scale = 30
+
+    def __init__(self, task):
+        self.task = task
+        if 'dataset' not in task:
+            raise Exception('No dataset found in task')
+            sys.exit(1)
+        self.dataset = task['dataset']
+        if 'geom' not in task:
+            raise Exception('No geometry found in task')
+            sys.exit(1)
+        if 'start' in task:
+            self.start = task['start']
+        if 'end' in task:
+            self.end = task['end']
+        if 'scale' in task:
+            self.scale = int(task['scale'])
+        if 'publishPath' not in task:
+            self.publishPath = os.getcwd()
+        else:
+            self.publishPath = task['publishPath']
+
+        self.connect()
+
+        self.geom = self.processGeometry(task['geom'])
+
+    def setPublishPath(self, path):
+        self.publishPath = path
+
+    def run(self):
+        #                   dataset     bands   geom        start       end     epsg        scale       crsTransform
+        data = self.getData(self.dataset, None, self.geom, self.start, self.end, self.epsg, self.scale, None)
+        return data
+
+    def ingest(self, format=None):
+        return self.run()

pyEOM/datasets/MODIS.py

@@ -35,17 +35,43 @@ def getTilesFromGeometry(self, geom):
 
 
 class GEE(GEESource):
+    dataset = None
+    start = '1900-01-01'
+    end = '2100-01-01'
+    epsg = 'SR-ORG:6974'
+
+    def __init__(self, task):
+        self.task = task
+        if 'dataset' not in task:
+            raise Exception('No dataset found in task')
+            sys.exit(1)
+        self.dataset = task['dataset']
+        if 'geom' not in task:
+            raise Exception('No geometry found in task')
+            sys.exit(1)
+        if 'start' in task:
+            self.start = start
+        if 'end' in task:
+            self.end = end
+        if 'publishPath' not in task:
+            self.publishPath = os.getcwd()
+        else:
+            self.publishPath = task['publishPath']
 
-    def __init__(self, task, dataset):
         self.connect()
-        pass
 
-    def process(self):
-        pass
+        self.geom = self.processGeometry(task['geom'])
+
+    def setPublishPath(self, path):
+        self.publishPath = path
 
-    def ingest(self, downloadPath=None):
+    def run(self):
+        #                   dataset     bands   geom        start       end     epsg        scale   crsTransform
+        data = self.getData(self.dataset, None, self.geom, self.start, self.end, self.epsg, None, True)
+        return data
 
-        pass
+    def ingest(self, format=None):
+        return self.run()
 
 
 def TestDatasetLPDAAC():

pyEOM/processing.py

@@ -317,6 +317,7 @@ def compress(self, type):
 
 class GDAL(object):
     outputfilext = {'GTiff': '.tif', 'netCDF': '.nc', 'HDF4Image': '.hdf','VRT': '.vrt', 'KEA': '.kea', }
+    path = '' # with a training slash
 
     resample_methods = ['near', 'bilinear', 'cubic', 'cubicspline', 'lanczos']
 
@@ -346,7 +347,7 @@ def execute(self, cmd, output):
     def gdal_translate(self, input, output, format="GTiff"):
         #if not os.path.isfile(input):
         #    raise Exception('Input file not found!')
-        return self.execute(['gdal_translate', '-of', format, input, output], output)
+        return self.execute([self.path+'gdal_translate', '-of', format, input, output], output)
 
 
     def gdal_merge(self, inputfiles, output, format="GTiff", nodata=""):
@@ -355,7 +356,7 @@ def gdal_merge(self, inputfiles, output, format="GTiff", nodata=""):
         #if len(inputfiles) <= 1:
         #    raise Exception('Inputfiles needs more than 1 file!')
 
-        return self.execute(['gdal_merge.py', '-o', output, '-of', format, inputfiles], output)
+        return self.execute([self.path+'gdal_merge.py', '-o', output, '-of', format, inputfiles], output)
 
 
     def gdalwarp(self, inputfile, output, format, dstnodata, epsg, resample, cutline):
@@ -378,14 +379,14 @@ def gdalwarp(self, inputfile, output, format, dstnodata, epsg, resample, cutline
         else:
             format = ''
 
-        return self.execute(['gdalwarp', format, srcnodata, dstnodata, reproject, cutline, inputfile, output], output)
+        return self.execute([self.path+'gdalwarp', format, srcnodata, dstnodata, reproject, cutline, inputfile, output], output)
 
 
     def gdal_compress(self, inputfile, output, type):
-        return self.execute(['gdal_translate', '-co COMPRESS='+type, inputfile, output], output)
+        return self.execute([self.path+'gdal_translate', '-co COMPRESS='+type, inputfile, output], output)
 
     def gdallocationinfo(self, inputfile, x, y):
-        return self.execute(['gdallocationinfo', '-valonly', '-geoloc', '-wgs84', inputfile, str(x), str(y)], None)
+        return self.execute([self.path+'gdallocationinfo', '-valonly', '-geoloc', '-wgs84', inputfile, str(x), str(y)], None)
 
     def gdalbuildvrt(self, inputs, output, separate=True):
         inputfilelist = ''
@@ -397,4 +398,4 @@ def gdalbuildvrt(self, inputs, output, separate=True):
         else:
             separate = ''
 
-        return self.execute(['gdalbuildvrt', separate, inputfilelist, output, inputs], output)
+        return self.execute([self.path+'gdalbuildvrt', separate, inputfilelist, output, inputs], output)

pyEOM/sources.py

@@ -3,6 +3,7 @@
 import pyEOM.gee_init as gee
 import ee
 import sys
+import ogr
 
 def TestGEESource():
     source = GEESource()
@@ -28,6 +29,16 @@ def connect(self, username=None, password=None):
         else:
             raise Exception('No valid user authentication found')
 
+    def processGeometry(self, wkt):
+        if 'POINT' in wkt:
+            geom = ogr.CreateGeometryFromWkt(wkt)
+            return self.getEEPoint(geom.GetX(), geom.GetY())
+        elif 'POLYGON' in wkt:
+            return self.getEEPolygon(wkt)
+        else:
+            raise Exception('Wrong geometry type (either POINT or POLYGON necessary)')
+            sys.exit(1)
+
     def getEEPoint(self, x, y):
         return ee.Geometry.Point(x, y)
 
@@ -45,23 +56,30 @@ def getData(self, dataset, bands, eeGeometry, start, end, epsg='SR-ORG:6974', sc
             if isinstance(bands, basisstring):
                 bands = [bands]
             collection = collection.select(bands)
+        if crsTransform == True:
+            image = collection.first().getInfo()
+            crsTransform = image['bands'][0]['crs_transform']
         try:
             data = collection.getRegion(eeGeometry, scale, epsg, crsTransform).getInfo()
         except Exception, e:
-            if e.message['code'] == 400 and "Too many values" in e.message['message']:
-                LOGGER.info('Too many pixels in polygon: '+e.message['message'])
-
-                msgAr = e.message['message'].split(' ')
-                points = int(msgAr[msgAr.index('points')-1])
-                bands = int(msgAr[msgAr.index('bands')-1])
-                images = int(msgAr[msgAr.index('images')-1])
-                limit = 1048576
-            else:
-                raise e
-                sys.exit(1)
+            raise e
+            #if e.message['code'] == 400 and "Too many values" in e.message['message']:
+            #    LOGGER.info('Too many pixels in polygon: '+e.message['message'])
+            #
+            #    msgAr = e.message['message'].split(' ')
+            #    points = int(msgAr[msgAr.index('points')-1])
+            #    bands = int(msgAr[msgAr.index('bands')-1])
+            #    images = int(msgAr[msgAr.index('images')-1])
+            #    limit = 1048576
+            #else:
+            #    raise e
+            #    sys.exit(1)
 
         return data
 
+    def reduceRegion(self):
+        pass
+
     def download(self, image, path):
         urls = []
         if isinstance(image, ee.Image):

pyEOM/tasks.py

@@ -142,8 +142,9 @@ def start(self):
             return self.processing['source'].ingest(self.format)
 
         elif "LANDSAT" == dataset[0]:
-            pass
-            #self.processing[dataset] = Landsat.Dataset(dataset)
+            from datasets import Landsat
+            self.processing['source'] = Landsat.GEE(task)
+            return self.processing['source'].run()
         elif "NOAA" == dataset[0]:
             pass