<jupyter_start><jupyter_text>Natural colour composites using Lansdat 7 imagery Annual basemap layers for the years 2013-2016 were produced with Landsat 8 imagery and are already available in the Global Forest Watch platform. In this notebook we'll explore different approaches to generate annual basemaps for the years 2000-2014 using Landsat 7. We will aim to replicate the natural-color tones of the 2013-2016 composites, and also explore ways to improve the rendering speed for high zoom level tiles.The main purpose of this notebook is just to give a preview of the look and feel of the resulting composites. A separate document will be prepared to describe the workflow for exporting tiles to GCS and into the microservice.If you're viewing this notebook form GH you'll need to wait a few minutes for all the images to show up.You can skip this introduction and head over sections 3, 4 and 5, which show the most promising results for rendering tiles at different z-resolutions<jupyter_code>import os
import ee
import json
import ipywidgets
from IPython.display import Image
import folium
import pprint
import datetime
from dateutil.relativedelta import relativedelta

ee.Initialize()<jupyter_output><empty_output><jupyter_text>Setup visualisation toolsWe will first define a few functions to help us visualise the results, these include:-  `tile_url`: generates tile urls from images computed in Earth Engine-  `dateRange`: returns a custom date range calculated around a target year -  `toImage`: retrieves images from composites computed in Earth Engine and plots them in 3 different zoom levels.<jupyter_code># Function to plot tiles (Images)
def tile_url(image, viz_params=None):
    """Create a target url for tiles for an image.
    e.g.
    im = ee.Image("LE7_TOA_1YEAR/" + year).select("B3","B2","B1")
    viz = {'opacity': 1, 'gain':3.5, 'bias':4, 'gamma':1.5}
    url = tile_url(image=im,viz_params=viz)
    """
    if viz_params:
        d = image.getMapId(viz_params)
    else:
        d = image.getMapId()
    base_url = 'https://earthengine.googleapis.com'
    url = (base_url + '/map/' + d['mapid'] + '/{z}/{x}/{y}?token=' + d['token'])
    return url
# Function to generate a date range
def dateRange(mid_year,period):
    '''
    Params:
        mid_year (assumes yyyy): Mid point of the time range 
        period (int): number of years
    Returns:
        dictionary with a 'start' and an 'end' date in a `yyyy-mm-dd` format.
        The reference year is also provided in the 'central' key
    '''
    # define central year, commencing in JAN-01
    mid = datetime.datetime.strptime(str(mid_year) + '-01-01', '%Y-%m-%d')
    
    # get half period
    half_period = int(period/2)
    
    # compute date range
    if period%2 ==1:
        end = mid + relativedelta(years= half_period,months=6)
        start = mid - relativedelta(years= half_period,months=6)
    else:
        end = mid + relativedelta(years= half_period)
        start = mid - relativedelta(years= half_period)
    
    dates = {'start':start.strftime('%Y-%m-%d'),'central':str(mid_year) ,'end':end.strftime('%Y-%m-%d')}
    
    print('Returning '+ str(period) + ' years range: from '+ dates['start'] + ' to ' + dates['end'])
    
    return(dates)
# Function to produce image
def toImage(composite,visParam,title):
    '''
    Retrieves Google Earth image's thumbUrls and plot the results at 3 different zoom levels. Note that the `region_z` objects are pulled from the global environment.
    composite: ee.Image (assumes 3 bands)
    visParam: image.vizualisation params e.g. {'opacity': 1, 'gain':3.5, 'bias':4, 'gamma':1.5}
    title: plot
    '''
    visParam['region'] = region_z10 
    print('Z10: '+ title + ': ' + dates['start'] + ' - ' + dates['end'])
    image1 = Image(url=composite.getThumbUrl(visParam))
    display(image1)
    
    visParam['region'] = region_z6 
    print('Z6: '+ title + ': ' + dates['start'] + ' - ' + dates['end'])
    image2= Image(url=composite.getThumbUrl(visParam))
    display(image2)
    
    visParam['region'] = region_z4 
    print('Z4: '+ title + ': ' + dates['start'] + ' - ' + dates['end'])
    image3= Image(url=composite.getThumbUrl(visParam))
    display(image3)<jupyter_output><empty_output><jupyter_text>Define ROIsBelow we will calculate buffer zones from a given lon-lat coordinate and generate 3 ROI objects. These will then be used in the `toImage` function to print Lansat images in the notebook. These ROIs are roughly equivalent to zoom levels 4, 6 and 10.Some sites of interest:|site|ISO|lon|lat||-----------|-------|---|---||El Triunfo|MEX|-92.7548|15.6151||Chimalapas|MEX|-94.3835|16.7604||Leticia|COL|-70.2727|-3.8094||Oslo|NOR|10.679|59.955||Okavango Delta|BWA|22.346|-18.792||Gran Canaria|ESP|-15.5786|27.9507|<jupyter_code># --- Coordinates
lon = -15.5786
lat = 27.9507
zstart = 10
# -- Set ROI for IPython.display
point = ee.Geometry.Point(lon,lat)

# set bufffers:units is meters 
# resulting zoom levels are an approximation
buffer_z10 = point.buffer(40000)
buffer_z6 = point.buffer(650000)
buffer_z4 = point.buffer(2500000)

# convert to FeatureCollection
bufferFeat_z10 = ee.FeatureCollection(buffer_z10).set('area',1)
bufferFeat_z6 = ee.FeatureCollection(buffer_z6).set('area',1)
bufferFeat_z4 = ee.FeatureCollection(buffer_z4).set('area',1)

# extract bounds
region_z10 = buffer_z10.bounds().getInfo()['coordinates']
region_z6 = buffer_z6.bounds().getInfo()['coordinates']
region_z4 = buffer_z4.bounds().getInfo()['coordinates']

#  -- preview ROI
# roi = tile_url(bufferFeat_z10,{'color': '#FF5733'})
# map = folium.Map(location=[lat,lon],zoom_start=zstart)
# map.add_tile_layer(tiles=roi, attr="na")
# map<jupyter_output><empty_output><jupyter_text>Landsat productsEarth Engine provides a variety of [Landsat products](https://developers.google.com/earth-engine/landsat) which are cataloged according to their data quality and processing methods. Some of these proved to be very useful at certain zoom levels when producing the basemaps for the 2013-2016 period. In the sections below we'll apply a similar compositing methodology with the equivalent Landsat 7 products. These include the following collections: Orthorectified**USGS Landsat 7 Raw Scenes (Orthorectified)**- id: LANDSAT/LE7_L1T- [description](https://code.earthengine.google.com/dataset/LANDSAT/LE7_L1T)- dates: Jan 1, 1999 - Apr 30, 2017 Collection 1 Tier 1Landsat scenes with the highest available data quality are placed into Tier 1 and are considered suitable for time-series processing analysis. The T1_RT collection contains both Tier 1 and Real-Time (RT) assets.**USGS Landsat 7 Collection 1 Tier 1 Raw Scenes**- id: LANDSAT/LE07/C01/T1- [description](https://code.earthengine.google.com/dataset/LANDSAT/LE07/C01/T1)- dates: Jan 1, 1999 - Present**USGS Landsat 7 Collection 1 Tier 1 and Real-Time data aw Scenes]**- id: LANDSAT/LE07/C01/T1_RT- [description](https://code.earthengine.google.com/dataset/LANDSAT/LE07/C01/T1_RT)- dates: Jan 1, 1999 - Present Calibrated TOAThese Landsat 7 composites are made from Level L1T orthorectified scenes, using the computed top-of-atmosphere (TOA) reflectance. As of May 1, 2017, the USGS is no longer producing Pre-Collection Landsat, and therefore these collections are on longer being updated.**USGS Landsat 7 Collection 1 Tier 1 and Real-Time data OA Reflectance]**- id: LANDSAT/LE07/C01/T1_RT_TOA- [description](https://code.earthengine.google.com/dataset/LANDSAT/LE07/C01/T1_RT_TOA)- dates: Jan 1, 1999 - Present**USGS Landsat 7 TOA Reflectance (Orthorectified)**- id: LANDSAT/LE7_L1T_TOA- [description](https://code.earthengine.google.com/dataset/LANDSAT/LE7_L1T_TOA)- dates: Jan 1, 1999 - Apr 30, 2017**Landsat 7 Annual Greenest-Pixel TOA Reflectance Composite**- id: LE7_L1T_ANNUAL_GREENEST_TOA- [description](https://code.earthengine.google.com/dataset/LANDSAT/LE7_L1T_ANNUAL_GREENEST_TOA)- dates: Jan 1, 1999 - Jan 1, 2018 ---------  1. USGS Landsat 7 Raw Scenes (Orthorectified) ` LANDSAT/LE7_L1T`<jupyter_code>dates = dateRange(2010,3)

# load collection and filter by dates
l7_l1t = ee.ImageCollection('LANDSAT/LE7_L1T').filterDate(dates['start'], dates['end'])<jupyter_output>Returning 3 years range: from 2008-07-01 to 2011-07-01<jupyter_text>Cloud free compositeReturns a cloud-free composite computed with the `ee.Algorithms.Landsat.simpleComposite()`. This function converts raw scenes to TOA reflectance, applies the simple cloud score and takes a given percentile of the least cloudy pixels.[Sample code](https://code.earthengine.google.com/e8832e9bf4d073cda0dcdb6638a065ae):```var collection = ee.ImageCollection('LANDSAT/LC8_L1T').filterDate('2015-01-01T00:00','2015-12-31T00:00');var composite = ee.Algorithms.Landsat.simpleComposite({collection: collection, percentile: 50, maxDepth:80, cloudScoreRange: 1, asFloat:true});var hsv2 = composite.select(['B4', 'B3', 'B2']).rgbToHsv();var median = hsv2.hsvToRgb().visualize({gain:1000, gamma: [1.15, 1.4, 1.15]});   ```<jupyter_code># cutsom composite
## - Playing with the percentile values results in different brightness intensities.
l7_l1t_composite = ee.Algorithms.Landsat.simpleComposite(
    collection= l7_l1t, 
    percentile= 30, 
    maxDepth=80, 
    cloudScoreRange= 1, asFloat=True)
    
# to HSV
hsv_l1t = l7_l1t_composite.select(['B3', 'B2', 'B1']).rgbToHsv()

# make image
rgb_l1t= hsv_l1t.hsvToRgb().visualize()
visParam = {
    'gain':3,
    'gamma':'1.15,1.4,1.15',
}

# plot
toImage(rgb_l1t,visParam,'simpleComposite (custom parameters)')
# Map tiles
# ee_tiles = tile_url(rgb_l1t,visParam)
# map = folium.Map(location=[lat,lon], zoom_start=5)
# map.add_tile_layer(tiles=ee_tiles, attr="NASA")
# map<jupyter_output><empty_output><jupyter_text>Cloud free composite + pan-sharpenNow we'll pansharpen the above composite.[Sample code](https://code.earthengine.google.com/cbe5822c7166f99fc414d59d8e70c7fc)<jupyter_code># Sharpen above compostie
l7_l1t_sharpened = ee.Image.cat([hsv_l1t.select('hue'), hsv_l1t.select('saturation'), 
                          l7_l1t_composite.select('B8')]).hsvToRgb() # .visualize(gain=600,gamma= [1.15, 1.4, 1.15])
visParam = {
    'gain':1000,
    'gamma':'1.15,1.4,1.15',
}

# plot
toImage(l7_l1t_sharpened,visParam,'simpleComposite (custom parameters) + pansharpen')
# Map tiles
# ee_tiles = tile_url(l7_l1t_sharpened,visParam)
# map = folium.Map(location=[lat,lon], zoom_start=11)
# map.add_tile_layer(tiles=ee_tiles, attr="NASA")
# map<jupyter_output><empty_output><jupyter_text>2. USGS Landsat 7 TOA Reflectance (Orthorectified)`LANDSAT/LE7_L1T_TOA`<jupyter_code>dates = dateRange(2010,4)

# load collection and filter by dates
l7_l1t_toa = ee.ImageCollection('LANDSAT/LE7_L1T_TOA').filterDate(dates['start'], dates['end'])<jupyter_output>Returning 4 years range: from 2008-01-01 to 2012-01-01<jupyter_text>Filter by cloud cover and median over time + pan-sharpenThis method selects a subset of scenes at a given location, applies a cloud cover filter and takes the median value of the resulting pixels. Since we're working with a TOA reflectance data we don't need to calculate this parameter. We'll also pan-sharpen these images so results can be used for high resolution tiles.The original workflow for processing Landsat 8 data included a `CLOUD_COVER_LAND` filter, which is used to calculate could scores of the [land pixels in a scene](https://landsat.usgs.gov/what-difference-between-scene-cloud-cover-and-land-cloud-cover-scores). This parameter is not available in the Landsat 7, so instead I'm applying a simple `CLOUD_COVER` score.The results are not of much use though...[Source](https://code.earthengine.google.com/2f72c2a15bd46646329e28c13ba2a558)```// Get the median over time, in each band, in each pixel.var median = ee.ImageCollection('LANDSAT/LC8_L1T_TOA')                  .filterDate('2013-04-01', '2016-12-31').filter(ee.Filter.lt('CLOUD_COVER_LAND', 35)).median();var hsv = median.select(['B4', 'B3', 'B2']).rgbToHsv();    var sharpened = ee.Image.cat([hsv.select('hue'), hsv.select('saturation'), median.select('B8')]).hsvToRgb().visualize({gain:1000, gamma: [1.15, 1.4, 1.15]});```<jupyter_code>l7_l1t_toa_med = l7_l1t_toa.filter(ee.Filter.lt('CLOUD_COVER', 7)).median()
hsv_l1t_toa = l7_l1t_toa_med.select(['B3', 'B2', 'B1']).rgbToHsv()
# Sharpen above compostie
l7_l1t_toa_sharpened = ee.Image.cat([hsv_l1t_toa.select('hue'), hsv_l1t_toa.select('saturation'),
                                     l7_l1t_toa_med.select('B8')]).hsvToRgb()
visParam= {
    'gain': 800, 
    'gamma':'1.15, 1.4, 1.15'
}

# plot
toImage(l7_l1t_toa_sharpened,visParam,'Filter by cloud & median values + Pansharpened')<jupyter_output>Z10: Filter by cloud & median values + Pansharpened: 2008-01-01 - 2012-01-01<jupyter_text>3. USGS Landsat 7 Collection 1 Tier 1 Raw Scenes `LANDSAT/LE07/C01/T1`<jupyter_code>dates = dateRange(2010,3)

# load collection and filter by dates
l7_t1 = ee.ImageCollection('LANDSAT/LE07/C01/T1').filterDate(dates['start'], dates['end'])<jupyter_output>Returning 3 years range: from 2008-07-01 to 2011-07-01<jupyter_text>Cloud free composite + pan-sharpenIn this example we will generate a composite with the `ee.Algorithms.Landsat.simpleComposite()` module. We'll provide with a custom cloud cover threshold and percentile value to filter the resulting pixels. This particular combination looks very useful for lower resolution tiles (Z6), and the color palette is similar to the Landsat basemaps in GFW for recent years[Source](https://developers.google.com/earth-engine/landsatsimple-composite)<jupyter_code># Create a cloud-free composite with custom parameters for
customComposite = ee.Algorithms.Landsat.simpleComposite(
  collection= l7_t1,
    percentile= 30,
    cloudScoreRange= 8,
    asFloat= True)
# Convert the RGB bands to the HSV color space.
hsv_custom = customComposite.select(['B3', 'B2', 'B1']).rgbToHsv()

# now from HSV to RGB
rgb_custom = hsv_custom.hsvToRgb()
# Define visualisation parameters
visParam= {
    'gain': 700, 
    'gamma':'1.15, 1.4, 1.15'
}

# plot
toImage(rgb_custom,visParam,'Simple composite (custom parameters)')
# Map tiles
# ee_tiles = tile_url(rgb_custom,visParam)  # hsv_custom looks cool
# map = folium.Map(location=[lat,lon], zoom_start=7)
# map.add_tile_layer(tiles=ee_tiles, attr="NASA")
# map<jupyter_output><empty_output><jupyter_text>Increasing the date query to 5 years and changing the visualisation parameters (see below) returns composites that could also be used for higher zoom levels.<jupyter_code>dates = dateRange(2010,5)  # <- increase date range

# load collection and filter by dates
l7_t1 = ee.ImageCollection('LANDSAT/LE07/C01/T1').filterDate(dates['start'], dates['end'])


# Create a cloud-free composite with custom parameters for
customComposite = ee.Algorithms.Landsat.simpleComposite(
  collection= l7_t1,
    percentile= 60,
    cloudScoreRange= 1,
    asFloat= True)

# Define visualisation parameters
visParam= {
    'gain': 800, 
    'gamma':'1.15, 1.40, 1.15'
}

# plot results
toImage(rgb_custom,visParam,'Simple composite (custom parameters)')<jupyter_output>Returning 5 years range: from 2007-07-01 to 2012-07-01
Z10: Simple composite (custom parameters): 2007-07-01 - 2012-07-01<jupyter_text>4. USGS Landsat 7 Collection 1 Tier 1 and Real-Time data Raw Scenes`LE07/C01/T1_RT/`Tier 1 collections contain images with the highest quality available. In this section we'll compute 2 composites using the  `simpleComposite` algorithm, one with default parameters and the other with a custom cloud score and percentile value[Source](https://code.earthengine.google.com/9c32b62c959c1f4121cb3218380afb7d)<jupyter_code>dates = dateRange(2010,4)

# load collection and filter by dates
l7_t1_rt = ee.ImageCollection('LANDSAT/LE07/C01/T1_RT').filterDate(dates['start'], dates['end'])
l7_t1_rt_simpleComp = ee.Algorithms.Landsat.simpleComposite(l7_t1_rt)
# Define visParams
visParam = {
    'bands':'B3,B2,B1',
    'gain':2.5,
    'gamma':'1.3, 1.3, 1.3',
#     'dimensions': 800
}

# plot
toImage(l7_t1_rt_simpleComp,visParam,'Simple composite (default params)')
# Create a cloud-free composite with custom parameters for
customComposite = ee.Algorithms.Landsat.simpleComposite(
  collection= l7_t1_rt,
    percentile= 50,
    cloudScoreRange= 1,
    asFloat= True)
# Convert the RGB bands to the HSV color space.
hsv_custom = customComposite.select(['B3', 'B2', 'B1']).rgbToHsv()

# now from HSV to RGB
rgb_custom = hsv_custom.hsvToRgb()
# Define visParams
visParam = {
    'gain':700, 
    'gamma':'1.15, 1.4, 1.15'
}

# plot
toImage(rgb_custom,visParam,'C1 Tier1 Pansharepened')
# Map tiles
# ee_tiles = tile_url(rgb_custom,visParam)
# map = folium.Map(location=[lat,lon], zoom_start=9)
# map.add_tile_layer(tiles=ee_tiles, attr="NASA")
# map.add_tile_layer(tiles=feat_tiles, attr="na")
# map<jupyter_output><empty_output><jupyter_text>5. USGS Landsat 7 Collection 1 Tier 1 and Real-Time data TOA Reflectance` LANDSAT/LE07/C01/T1_RT_TOA`<jupyter_code>dates = dateRange(2010,4)
# load collection and filter by dates
l7_t1_rt_toa = ee.ImageCollection('LANDSAT/LE07/C01/T1_RT_TOA').filterDate(dates['start'], dates['end'])<jupyter_output>Returning 4 years range: from 2008-01-01 to 2012-01-01<jupyter_text>Filter by cloud cover and median over time + pan-sharpenThis method selects a subset of scenes at a given location, applies a cloud cover filter and takes the median value of the resulting pixels. Since we're working with a TOA reflectance data we don't need to calculate this parameter. Finally we'll also pan-sharpen these imagesThe results could be of potential use for high resolution tiles<jupyter_code>l7_t1_rt_toa_med = l7_t1_rt_toa.filter(ee.Filter.lt('CLOUD_COVER', 10)).median()
# Define visParams
visParam = {
    'bands':'B3,B2,B1',
    'min':0,
    'max':0.25,
#     'dimensions': 800
}


# plot
toImage(l7_t1_rt_toa_med,visParam,'Filter by cloud cover | median')<jupyter_output>Z10: Filter by cloud cover | median: 2008-01-01 - 2012-01-01<jupyter_text>6. Landsat 7 Annual Greenest-Pixel TOA Reflectance Composite`LE7_L1T_ANNUAL_GREENEST_TOA`<jupyter_code>dates = dateRange(2010,5)

# load collection and filter by dates
l7_toa_greenest_annual = ee.ImageCollection('LE7_L1T_ANNUAL_GREENEST_TOA').filterDate(dates['start'], dates['end'])<jupyter_output>Returning 5 years range: from 2007-07-01 to 2012-07-01<jupyter_text>Greenest pixel, and Greenest + PansharpLastly we'll experiment with the greenest-pixel collection. We'll compute a simple composite with the median values of a given period, and will also generate a pan-sharpened version of this composite[Source](https://code.earthengine.google.com/3b230f4c217c589ab267e438798b5e28)    ```//--- Greenest pixel// NDVI Functionvar addNDVI = function(image) {var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');  return image.addBands(ndvi);};var withNDVI = l8.map(addNDVI);// Make a "greenest" pixel composite.var greenest = withNDVI.qualityMosaic('NDVI');// Display the result.var visParams = {bands: ['B3', 'B2', 'B1'], max: 0.3};Map.addLayer(greenest, visParams, 'Greenest pixel composite');// ---- Greenest + Pansharpvar hsv = greenest.select(['B3', 'B2', 'B1']).rgbToHsv();    var sharpened = ee.Image.cat([hsv.select('hue'), hsv.select('saturation'), greenest.select('B8')]).hsvToRgb().visualize({gain:1000, gamma: [1.15, 1.4, 1.15]});Map.addLayer(sharpened,null,'pan+green');```<jupyter_code># get media (filtering by CLOUD_COVER is not available in this collection)
l7_toa_greenest_annual = l7_toa_greenest_annual.median()

# Plot
visParam = {
    'bands':'B3,B2,B1',
    'max':0.3
}

# plot
toImage(l7_toa_greenest_annual,visParam,'Greenest-pixel')
# calculate HSV
hsv_greenest = l7_toa_greenest_annual.select(['B3', 'B2', 'B1']).rgbToHsv()

# Sharpen above compostie
l7_toa_greenest_sharpened = ee.Image.cat([hsv_greenest.select('hue'), hsv_greenest.select('saturation'),
                                     l7_toa_greenest_annual.select('B8')]).hsvToRgb()
# Plot
visParam = {
    'gain':1000, 
    'gamma': '1.15, 1.4, 1.15'
}

# plot
toImage(l7_toa_greenest_sharpened,visParam,'Greenest-pixel + pansharpen')<jupyter_output>Z10: Greenest-pixel + pansharpen: 2007-07-01 - 2012-07-01
{'fraction_non_alphanumeric': 0.1251485736925515, 'fraction_numerical': 0.05219889064976228, 'mean_word_length': 5.063963963963964, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 56, 'https://': 15, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '465137', 'license': 'MIT', 'n_tokens_mistral': 7730, 'n_tokens_neox': 6946, 'n_words': 1790, 'path': 'work/Landsat 7 annual basemaps - Pt 1.ipynb', 'repo_name': 'Vizzuality/data_sci_tutorials', 'size': '20192'}