Add COG and STAC support

giswqs · giswqs · commit 6fcaa195b64b · 2020-12-06T12:15:35.000-05:00
diff --git a/eefolium/common.py b/eefolium/common.py
@@ -2332,6 +2332,299 @@ def load_GeoTIFFs(URLs):
     return ee.ImageCollection(collection)
 
 
+def get_COG_tile(url, titiler_endpoint="https://api.cogeo.xyz/", **kwargs):
+    """Get a tile layer from a Cloud Optimized GeoTIFF (COG).
+        Source code adapted from https://developmentseed.org/titiler/examples/Working_with_CloudOptimizedGeoTIFF_simple/
+
+    Args:
+        url (str): HTTP URL to a COG, e.g., https://opendata.digitalglobe.com/events/mauritius-oil-spill/post-event/2020-08-12/105001001F1B5B00/105001001F1B5B00.tif
+        titiler_endpoint (str, optional): Titiler endpoint. Defaults to "https://api.cogeo.xyz/".
+
+    Returns:
+        tuple: Returns the COG Tile layer URL and bounds. 
+    """
+    import requests
+
+    params = {"url": url}
+
+    TileMatrixSetId='WebMercatorQuad'
+    if "TileMatrixSetId" in kwargs.keys():
+        TileMatrixSetId = kwargs["TileMatrixSetId"]
+    if "tile_format" in kwargs.keys():
+        params["tile_format"] = kwargs["tile_format"]
+    if "tile_scale" in kwargs.keys():
+        params["tile_scale"] = kwargs["tile_scale"]
+    if "minzoom" in kwargs.keys():
+        params["minzoom"] = kwargs["minzoom"]
+    if "maxzoom" in kwargs.keys():
+        params["maxzoom"] = kwargs["maxzoom"]
+
+    r = requests.get(
+        f"{titiler_endpoint}/cog/{TileMatrixSetId}/tilejson.json",
+        params = params
+    ).json()
+
+    return r["tiles"][0]
+
+
+def get_COG_bounds(url, titiler_endpoint="https://api.cogeo.xyz/"):
+    """Get the bounding box of a Cloud Optimized GeoTIFF (COG).
+
+    Args:
+        url (str): HTTP URL to a COG, e.g., https://opendata.digitalglobe.com/events/mauritius-oil-spill/post-event/2020-08-12/105001001F1B5B00/105001001F1B5B00.tif
+        titiler_endpoint (str, optional): Titiler endpoint. Defaults to "https://api.cogeo.xyz/".
+
+    Returns:
+        list: A list of values representing [left, bottom, right, top]
+    """    
+    import requests
+
+    r = requests.get(
+        f"{titiler_endpoint}/cog/bounds",
+        params = {
+            "url": url
+        }
+    ).json()   
+
+    bounds = r["bounds"] 
+    return bounds
+
+
+def get_COG_center(url, titiler_endpoint="https://api.cogeo.xyz/"):
+    """Get the centroid of a Cloud Optimized GeoTIFF (COG).
+
+    Args:
+        url (str): HTTP URL to a COG, e.g., https://opendata.digitalglobe.com/events/mauritius-oil-spill/post-event/2020-08-12/105001001F1B5B00/105001001F1B5B00.tif
+        titiler_endpoint (str, optional): Titiler endpoint. Defaults to "https://api.cogeo.xyz/".
+
+    Returns:
+        tuple: A tuple representing (longitude, latitude)
+    """    
+    bounds = get_COG_bounds(url, titiler_endpoint)
+    center=((bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2)  # (lat, lon)
+    return center
+
+
+def get_COG_bands(url, titiler_endpoint="https://api.cogeo.xyz/"):
+    """Get band names of a Cloud Optimized GeoTIFF (COG).
+
+    Args:
+        url (str): HTTP URL to a COG, e.g., https://opendata.digitalglobe.com/events/mauritius-oil-spill/post-event/2020-08-12/105001001F1B5B00/105001001F1B5B00.tif
+        titiler_endpoint (str, optional): Titiler endpoint. Defaults to "https://api.cogeo.xyz/".
+
+    Returns:
+        list: A list of band names
+    """    
+    import requests
+
+    r = requests.get(
+        f"{titiler_endpoint}/cog/info",
+        params = {
+            "url": url,
+        }
+    ).json()
+
+    bands = [b[1] for b in r['band_descriptions']]
+    return bands
+ 
+
+def get_STAC_tile(url, bands=None, titiler_endpoint="https://api.cogeo.xyz/", **kwargs):
+    """Get a tile layer from a single SpatialTemporal Asset Catalog (STAC) item.
+
+    Args:
+        url (str): HTTP URL to a STAC item, e.g., https://canada-spot-ortho.s3.amazonaws.com/canada_spot_orthoimages/canada_spot5_orthoimages/S5_2007/S5_11055_6057_20070622/S5_11055_6057_20070622.json
+        titiler_endpoint (str, optional): Titiler endpoint. Defaults to "https://api.cogeo.xyz/".
+
+    Returns:
+        tuple: Returns the COG Tile layer URL and bounds. 
+    """
+    import requests
+
+    params = {"url": url}
+
+    TileMatrixSetId='WebMercatorQuad'
+    if "TileMatrixSetId" in kwargs.keys():
+        TileMatrixSetId = kwargs["TileMatrixSetId"]
+    if "expression" in kwargs.keys():
+        params["expression"] = kwargs["expression"]
+    if "tile_format" in kwargs.keys():
+        params["tile_format"] = kwargs["tile_format"]
+    if "tile_scale" in kwargs.keys():
+        params["tile_scale"] = kwargs["tile_scale"]
+    if "minzoom" in kwargs.keys():
+        params["minzoom"] = kwargs["minzoom"]
+    if "maxzoom" in kwargs.keys():
+        params["maxzoom"] = kwargs["maxzoom"]
+
+    allowed_bands = get_STAC_bands(url, titiler_endpoint)
+
+    if bands is None:
+        bands = [allowed_bands[0]]
+    elif len(bands) <= 3 and all(x in allowed_bands for x in bands):
+        pass
+    else:
+        raise Exception('You can only select 3 bands from the following: {}'.format(
+            ', '.join(allowed_bands))) 
+        
+    assets = ','.join(bands)
+    params["assets"] = assets
+
+    r = requests.get(
+        f"{titiler_endpoint}/stac/{TileMatrixSetId}/tilejson.json",
+        params = params
+    ).json()
+
+    return r["tiles"][0]
+
+
+def get_STAC_bounds(url, titiler_endpoint="https://api.cogeo.xyz/"):
+    """Get the bounding box of a single SpatialTemporal Asset Catalog (STAC) item.
+
+    Args:
+        url (str): HTTP URL to a STAC item, e.g., https://canada-spot-ortho.s3.amazonaws.com/canada_spot_orthoimages/canada_spot5_orthoimages/S5_2007/S5_11055_6057_20070622/S5_11055_6057_20070622.json
+        titiler_endpoint (str, optional): Titiler endpoint. Defaults to "https://api.cogeo.xyz/".
+
+    Returns:
+        list: A list of values representing [left, bottom, right, top]
+    """    
+    import requests
+
+    r = requests.get(
+        f"{titiler_endpoint}/stac/bounds",
+        params = {
+            "url": url
+        }
+    ).json()   
+
+    bounds = r["bounds"] 
+    return bounds
+
+
+def get_STAC_center(url, titiler_endpoint="https://api.cogeo.xyz/"):
+    """Get the centroid of a single SpatialTemporal Asset Catalog (STAC) item.
+
+    Args:
+        url (str): HTTP URL to a STAC item, e.g., https://canada-spot-ortho.s3.amazonaws.com/canada_spot_orthoimages/canada_spot5_orthoimages/S5_2007/S5_11055_6057_20070622/S5_11055_6057_20070622.json
+        titiler_endpoint (str, optional): Titiler endpoint. Defaults to "https://api.cogeo.xyz/".
+
+    Returns:
+        tuple: A tuple representing (longitude, latitude)
+    """    
+    bounds = get_STAC_bounds(url, titiler_endpoint)
+    center=((bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2)  # (lat, lon)
+    return center
+
+
+def get_STAC_bands(url, titiler_endpoint="https://api.cogeo.xyz/"):
+    """Get band names of a single SpatialTemporal Asset Catalog (STAC) item.
+
+    Args:
+        url (str): HTTP URL to a STAC item, e.g., https://canada-spot-ortho.s3.amazonaws.com/canada_spot_orthoimages/canada_spot5_orthoimages/S5_2007/S5_11055_6057_20070622/S5_11055_6057_20070622.json
+        titiler_endpoint (str, optional): Titiler endpoint. Defaults to "https://api.cogeo.xyz/".
+
+    Returns:
+        list: A list of band names
+    """    
+    import requests
+
+    r = requests.get(
+        f"{titiler_endpoint}/stac/info",
+        params = {
+            "url": url,
+        }
+    ).json()
+
+    return r
+
+
+def explode(coords):
+    """Explode a GeoJSON geometry's coordinates object and yield
+    coordinate tuples. As long as the input is conforming, the type of
+    the geometry doesn't matter.  From Fiona 1.4.8"""
+    for e in coords:
+        if isinstance(e, (float, int)):
+            yield coords
+            break
+        else:
+            for f in explode(e):
+                yield f
+
+
+def get_bounds(geometry, north_up=True, transform=None):
+    """Bounding box of a GeoJSON geometry, GeometryCollection, or FeatureCollection.
+    left, bottom, right, top
+    *not* xmin, ymin, xmax, ymax
+    If not north_up, y will be switched to guarantee the above.
+    Source code adapted from https://github.com/mapbox/rasterio/blob/master/rasterio/features.py#L361
+    """
+
+    if 'bbox' in geometry:
+        return tuple(geometry['bbox'])
+
+    geometry = geometry.get('geometry') or geometry  
+
+    # geometry must be a geometry, GeometryCollection, or FeatureCollection
+    if not ('coordinates' in geometry or 'geometries' in geometry or 'features' in geometry):
+        raise ValueError(
+            "geometry must be a GeoJSON-like geometry, GeometryCollection, "
+            "or FeatureCollection"
+        )
+
+    if 'features' in geometry:
+        # Input is a FeatureCollection
+        xmins = []
+        ymins = []
+        xmaxs = []
+        ymaxs = []
+        for feature in geometry['features']:
+            xmin, ymin, xmax, ymax = get_bounds(feature['geometry'])
+            xmins.append(xmin)
+            ymins.append(ymin)
+            xmaxs.append(xmax)
+            ymaxs.append(ymax)
+        if north_up:
+            return min(xmins), min(ymins), max(xmaxs), max(ymaxs)
+        else:
+            return min(xmins), max(ymaxs), max(xmaxs), min(ymins)
+
+    elif 'geometries' in geometry:
+        # Input is a geometry collection
+        xmins = []
+        ymins = []
+        xmaxs = []
+        ymaxs = []
+        for geometry in geometry['geometries']:
+            xmin, ymin, xmax, ymax = get_bounds(geometry)
+            xmins.append(xmin)
+            ymins.append(ymin)
+            xmaxs.append(xmax)
+            ymaxs.append(ymax)
+        if north_up:
+            return min(xmins), min(ymins), max(xmaxs), max(ymaxs)
+        else:
+            return min(xmins), max(ymaxs), max(xmaxs), min(ymins)
+
+    elif 'coordinates' in geometry:
+        # Input is a singular geometry object
+        if transform is not None:
+            xyz = list(explode(geometry['coordinates']))
+            xyz_px = [transform * point for point in xyz]
+            xyz = tuple(zip(*xyz_px))
+            return min(xyz[0]), max(xyz[1]), max(xyz[0]), min(xyz[1])
+        else:
+            xyz = tuple(zip(*list(explode(geometry['coordinates']))))
+            if north_up:
+                return min(xyz[0]), min(xyz[1]), max(xyz[0]), max(xyz[1])
+            else:
+                return min(xyz[0]), max(xyz[1]), max(xyz[0]), min(xyz[1])
+
+    # all valid inputs returned above, so whatever falls through is an error
+    raise ValueError(
+            "geometry must be a GeoJSON-like geometry, GeometryCollection, "
+            "or FeatureCollection"
+        )
+
+        
 def image_props(img, date_format='YYYY-MM-dd'):
     """Gets image properties.
 
diff --git a/eefolium/eefolium.py b/eefolium/eefolium.py
@@ -457,6 +457,57 @@ def add_tile_layer(
         except:
             print("Failed to add the specified TileLayer.")
 
+    def add_COG_layer(
+        self,
+        url,
+        name='Untitled',
+        attribution=".",
+        opacity=1.0,
+        shown=True,
+        titiler_endpoint = "https://api.cogeo.xyz/",
+        **kwargs
+    ):
+        """Adds a COG TileLayer to the map.
+
+        Args:
+            url (str): The URL of the COG tile layer. 
+            name (str, optional): The layer name to use for the layer. Defaults to 'Untitled'.
+            attribution (str, optional): The attribution to use. Defaults to '.'.
+            opacity (float, optional): The opacity of the layer. Defaults to 1.
+            shown (bool, optional): A flag indicating whether the layer should be on by default. Defaults to True.
+            titiler_endpoint (str, optional): Titiler endpoint. Defaults to "https://api.cogeo.xyz/".
+        """
+        tile_url = get_COG_tile(url, titiler_endpoint, **kwargs)
+        center= get_COG_center(url, titiler_endpoint)  # (lon, lat)
+        self.add_tile_layer(tiles=tile_url, name=name, attribution=attribution, opacity=opacity, shown=shown)
+        self.set_center(lon=center[0], lat=center[1], zoom=10)
+
+    def add_STAC_layer(
+        self,
+        url,
+        bands=None,
+        name='Untitled',
+        attribution=".",
+        opacity=1.0,
+        shown=True,
+        titiler_endpoint = "https://api.cogeo.xyz/",
+        **kwargs
+    ):
+        """Adds a STAC TileLayer to the map.
+
+        Args:
+            url (str): The URL of the COG tile layer. 
+            name (str, optional): The layer name to use for the layer. Defaults to 'Untitled'.
+            attribution (str, optional): The attribution to use. Defaults to '.'.
+            opacity (float, optional): The opacity of the layer. Defaults to 1.
+            shown (bool, optional): A flag indicating whether the layer should be on by default. Defaults to True.
+            titiler_endpoint (str, optional): Titiler endpoint. Defaults to "https://api.cogeo.xyz/".
+        """
+        tile_url = get_STAC_tile(url, bands, titiler_endpoint, **kwargs)
+        center= get_STAC_center(url, titiler_endpoint)
+        self.add_tile_layer(tiles=tile_url, name=name, attribution=attribution, opacity=opacity, shown=shown)
+        self.set_center(lon=center[0], lat=center[1], zoom=10)
+        
     def publish(
         self,
         name=None,
