better-satcen-00002 data pipeline results (Sentinel-1 Multi-temporal SLC stack and Multi-temporal Coherence stack data pipeline): loading from pickle file and exploiting

This Notebook shows how to: * reload bbetter-satcen-00002 data pipeline results previously stored in a pickle file * Create a new dataframe with data values cropped over the area of interest * Plot an RGB * Create a geotif

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-00002_validation_NW_Niger.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-00002... 2018-09-22 05:19:19.147 0ecaac621b804f8362ca89b13025f218424aed16 https://catalog.terradue.com//better-satcen-00... 2018-09-22 05:19:03.036 S1B_S1_SLC__1SDV_20180910T051902_20180910T0519... POLYGON((10.6534444444444 12.1762222222222
1 https://store.terradue.com/better-satcen-00002... 2018-09-22 05:19:19.147 1c04a6d1f769d6d6c1842b5f833e2bcf83a9e9bf https://catalog.terradue.com//better-satcen-00... 2018-09-22 05:19:03.036 SLC_S1B_S1_SLC__1SDV_20180910T051902_20180910T... POLYGON((10.6534444444444 12.1762222222222
2 https://store.terradue.com/better-satcen-00002... 2018-09-22 05:19:19.147 31bba41907a285824f472b2569cd2ecfdda7c5a7 https://catalog.terradue.com//better-satcen-00... 2018-09-22 05:19:03.036 MTC_S1B_S1_SLC__1SDV_20180910T051902_20180910T... POLYGON((10.6534444444444 12.1762222222222
3 https://store.terradue.com/better-satcen-00002... 2018-09-22 05:19:19.147 519f895505caadd1d45a156b38627ef3b2260027 https://catalog.terradue.com//better-satcen-00... 2018-09-22 05:19:03.036 SLC_S1B_S1_SLC__1SDV_20180910T051902_20180910T... POLYGON((10.6534444444444 12.1762222222222
4 https://store.terradue.com/better-satcen-00002... 2018-09-22 05:19:19.147 5d629460a65b06ea9818970d0d1af340d19bbd54 https://catalog.terradue.com//better-satcen-00... 2018-09-22 05:19:03.036 SLC_S1B_S1_SLC__1SDV_20180910T051902_20180910T... POLYGON((10.6534444444444 12.1762222222222

Credentials for ellip platform access

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

import getpass

user = 'msagona'#getpass.getpass("Enter you ellip username:")
api_key = 'AKCp2V6wov45FKVuto1g558a3QKv5ZX4aD4vieCoZP8PtBU2TSFF2pBx57grLHp6zJuTY67S1'#getpass.getpass("Enter the API key:")

Define filter parameters

Time of interest

In [6]:

#mydates = pd.to_datetime(['2018-10-01 17:22:23.7090000', '2018-09-22 05:19:03.036000', '2018-09-10 05:19:02.7223870Z', '2018-08-29 05:19:02.107153', '2018-07-12 05:18:59.285207', '2018-07-09 17:22:19.676794', '2018-06-27 17:22:18.9696100Z', '2018-06-18 05:18:58.0635280Z', '2018-06-06 05:18:57.166982'])
#'2018-09-22 05:19:03.036000', '2018-08-29 05:19:02.107153',
mydates = pd.to_datetime(['2018-06-27 17:22:18.9696100Z', '2018-06-18 05:18:58.0635280Z', '2018-06-06 05:18:57.166982'])

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(9.838, 13.87)

In [8]:

buffer_size = 0.5

In [9]:

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

In [10]:

aoi_wkt1.wkt

Out[10]:

'POLYGON ((10.338 13.37, 10.338 14.37, 9.337999999999999 14.37, 9.337999999999999 13.37, 10.338 13.37))'
  • Or use the reference AOI (in this case the NW Niger area)
In [11]:

geom = 'MULTIPOLYGON (((10.65344444444444 12.17622222222222, 10.65344444444444 15.48261111111111, 8.398666666666667 15.48261111111111, 8.398666666666667 12.17622222222222, 10.65344444444444 12.17622222222222)), ((15.42594444444445 18.89533333333333, 15.42594444444445 23.45238888888889, 12.93755555555556 23.45238888888889, 12.93755555555556 18.89533333333333, 15.42594444444445 18.89533333333333)), ((6.452166666666667 18.3855, 6.452166666666667 20.24616666666667, 5.364944444444444 20.24616666666667, 5.364944444444444 18.3855, 6.452166666666667 18.3855)))'

aoi_wkt2 = loads(geom)[0]
#aoi_wkt2 = loads('POLYGON((9.267 14.695,10.392 14.69,10.393 13.186,9.27 13.188,9.267 14.695))')

In [12]:

aoi_wkt2.wkt

Out[12]:

'POLYGON ((10.65344444444444 12.17622222222222, 10.65344444444444 15.48261111111111, 8.398666666666667 15.48261111111111, 8.398666666666667 12.17622222222222, 10.65344444444444 12.17622222222222))'

Create a new dataframe with data values

  • Get the subframe metadata values related to selected bands and add the related data extended info: (cropped) data values and band name.

Define auxiliary methods to create new dataframe

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

def get_cropped_image_as_array(url, crop, bbox):

    #[-projwin ulx uly lrx lry]
    ulx, uly, lrx, lry = bbox[0], bbox[3], bbox[2], bbox[1]

    output = '/vsimem/clip.tif'
    gdal.UseExceptions()

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

    if crop==True:
        ds = gdal.Translate(output, ds, projWin = [ulx, uly, lrx, lry], projWinSRS = 'EPSG:4326')
        ds = None
    else:
        ds = gdal.Translate(output, ds)
        ds=None

    ds = gdal.Open('/vsimem/clip.tif')
    band = ds.GetRasterBand(1)

    w = band.XSize
    h = band.YSize

    data = band.ReadAsArray(0, 0, w, h).astype(np.float)

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

    ds = None
    band = None

    return data
  • The second one selects the data to be managed (only the actual results related to the filtered metadata) and returns a new GeoDataFrame containing the requested bands extended info with the original metadata set
In [14]:

def get_info(row, user, api_key, crop=True, bbox=[]):

    extended_info = dict()

    # add the band info
    extended_info['band'] = '_'.join(row['title'].rsplit('_',4)[1:])
    print 'Managing %s - band %s' %(row['title'].rsplit('_',4)[0], extended_info['band'])

    # add the data for the AOI
    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['aoi_data'] = get_cropped_image_as_array(url, crop, bbox)

    # create the series
    series = pd.Series(extended_info)

    return series
  • Select from metadataframe a subframe containing the coherence geotiff results related to the selected datetimes of interest
In [15]:

#mydate_results = results[((results.apply(lambda row: 'MTC_' in row['title'], axis=1)) |
#                          (results.apply(lambda row: 'SLC_' in row['title'], axis=1))) &
#                          (results.apply(lambda row: '.tif' in row['enclosure'], axis=1)) &
#                          (results.apply(lambda row: row['startdate'] == mydates[-1] , axis=1))]

mydate_results = results[(results.apply(lambda row: 'coh_' in row['title'], axis=1)) &
                         (results.apply(lambda row: row['startdate'] in mydates , axis=1))]


In [16]:

mydate_results

Out[16]:
enclosure enddate identifier self startdate title wkt
33 https://store.terradue.com/better-satcen-00002... 2018-06-27 17:22:35.073582 b18961af52c6fad3213feb5bb9d54094e4af230c https://catalog.terradue.com//better-satcen-00... 2018-06-27 17:22:18.969610 MTC_S1B_S1_SLC__1SDV_20180615T172218_20180615T... POLYGON((10.6534444444444 12.1762222222222
44 https://store.terradue.com/better-satcen-00002... 2018-06-18 05:19:14.166966 bce1a34e4a7a0f4feb35f468c96ea81c50c07964 https://catalog.terradue.com//better-satcen-00... 2018-06-18 05:18:58.063528 MTC_S1B_S1_SLC__1SDV_20180606T051857_20180606T... POLYGON((10.6534444444444 12.1762222222222
56 https://store.terradue.com/better-satcen-00002... 2018-06-06 05:19:13.262935 668a1aeeae2043b829312a825f82c75db804fbac https://catalog.terradue.com//better-satcen-00... 2018-06-06 05:18:57.166982 MTC_S1B_S1_SLC__1SDV_20180525T051856_20180525T... POLYGON((10.6534444444444 12.1762222222222

*PARAMETERS*

  • crop = True|False

if crop=False the bbox can be omitted

In [17]:

aoi_wkt = aoi_wkt1
#aoi_wkt = aoi_wkt2

mydate_results = mydate_results.merge(mydate_results.apply(lambda row: get_info(row, user, api_key, True, list(aoi_wkt.bounds)), axis=1), left_index=True, right_index=True)


Managing MTC_S1B_S1_SLC__1SDV_20180615T172218_20180615T172234_011387_014E95_4DF7_S1B_S1_SLC__1SDV_20180627T172218_20180627T172235_011562_01540A_C8D5 - band coh_VV_15Jun2018_27Jun2018
Managing MTC_S1B_S1_SLC__1SDV_20180606T051857_20180606T051913_011248_014A49_D5FE_S1B_S1_SLC__1SDV_20180618T051858_20180618T051914_011423_014FAE_9A00 - band coh_VV_06Jun2018_18Jun2018
Managing MTC_S1B_S1_SLC__1SDV_20180525T051856_20180525T051912_011073_0144B9_1A53_S1B_S1_SLC__1SDV_20180606T051857_20180606T051913_011248_014A49_D5FE - band coh_VV_25May2018_06Jun2018

In [18]:

mydate_results

Out[18]:
enclosure enddate identifier self startdate title wkt aoi_data band
33 https://store.terradue.com/better-satcen-00002... 2018-06-27 17:22:35.073582 b18961af52c6fad3213feb5bb9d54094e4af230c https://catalog.terradue.com//better-satcen-00... 2018-06-27 17:22:18.969610 MTC_S1B_S1_SLC__1SDV_20180615T172218_20180615T... POLYGON((10.6534444444444 12.1762222222222 [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... coh_VV_15Jun2018_27Jun2018
44 https://store.terradue.com/better-satcen-00002... 2018-06-18 05:19:14.166966 bce1a34e4a7a0f4feb35f468c96ea81c50c07964 https://catalog.terradue.com//better-satcen-00... 2018-06-18 05:18:58.063528 MTC_S1B_S1_SLC__1SDV_20180606T051857_20180606T... POLYGON((10.6534444444444 12.1762222222222 [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... coh_VV_06Jun2018_18Jun2018
56 https://store.terradue.com/better-satcen-00002... 2018-06-06 05:19:13.262935 668a1aeeae2043b829312a825f82c75db804fbac https://catalog.terradue.com//better-satcen-00... 2018-06-06 05:18:57.166982 MTC_S1B_S1_SLC__1SDV_20180525T051856_20180525T... POLYGON((10.6534444444444 12.1762222222222 [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,... coh_VV_25May2018_06Jun2018
  • Show data bands name
In [19]:

print list(mydate_results.band)


['coh_VV_15Jun2018_27Jun2018', 'coh_VV_06Jun2018_18Jun2018', 'coh_VV_25May2018_06Jun2018']

Plotting data and exporting geotif

Plot data

Choose bands to be plotted

In [20]:

bands = None

bands =  mydate_results['aoi_data'].values
titles = mydate_results['band'].values

numbands = len(bands)

for i in range(numbands):
    bmin = bands[i].min()
    bmax = bands[i].max()
    if bmin != bmax:
        bands[i] = (bands[i] - bmin)/(bmax - bmin) * 255
    else:
        print 'WARNING bmin %s = bmax %s (%i)!' %(bmin, bmax,i)
  • Plot the bands separately
In [21]:

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


for i in range(numbands):
    a = fig.add_subplot(3, 4, i+1)
    #
    imgplot = plt.imshow(bands[i].astype(np.uint8),
                         cmap='gray')
    a.set_title(titles[i])

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

fig.clf()
plt.close()




MemoryErrorTraceback (most recent call last)
/opt/anaconda/lib/python2.7/site-packages/IPython/core/formatters.pyc in __call__(self, obj)
    305                 pass
    306             else:
--> 307                 return printer(obj)
    308             # Finally look for special method names
    309             method = get_real_method(obj, self.print_method)

/opt/anaconda/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in <lambda>(fig)
    225
    226     if 'png' in formats:
--> 227         png_formatter.for_type(Figure, lambda fig: print_figure(fig, 'png', **kwargs))
    228     if 'retina' in formats or 'png2x' in formats:
    229         png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

/opt/anaconda/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in print_figure(fig, fmt, bbox_inches, **kwargs)
    117
    118     bytes_io = BytesIO()
--> 119     fig.canvas.print_figure(bytes_io, **kw)
    120     data = bytes_io.getvalue()
    121     if fmt == 'svg':

/opt/anaconda/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
   2198                     orientation=orientation,
   2199                     dryrun=True,
-> 2200                     **kwargs)
   2201                 renderer = self.figure._cachedRenderer
   2202                 bbox_inches = self.figure.get_tightbbox(renderer)

/opt/anaconda/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in print_png(self, filename_or_obj, *args, **kwargs)
    543
    544     def print_png(self, filename_or_obj, *args, **kwargs):
--> 545         FigureCanvasAgg.draw(self)
    546         renderer = self.get_renderer()
    547         original_dpi = renderer.dpi

/opt/anaconda/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in draw(self)
    462
    463         try:
--> 464             self.figure.draw(self.renderer)
    465         finally:
    466             RendererAgg.lock.release()

/opt/anaconda/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65

/opt/anaconda/lib/python2.7/site-packages/matplotlib/figure.pyc in draw(self, renderer)
   1142
   1143             mimage._draw_list_compositing_images(
-> 1144                 renderer, self, dsu, self.suppressComposite)
   1145
   1146             renderer.close_group('figure')

/opt/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in _draw_list_compositing_images(renderer, parent, dsu, suppress_composite)
    137     if not_composite or not has_images:
    138         for zorder, a in dsu:
--> 139             a.draw(renderer)
    140     else:
    141         # Composite any adjacent images together

/opt/anaconda/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65

/opt/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe)
   2424             renderer.stop_rasterizing()
   2425
-> 2426         mimage._draw_list_compositing_images(renderer, self, dsu)
   2427
   2428         renderer.close_group('axes')

/opt/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in _draw_list_compositing_images(renderer, parent, dsu, suppress_composite)
    137     if not_composite or not has_images:
    138         for zorder, a in dsu:
--> 139             a.draw(renderer)
    140     else:
    141         # Composite any adjacent images together

/opt/anaconda/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
     61     def draw_wrapper(artist, renderer, *args, **kwargs):
     62         before(artist, renderer)
---> 63         draw(artist, renderer, *args, **kwargs)
     64         after(artist, renderer)
     65

/opt/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in draw(self, renderer, *args, **kwargs)
    541         else:
    542             im, l, b, trans = self.make_image(
--> 543                 renderer, renderer.get_image_magnification())
    544             if im is not None:
    545                 renderer.draw_image(gc, l, b, im)

/opt/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in make_image(self, renderer, magnification, unsampled)
    768         return self._make_image(
    769             self._A, bbox, transformed_bbox, self.axes.bbox, magnification,
--> 770             unsampled=unsampled)
    771
    772     def _check_unsampled_image(self, renderer):

/opt/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in _make_image(self, A, in_bbox, out_bbox, clip_bbox, magnification, unsampled, round_to_pixel_border)
    371                     # clipping to [0, 1] and we use out-of-bounds
    372                     # values to carry the over/under/bad information
--> 373                     rgba = np.empty((A.shape[0], A.shape[1], 4), dtype=A.dtype)
    374                     rgba[..., 0] = A  # normalized data
    375                     # this is to work around spurious warnings coming

MemoryError:

<matplotlib.figure.Figure at 0x7f0b00af57d0>
  • Select 3 bands and plot them as an RGB
In [25]:

rgbbands=(bands[0],bands[1],bands[2])

rgb_uint8 = np.dstack(rgbbands).astype(np.uint16)

width = 10
height = 10
plt.figure(figsize=(width, height))
img = Image.fromarray(rgb_uint8)
imgplot = plt.imshow(img)



TypeErrorTraceback (most recent call last)
<ipython-input-25-0d1a1fb357d6> in <module>()
      6 height = 10
      7 plt.figure(figsize=(width, height))
----> 8 img = Image.fromarray(rgb_uint8)
      9 imgplot = plt.imshow(img)

/opt/anaconda/lib/python2.7/site-packages/PIL/Image.pyc in fromarray(obj, mode)
   2166         except KeyError:
   2167             # print typekey
-> 2168             raise TypeError("Cannot handle this data type")
   2169     else:
   2170         rawmode = mode

TypeError: Cannot handle this data type

<matplotlib.figure.Figure at 0x7ff96bd22690>

Export the Geotiff

In [22]:

band_number = 3
cols = myGDF['xsize'].values[0]
rows = myGDF['ysize'].values[0]

print (cols,rows)

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

drv = gdal.GetDriverByName('GTiff')

ds = drv.Create('export.tif', cols, rows, band_number, gdal.GDT_Float32)

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

ds.GetRasterBand(1).WriteArray(rgbbands[0], 0, 0)
ds.GetRasterBand(2).WriteArray(rgbbands[1], 0, 0)
ds.GetRasterBand(3).WriteArray(rgbbands[2], 0, 0)
ds.FlushCache()

(1241, 1567)

Download functionalities

Download a product

  • Define the download function
In [23]:

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 the first date
In [24]:

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

enclosure

Out[24]:

'https://store.terradue.com/better-wfp-00005/_results/workflows/ec_better_ewf_wfp_01_02_03_ewf_wfp_01_02_03_0_5/run/0fe626d4-19bf-11e9-adb5-0242ac11000f/0001333-181221095105003-oozie-oozi-W/ff52d47e-a750-4368-a476-55a197fa1980/S1A_IW_GRDH_1SDV_20170929T132959_20170929T133024_018591_01F571_F9F4_Sigma0_all_bands.tif'

In [25]:

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

output_name

Out[25]:

'S1A_IW_GRDH_1SDV_20170929T132959_20170929T133024_018591_01F571_F9F4_Sigma0_all_bands.tif'

In [26]:

get_product(enclosure,
            output_name,
            api_key)

Out[26]:

200

Bulk Download

  • Define the bulk download function
In [27]:

def get_product_bulk(row, api_key):

    return get_product(row['enclosure'],
                       row['title'],
                       api_key)
  • Download all the products related to the chosen dates
In [28]:

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

Out[28]:

0     200
1     200
2     200
3     200
4     200
5     200
6     200
7     200
8     200
9     200
10    200
11    200
dtype: int64