better-satcen-00004 data pipeline results (Sentinel-2 DVI, composites and BOA reflectances): loading from pickle file and exploiting

This Notebook shows how to: * reload better-satcen-00004 data pipeline results previously stored in a pickle file * Create a new dataframe with data values cropped over the area of interest * Plot data * Create a geotiff * Download data

Import needed modules

In [1]:

import pandas as pd
from geopandas import GeoDataFrame
import gdal
import numpy as np
from shapely.geometry import Point
from shapely.geometry import Polygon
import matplotlib
import matplotlib.pyplot as plt
from PIL import Image
%matplotlib inline
from shapely.wkt import loads
from shapely.geometry import box
from urlparse import urlparse
import requests

Read picke data file

In [2]:

pickle_filename = 'better-satcen-00004_Jun-Jul_2018_Afghanistan.pkl'

In [3]:

results = GeoDataFrame(pd.read_pickle(pickle_filename))

In [4]:

results.head()


Out[4]:
enclosure enddate identifier self startdate title wkt
0 https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 01be9eec66cb994310e469e012a3a2bc670f92e8 https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_2... POLYGON((64.0648798642054 31.5308017399005,65....
1 https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 037133bb420b2f5407a223e2f2709e7c6eefbdf0 https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41RNP_2... POLYGON((63.9421703624035 30.4373158799178,64....
2 https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 03aa38a9b430af9f2b2b2cb47cae31788c0fd59c https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_2... POLYGON((63.9421873560077 30.6397709351317,64....
3 https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 0949c278e99a1b2244ddd0edb3be4e38c7b10e81 https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41SNR_2... POLYGON((63.9422209286482 31.5415371135011,64....
4 https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 0cccedf5d105c98b1bce071eebd3639084ebc43d https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_2... POLYGON((63.9421873560077 30.6397709351317,64....

Credentials for ellip platform access

  • Provide here your ellip username and api key for platform access
In [5]:

import getpass

user = getpass.getpass("Enter your username:")
api_key = getpass.getpass("Enter the API key:")

Define filter parameters

Time of interest

In [6]:

mydates = ['2018-07-29 06:16:31.024', '2018-07-27 06:26:29.024','2018-07-24 06:16:29.024', '2018-07-22 06:26:31.024','2018-07-19 06:16:31.024']

Area of interest

  • The user can choose to define an AOI selecting a Point and a buffer size used to build a squared polygon around that point
In [7]:

point_of_interest = Point(64.108, 31.586)

In [8]:

buffer_size = 0.07

In [9]:

aoi_wkt1 = box(*point_of_interest.buffer(buffer_size).bounds)

In [10]:

aoi_wkt1.wkt

Out[10]:

'POLYGON ((64.178 31.516, 64.178 31.656, 64.03800000000001 31.656, 64.03800000000001 31.516, 64.178 31.516))'
  • Or define a Polygon (in this case this is exactly the Afghanistan reference AOI) from a WKT geometry
In [11]:

geom = 'MULTIPOLYGON (((-77.24516666666666 0.9424444444444445, -77.24516666666666 2.753833333333333, -79.0253888888889 2.753833333333333, -79.0253888888889 0.9424444444444445, -77.24516666666666 0.9424444444444445)), ((21.29611111111111 39.58638888888889, 21.29611111111111 41.032, 19.89788888888889 41.032, 19.89788888888889 39.58638888888889, 21.29611111111111 39.58638888888889)), ((65.02055555555556 30.43894444444445, 65.02055555555556 33.39566666666666, 63.94222222222222 33.39566666666666, 63.94222222222222 30.43894444444445, 65.02055555555556 30.43894444444445)))'

# Afghanistan area
aoi_wkt2 = loads(geom)[2]


In [12]:

aoi_wkt2.wkt

Out[12]:

'POLYGON ((65.02055555555556 30.43894444444445, 65.02055555555556 33.39566666666666, 63.94222222222222 33.39566666666666, 63.94222222222222 30.43894444444445, 65.02055555555556 30.43894444444445))'

Create a new dataframe with data values

  • Get the subframe data values related to selected tiles and bands, and create a new GeoDataFrame having the original metadata stack with the specific bands data extended info (data values, size, and geo projection and transformation).

Define auxiliary methods to create new dataframe

  • The first method gets data url reference, bands list and names, cropping flag, AOI cropping area as inputs and returns the related n bands data array with related geospatial infos
  • If crop = False the original bands extent is returned
In [13]:

def get_bands_as_array(url,band_indices,bands_name,crop,bbox):

    output = '/vsimem/clip.tif'
    bands=[]

    ds = gdal.Open('/vsicurl/%s' % url)

    bands_count = ds.RasterCount

    if bands_count < len(band_indices):
        raise('ERROR: Selected bands are more than the data ones')

    if crop == True:
        ulx, uly, lrx, lry = bbox[0], bbox[3], bbox[2], bbox[1]
        ds = gdal.Translate(output, ds, bandList=band_indices, projWin = [ulx, uly, lrx, lry], projWinSRS = 'EPSG:4326')

    else:
        ds = gdal.Translate(output, ds, bandList=band_indices)

    ds = None
    ds = gdal.Open(output)
    w = ds.GetRasterBand(1).XSize
    h = ds.GetRasterBand(1).YSize
    geo_transform = ds.GetGeoTransform()
    projection = ds.GetProjection()

    for i in range(len(band_indices)):

        data = ds.GetRasterBand(i+1).ReadAsArray(0, 0, w, h).astype(np.float)

        data[data==-9999] = np.nan

        bands.append(data)

        band = None
        data = None


    ds = None

    return bands,geo_transform,projection,w,h
  • The second method gets reference url (for composites data), cropping flag and AOI cropping area as inputs and returns the related n bands data array with geospatial infos
  • If crop = False the original bands extent is returned
In [14]:

def get_composites_bands_as_array(url,crop,bbox):

    output = '/vsimem/clip.tif'

    ds = gdal.Open('/vsicurl/%s' % url)


    if crop == True:
        ulx, uly, lrx, lry = bbox[0], bbox[3], bbox[2], bbox[1]
        ds = gdal.Translate(output, ds, projWin = [ulx, uly, lrx, lry], projWinSRS = 'EPSG:4326')
        ds = None
        ds = gdal.Open(output)

    geo_transform = ds.GetGeoTransform()
    projection = ds.GetProjection()

    data = ds.ReadAsArray().astype(np.float)

    h = data.shape[1]
    w = data.shape[2]

    data[data==-9999] = np.nan

    ds = None
    band = None

    return data,geo_transform,projection,w,h
  • The third one selects the data to be managed (only the actual result related to the selected tile) and returns a new GeoDataFrame containing the requested bands extended info with the original metadata set
In [15]:

def get_GDF_with_datavalues(row, user, api_key, bands_list='all', crop=False, bbox=None):


    bands_name = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B9', 'B11', 'B12','dvi', 'quality_cloud_confidence', 'quality_snow_confidence', 'quality_scene_classification']
    band_indices = range(1,len(bands_name)+1)

    if bands_list is not 'all':
        band_indices = [bands_name.index(x)+1 for x in bands_list ]

    # get data values: everything but the mask tif result

    parsed_url = urlparse(row['enclosure'])

    url = '%s://%s:%s@%s/api%s' % (list(parsed_url)[0], user, api_key, list(parsed_url)[1], list(parsed_url)[2])

    extended_info = list()

    if 'DVI_BOA.tif' in row['title']:
        print 'Getting bands %s for %s' %(band_indices, row['title'])
        bands,geo_transform,projection,w,h = get_bands_as_array(url,band_indices,bands_name,crop,bbox)

        for i,bi in enumerate(band_indices):
            extended_info.append(dict(row))
            extended_info[i]['band'] = bands_name[bi-1]
            extended_info[i]['data'] = bands[i]
            extended_info[i]['geo_transform'] = geo_transform
            extended_info[i]['projection'] = projection
            extended_info[i]['xsize'] = w
            extended_info[i]['ysize'] = h
    else:
        print 'Getting composite bands %s' %row['title']
        data,geo_transform,projection,w,h = get_composites_bands_as_array(url,crop,bbox)

        extended_info.append(dict(row))
        extended_info[0]['band'] = '_'.join(row['title'].split('.')[0].split('_')[7:])
        extended_info[0]['data'] = data
        extended_info[0]['geo_transform'] = geo_transform
        extended_info[0]['projection'] = projection
        extended_info[0]['xsize'] = w
        extended_info[0]['ysize'] = h

    #print 'Data values for %s : GOT' %row['title']
    return extended_info
  • Define the reference tiles
In [16]:

ref_tiles = ['T41SPR', 'T41RNQ', 'T41RNP', 'T41SNS', 'T41SPT', 'T41SPS', 'T41RPQ', 'T41SNR', 'T41RPP', 'T41SNT']
  • Define the reference bands list (we want to select)
In [17]:

b_list = ['B4','dvi']
  • Select from metadataframe a subframe of results related to more tiles on the different dates:
In [18]:

mydate_results = results[(results.apply(lambda row: row['startdate'] in pd.to_datetime(mydates[0:2]) , axis=1)) &
                         ((results.apply(lambda row: any(ref_tile in row['title'] for ref_tile in ref_tiles[0:2]), axis=1)) |
                          (results.apply(lambda row: any(b in row['title'] for b in b_list), axis=1)))]

mydate_results.head()

Out[18]:
enclosure enddate identifier self startdate title wkt
0 https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 01be9eec66cb994310e469e012a3a2bc670f92e8 https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_2... POLYGON((64.0648798642054 31.5308017399005,65....
2 https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 03aa38a9b430af9f2b2b2cb47cae31788c0fd59c https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_2... POLYGON((63.9421873560077 30.6397709351317,64....
3 https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 0949c278e99a1b2244ddd0edb3be4e38c7b10e81 https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41SNR_2... POLYGON((63.9422209286482 31.5415371135011,64....
4 https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 0cccedf5d105c98b1bce071eebd3639084ebc43d https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_2... POLYGON((63.9421873560077 30.6397709351317,64....
5 https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 13f7b9c57f80b8239602fa8643e850e2810dbe86 https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_2... POLYGON((63.9421873560077 30.6397709351317,64....
  • Create a new GeoDataFrame with the selected metadata plus the data values of the selected bands (set ‘all’ in order to select all the bands) cropped over the aoi_wkt1 area
In [19]:

myGDF = GeoDataFrame()
aoi_wkt = aoi_wkt1

for row in mydate_results.iterrows():
    myGDF = myGDF.append(get_GDF_with_datavalues(row[1], user, api_key, b_list, True, bbox=aoi_wkt.bounds), ignore_index=True)
    #myGDF = myGDF.append(get_GDF_with_datavalues(row[1], user, api_key, 'all', True, bbox=aoi_wkt.bounds), ignore_index=True)


Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_20180729T082519_DVI_BOA_B8A_B11_B2.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_20180729T082519_DVI_BOA_B8A_B11_B2.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNR_20180729T082519_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_20180729T082519_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_20180729T082519_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RPQ_20180729T082519_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNS_20180729T082519_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNT_20180729T082519_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPS_20180729T082519_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RPQ_20180729T082519_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNP_20180729T082519_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNR_20180729T082519_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_20180729T082519_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNP_20180729T082519_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_20180729T082519_DVI_BOA_B12_B4_B3.tif
Getting bands [4, 13] for S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_20180729T082519_DVI_BOA.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPT_20180729T082519_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RPQ_20180729T082519_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNP_20180729T082519_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNT_20180729T082519_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNR_20180729T082519_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNS_20180729T082519_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNS_20180729T082519_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RPP_20180729T082519_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNR_20180729T082519_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPT_20180729T082519_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RPP_20180729T082519_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPS_20180729T082519_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPT_20180729T082519_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPS_20180729T082519_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNP_20180729T082519_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_20180729T082519_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPT_20180729T082519_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNS_20180729T082519_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RPP_20180729T082519_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPT_20180729T082519_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPT_20180729T082519_DVI_BOA_B8_B4_B3.tif
Getting bands [4, 13] for S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_20180729T082519_DVI_BOA.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_20180729T082519_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_20180729T082519_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNT_20180729T082519_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNP_20180729T082519_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RPP_20180729T082519_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_20180729T082519_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RPQ_20180729T082519_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNS_20180729T082519_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNR_20180729T082519_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPS_20180729T082519_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNP_20180729T082519_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_20180729T082519_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPS_20180729T082519_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNT_20180729T082519_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RPQ_20180729T082519_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNR_20180729T082519_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_20180729T082519_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNT_20180729T082519_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNT_20180729T082519_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RPP_20180729T082519_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RPQ_20180729T082519_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SNS_20180729T082519_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_20180729T082519_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPS_20180729T082519_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_20180729T082519_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2A_MSIL2A_20180729T061631_N0206_R034_T41RPP_20180729T082519_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RPQ_20180727T091520_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNR_20180727T091520_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNR_20180727T091520_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPR_20180727T091520_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RPQ_20180727T091520_DVI_BOA_B8_B4_B3.tif
Getting bands [4, 13] for S2B_MSIL2A_20180727T062629_N0206_R077_T41SPR_20180727T091520_DVI_BOA.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNR_20180727T091520_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNR_20180727T091520_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPR_20180727T091520_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RPQ_20180727T091520_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPR_20180727T091520_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPR_20180727T091520_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNR_20180727T091520_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPR_20180727T091520_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPR_20180727T091520_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNR_20180727T091520_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPR_20180727T091520_DVI_BOA_B8A_B11_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPT_20180727T083305_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNS_20180727T083305_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNP_20180727T091520_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPT_20180727T083305_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPS_20180727T091520_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNS_20180727T091520_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNS_20180727T091520_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RPQ_20180727T091520_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNT_20180727T083305_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNS_20180727T091520_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNQ_20180727T091520_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNT_20180727T083305_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RPQ_20180727T091520_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNQ_20180727T091520_DVI_BOA_B8A_B11_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPS_20180727T083305_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNS_20180727T083305_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPS_20180727T083305_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNS_20180727T091520_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPT_20180727T083305_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNP_20180727T091520_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNQ_20180727T091520_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNS_20180727T083305_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNQ_20180727T091520_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNP_20180727T091520_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNT_20180727T083305_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPT_20180727T083305_DVI_BOA_B12_B4_B3.tif
Getting bands [4, 13] for S2B_MSIL2A_20180727T062629_N0206_R077_T41RNQ_20180727T091520_DVI_BOA.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPT_20180727T083305_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNS_20180727T091520_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPS_20180727T083305_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNS_20180727T083305_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNP_20180727T091520_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNQ_20180727T091520_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNP_20180727T091520_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPS_20180727T083305_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNS_20180727T083305_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPS_20180727T091520_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNT_20180727T083305_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPS_20180727T083305_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPS_20180727T091520_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNS_20180727T083305_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNT_20180727T083305_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPS_20180727T091520_DVI_BOA_B11_B8A_B4.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPS_20180727T091520_DVI_BOA_B12_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNQ_20180727T091520_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNQ_20180727T091520_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNT_20180727T083305_DVI_BOA_B8_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPT_20180727T083305_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPS_20180727T083305_DVI_BOA_B8A_B4_B12.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RPQ_20180727T091520_DVI_BOA_B4_B3_B2.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41RNP_20180727T091520_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SNS_20180727T091520_DVI_BOA_B8A_B4_B3.tif
Getting composite bands S2B_MSIL2A_20180727T062629_N0206_R077_T41SPS_20180727T091520_DVI_BOA_B8A_B4_B3.tif

*Or*

  • Select from metadataframe a subframe of Vegetation Index results related to a single georeferenced tiles of interest and multiple datetimes

mydate_results = results[(results.apply(lambda row: row[‘startdate’] in pd.to_datetime(mydates[0:4]) , axis=1)) & (results.apply(lambda row: ref_tiles[0] in row[‘title’], axis=1))]

mydate_results.head()

  • Create a new GeoDataFrame with the selected metadata plus the data values of the selected bands (set ‘all’ in order to select all the bands) cropped over the aoi_wkt2 area

myGDF = GeoDataFrame()

aoi_wkt = aoi_wkt2

for row in mydate_results.iterrows():

myGDF = myGDF.append(get_GDF_with_datavalues(row[1], user, api_key, ['dvi'], True, bbox=aoi_wkt.bounds), ignore_index=True)

#myGDF = myGDF.append(get_GDF_with_datavalues(row[1], user, api_key, 'all', True, bbox=aoi_wkt.bounds), ignore_index=True)

In [20]:

myGDF.head()

Out[20]:
band data enclosure enddate geo_transform identifier projection self startdate title wkt xsize ysize
0 DVI_BOA_B8A_B11_B2 [[[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0... https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 (598410.0, 10.0, 0.0, 3502780.0, 0.0, -10.0) 01be9eec66cb994310e469e012a3a2bc670f92e8 PROJCS["WGS 84 / UTM zone 41N",GEOGCS["WGS 84"... https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_2... POLYGON((64.0648798642054 31.5308017399005,65.... 1344 1538
1 DVI_BOA_B8A_B11_B2 [[[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0... https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 (598410.0, 10.0, 0.0, 3502780.0, 0.0, -10.0) 03aa38a9b430af9f2b2b2cb47cae31788c0fd59c PROJCS["WGS 84 / UTM zone 41N",GEOGCS["WGS 84"... https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_2... POLYGON((63.9421873560077 30.6397709351317,64.... 1344 1538
2 DVI_BOA_B8A_B4_B12 [[[149.0, 165.0, 165.0, 171.0, 171.0, 161.0, 1... https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 (598410.0, 10.0, 0.0, 3502780.0, 0.0, -10.0) 0949c278e99a1b2244ddd0edb3be4e38c7b10e81 PROJCS["WGS 84 / UTM zone 41N",GEOGCS["WGS 84"... https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41SNR_2... POLYGON((63.9422209286482 31.5415371135011,64.... 1344 1538
3 DVI_BOA_B8A_B4_B12 [[[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0... https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 (598410.0, 10.0, 0.0, 3502780.0, 0.0, -10.0) 0cccedf5d105c98b1bce071eebd3639084ebc43d PROJCS["WGS 84 / UTM zone 41N",GEOGCS["WGS 84"... https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_2... POLYGON((63.9421873560077 30.6397709351317,64.... 1344 1538
4 DVI_BOA_B11_B8A_B4 [[[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0... https://store.terradue.com/better-satcen-00004... 2018-07-29 06:16:31.024 (598410.0, 10.0, 0.0, 3502780.0, 0.0, -10.0) 13f7b9c57f80b8239602fa8643e850e2810dbe86 PROJCS["WGS 84 / UTM zone 41N",GEOGCS["WGS 84"... https://catalog.terradue.com//better-satcen-00... 2018-07-29 06:16:31.024 S2A_MSIL2A_20180729T061631_N0206_R034_T41RNQ_2... POLYGON((63.9421873560077 30.6397709351317,64.... 1344 1538

Plotting data and exporting geotif

Plot data

Choose DVI_BOA bands to plot: we are selecting the DVI_BOA data related to the first two reference tiles

In [21]:

bands = list()

mybands = myGDF[myGDF.apply(lambda row: any(ref_tile in row['title'] for ref_tile in ref_tiles[0:2]), axis=1) &
                    myGDF.apply(lambda row: 'DVI_BOA.tif' in row['title'], axis=1)]

for i,row in enumerate(mybands.iterrows()):
    title = '%s (band %s)'%(row[1]['title'],row[1]['band'])
    bmin = row[1]['data'].min()
    bmax = row[1]['data'].max()
    b = dict()
    if bmin != bmax:
        b[title] = (row[1]['data'] - bmin)/(bmax - bmin) * 255
    else:
        b[title] = row[1]['data']
    bands.append(b)

numbands = len(bands)
print numbands


8
  • Plot the selected bands separately
In [22]:

fig = plt.figure(figsize=(20,20))


m = len(b_list)
n = int(numbands/m)

j=0
for i,x in enumerate(bands):
    for k,v in x.iteritems():
        j=j+1
        a = fig.add_subplot(n, m, j)
        imgplot = plt.imshow(v.astype(np.uint8),
                             cmap='gray')
        a.set_title(k)

plt.tight_layout()
fig = plt.gcf()
plt.show()

fig.clf()
plt.close()
../../../../_images/pipelines_SATCEN_satcen-01-03-01_exploitation_satcen-00004_exploitation_from_pickle_48_0.png

Export a Geotiff

Export the band4 and dvi band related to the date time mydates[0] and ref_tiles[0,1] in a single Geotif

In [24]:

bands_number = len(b_list)

for i in range(2):
    A = mybands[mybands.apply(lambda row: ref_tiles[i] in row['title'], axis=1) &
                   mybands.apply(lambda row:  pd.to_datetime(mydates[0]) == row['startdate'], axis=1)]

    cols = A['xsize'].values[0]
    rows = A['ysize'].values[0]

    geo_transform = A['geo_transform'].values[0]
    projection = A['projection'].values[0]

    drv = gdal.GetDriverByName('GTiff')
    tiff_name = A['title'].values[0].split('.')[0]+'_'+'_'.join(b_list)+'.tif'

    ds = drv.Create(tiff_name, cols, rows, bands_number, gdal.GDT_Float32)

    ds.SetGeoTransform(geo_transform)
    ds.SetProjection(projection)

    ds.GetRasterBand(1).WriteArray(A['data'].values[0], 0, 0)
    ds.GetRasterBand(2).WriteArray(A['data'].values[1], 0, 0)
    ds.FlushCache()
  • Select and Plot the composites bands
In [25]:

mybands_composistes = myGDF[myGDF.apply(lambda row: pd.to_datetime(mydates[0]) == row['startdate'], axis=1) &
                            myGDF.apply(lambda row: ref_tiles[0] in row['title'], axis=1) &
                            myGDF.apply(lambda row: 'DVI_BOA.tif' not in row['title'], axis=1)]

In [26]:

fig = plt.figure(figsize=(20,20))

title = '_'.join(mybands_composistes['title'][0].split('_')[0:7])

fig.suptitle('%s : DVI composites bands' %title, fontsize=16)

for index,row in enumerate(mybands_composistes.iterrows()):
    index=index+1
    rgb_uint8 = np.dstack(row[1]['data']).astype(np.uint8)
    sub_fig = fig.add_subplot(2, 4, index)
    img = Image.fromarray(rgb_uint8)
    imgplot = plt.imshow(img)
    sub_fig.set_title(row[1]['band'])
../../../../_images/pipelines_SATCEN_satcen-01-03-01_exploitation_satcen-00004_exploitation_from_pickle_53_0.png

Download functionalities

Download a product

  • Define the download function
In [27]:

def get_product(url, dest, api_key):

    request_headers = {'X-JFrog-Art-Api': api_key}

    r = requests.get(url, headers=request_headers)

    open(dest, 'wb').write(r.content)

    return r.status_code
  • Get the reference download endpoint for the product related to a selected date (second date)
In [28]:

enclosure = myGDF[(myGDF['startdate'] == mydates[0])]['enclosure'].values[0]

enclosure

Out[28]:

'https://store.terradue.com/better-satcen-00004/_results/workflows/ec_better_ewf_satcen_01_03_01_ewf_satcen_01_03_01_0_7/run/26eddb46-20a5-11e9-93f5-0242ac11000f/0002389-181221095105003-oozie-oozi-W/c0f7cfbd-60e6-4842-8c2b-2c8ec3211f61/S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_20180729T082519_DVI_BOA_B8A_B11_B2.tif'

In [29]:

output_name = myGDF[(myGDF['startdate'] == mydates[0])]['title'].values[0]

output_name

Out[29]:

'S2A_MSIL2A_20180729T061631_N0206_R034_T41SPR_20180729T082519_DVI_BOA_B8A_B11_B2.tif'

In [30]:

get_product(enclosure,
            output_name,
            api_key)

Out[30]:

200

Bulk Download

  • Define the bulk download function
In [31]:

def get_product_bulk(row, api_key):

    return get_product(row['enclosure'],
                       row['title'],
                       api_key)
  • Download all the products related to a chosen date (first date) and a chosen reference tile
In [32]:

date_results = myGDF[(myGDF['startdate'] == mydates[0]) &
                     (myGDF.apply(lambda row: ref_tiles[0] in row['title'], axis=1))]


date_results.apply(lambda row: get_product_bulk(row, api_key), axis=1)

Out[32]:

0     200
12    200
14    200
32    200
38    200
39    200
41    200
45    200
64    200
dtype: int64