Discover Sentinel-1 pairs for pre-seismic and coseismic coherence change analysis¶

Import the Python packages

[1]:

%matplotlib inline

import warnings
warnings.filterwarnings("ignore")
import os
import sys
import glob

import cioppy
ciop = cioppy.Cioppy()

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as colors

from snappy import jpy
from snappy import ProductIO
from snappy import GPF
from snappy import HashMap

import gc

from shapely.geometry import box
from shapely.wkt import loads

import py_earthquakes

from datetime import datetime, timedelta
import dateutil.parser

import geopandas as gp

import folium

Use an earthquake

[2]:

bbox=[21.944860, 36.940880,23.944860, 38.940880]

min_mag = 4.1

start_date = '2016-10-01'

stop_date = '2016-11-30'

[3]:

eq_search = py_earthquakes.EarthQuakes(start_date,
                                       stop_date,
                                       min_mag = min_mag,
                                       bbox=bbox)

[4]:

eq_search.earthquakes[0].title

[4]:

'M 4.2 - 1km SSW of Makrakomi, Greece'

[5]:

eq_search.earthquakes[0].id

[5]:

'us20007gi7'

[6]:

eq_search.earthquakes[0].date

[6]:

'2016-10-25T00:09:56.590000Z'

Create a buffer of 0.5 degrees

[7]:

eq_search.earthquakes[0].wkt

[7]:

'POINT(22.1115 38.9244)'

[8]:

aoi_wkt = box(*loads(eq_search.earthquakes[0].wkt).buffer(0.5).bounds).wkt

[9]:

aoi_wkt

[9]:

'POLYGON ((22.6115 38.4244, 22.6115 39.4244, 21.6115 39.4244, 21.6115 38.4244, 22.6115 38.4244))'

Search parameters

Set the catalogue endpoint to Sentinel-1:

[10]:

series = 'https://catalog.terradue.com/sentinel1/search'

Define the end of the time of interest:

[11]:

slave_search_stop_date = (dateutil.parser.parse(eq_search.earthquakes[0].date) + timedelta(days=6)).isoformat()

Build and submit the catalog search

[12]:

search_params = dict([('geom', aoi_wkt),
                     ('start', eq_search.earthquakes[0].date),
                     ('stop', slave_search_stop_date),
                      ('pt', 'SLC')])

[13]:

slave_search = ciop.search(end_point = series,
                     params = search_params,
                     output_fields='self,productType,track,enclosure,identifier,wkt,startdate',
                     model='EOP')

Put all slaves in a geodataframe and plot the Sentinel-1 slave candidates

[14]:

aoi = loads(aoi_wkt)

[15]:

result = []

locations = []

for index, elem in enumerate(slave_search):

    locations.append([t[::-1] for t in list(loads(elem['wkt']).exterior.coords)])

    slave_wkt = loads(elem['wkt'])

    result.append({'self' : elem['self'],
                   'identifier' : elem['identifier'],
                   'enclosure' : elem['enclosure'],
                   'date' : elem['startdate'],
                   'wkt': loads(elem['wkt']),
                   'aoi_intersec' : (slave_wkt.intersection(aoi).area/aoi.area) * 100,
                   'contains': slave_wkt.contains(aoi)
                  })

slaves = gp.GeoDataFrame(result)

[16]:

lat = loads(eq_search.earthquakes[0].wkt).y
lon = loads(eq_search.earthquakes[0].wkt).x

zoom_start = 7

m = folium.Map(location=[lat, lon], zoom_start=zoom_start)

radius = 4
folium.CircleMarker(
    location=[lat, lon],
    radius=radius,
    color='#FF0000',
    stroke=False,
    fill=True,
    fill_opacity=0.6,
    opacity=1,
    popup='{} pixels'.format(radius),
    tooltip='I am in pixels',
).add_to(m)


folium.PolyLine(
    locations=np.asarray([t[::-1] for t in list(loads(aoi_wkt).exterior.coords)]).tolist(),
    color='#FF0000',
    weight=2,
    tooltip='Japan flooding',
).add_to(m)

folium.PolyLine(
    locations=locations,
    color='orange',
    weight=1,
    opacity=1,
    smooth_factor=0,
).add_to(m)

m.save(os.path.join('results', '%s_search.html' % eq_search.earthquakes[0].id))

m

[16]:

[17]:

slaves

[17]:

	aoi_intersec	contains	date	enclosure	identifier	self	wkt
0	100.000000	True	2016-10-30T16:31:41.3695830Z	https://store.terradue.com/download/sentinel1/...	S1A_IW_SLC__1SDV_20161030T163141_20161030T1632...	https://catalog.terradue.com/sentinel1/search?...	POLYGON ((19.736723 39.341087, 22.701288 39.74...
1	97.711574	False	2016-10-30T04:38:48.2069410Z	https://store.terradue.com/download/sentinel1/...	S1B_IW_SLC__1SDV_20161030T043848_20161030T0439...	https://catalog.terradue.com/sentinel1/search?...	POLYGON ((22.34882 37.795399, 19.450127 38.198...
2	5.551913	False	2016-10-30T04:38:22.4261540Z	https://store.terradue.com/download/sentinel1/...	S1B_IW_SLC__1SDV_20161030T043822_20161030T0438...	https://catalog.terradue.com/sentinel1/search?...	POLYGON ((22.728542 39.285812, 19.762106 39.68...
3	71.287095	False	2016-10-25T16:23:37.9832800Z	https://store.terradue.com/download/sentinel1/...	S1A_IW_SLC__1SDV_20161025T162337_20161025T1624...	https://catalog.terradue.com/sentinel1/search?...	POLYGON ((21.63916 39.960037, 24.634281 40.361...
4	5.062699	False	2016-10-25T16:23:13.1686030Z	https://store.terradue.com/download/sentinel1/...	S1A_IW_SLC__1SDV_20161025T162313_20161025T1623...	https://catalog.terradue.com/sentinel1/search?...	POLYGON ((22.016951 38.469093, 24.94478 38.871...
5	93.929883	False	2016-10-25T04:30:36.4692880Z	https://store.terradue.com/download/sentinel1/...	S1B_IW_SLC__1SDV_20161025T043036_20161025T0431...	https://catalog.terradue.com/sentinel1/search?...	POLYGON ((24.384764 37.672104, 21.492056 38.07...

Select the slave that ‘best’ covers the AOI

[18]:

slave = slave_search[slaves['aoi_intersec'].idxmax()]

slave

[18]:

{'enclosure': 'https://store.terradue.com/download/sentinel1/files/v1/S1A_IW_SLC__1SDV_20161030T163141_20161030T163208_013722_01603F_4094',
 'identifier': 'S1A_IW_SLC__1SDV_20161030T163141_20161030T163208_013722_01603F_4094',
 'productType': 'SLC',
 'self': 'https://catalog.terradue.com/sentinel1/search?format=atom&uid=S1A_IW_SLC__1SDV_20161030T163141_20161030T163208_013722_01603F_4094',
 'startdate': '2016-10-30T16:31:41.3695830Z',
 'track': '175',
 'wkt': 'POLYGON((19.736723 39.341087,22.701288 39.742973,23.034481 38.123413,20.138891 37.720314,19.736723 39.341087))'}

[19]:

slave['startdate']

[19]:

'2016-10-30T16:31:41.3695830Z'

Search for the pre-event masters

[20]:

master_search_start_date = (dateutil.parser.parse(slave['startdate']) + timedelta(days=-24)).isoformat()

[21]:

master_search_stop_date = (dateutil.parser.parse(slave['startdate']) + timedelta(days=-1)).isoformat()

[22]:

master_search_params = dict([('geom', slave['wkt']),
                           ('track', slave['track']),
                            ('pt',slave['productType']),
                            ('start', master_search_start_date),
                            ('stop', master_search_stop_date)])

[23]:

try:
    master_search = ciop.search(end_point=series,
                            params=master_search_params,
                            output_fields='identifier,enclosure,self,startdate,wkt',
                            model='EOP')
except IndexError:
    print('no masters')

[24]:

result = []

for index, elem in enumerate(master_search):

    master_wkt = loads(elem['wkt'])

    result.append({'self' : elem['self'],
                   'identifier' : elem['identifier'],
                   'enclosure' : elem['enclosure'],
                   'wkt': loads(elem['wkt']),
                   'aoi_intersec' : (master_wkt.intersection(aoi).area/aoi.area) * 100,
                   'slave_intersec' : (master_wkt.intersection(loads(slave['wkt']))).area / loads(slave['wkt']).area * 100,
                   'contains': master_wkt.contains(aoi),
                   'days': (dateutil.parser.parse(slave['startdate']) - dateutil.parser.parse(elem['startdate'])).days
                  })

masters = gp.GeoDataFrame(result)

[25]:

masters

[25]:

	aoi_intersec	contains	days	enclosure	identifier	self	slave_intersec	wkt
0	68.636837	False	6	https://store.terradue.com/download/sentinel1/...	S1B_IW_SLC__1SDV_20161024T163111_20161024T1631...	https://catalog.terradue.com/sentinel1/search?...	55.755995	POLYGON ((19.520258 40.118217, 22.542038 40.52...
1	44.640286	False	6	https://store.terradue.com/download/sentinel1/...	S1B_IW_SLC__1SDV_20161024T163046_20161024T1631...	https://catalog.terradue.com/sentinel1/search?...	52.199313	POLYGON ((19.906658 38.567837, 22.860817 38.97...
2	0.000000	False	11	https://store.terradue.com/download/sentinel1/...	S1A_IW_SLC__1SDV_20161018T163206_20161018T1632...	https://catalog.terradue.com/sentinel1/search?...	7.893120	POLYGON ((19.368031 40.835575, 22.396507 41.23...
3	100.000000	True	12	https://store.terradue.com/download/sentinel1/...	S1A_IW_SLC__1SDV_20161018T163141_20161018T1632...	https://catalog.terradue.com/sentinel1/search?...	99.881697	POLYGON ((19.736567 39.342892, 22.70108 39.744...
4	0.000000	False	12	https://store.terradue.com/download/sentinel1/...	S1A_IW_SLC__1SDV_20161018T163116_20161018T1631...	https://catalog.terradue.com/sentinel1/search?...	7.974487	POLYGON ((20.105968 37.850845, 23.007288 38.25...
5	0.000000	False	23	https://store.terradue.com/download/sentinel1/...	S1A_IW_SLC__1SDV_20161006T163206_20161006T1632...	https://catalog.terradue.com/sentinel1/search?...	7.933108	POLYGON ((19.367752 40.834766, 22.396078 41.23...
6	100.000000	True	24	https://store.terradue.com/download/sentinel1/...	S1A_IW_SLC__1SDV_20161006T163141_20161006T1632...	https://catalog.terradue.com/sentinel1/search?...	99.917743	POLYGON ((19.736254 39.342201, 22.700668 39.74...
7	0.000000	False	24	https://store.terradue.com/download/sentinel1/...	S1A_IW_SLC__1SDV_20161006T163116_20161006T1631...	https://catalog.terradue.com/sentinel1/search?...	7.912009	POLYGON ((20.105749 37.849777, 23.00696 38.252...

Select the two masters according to the ranking of AOI coverage and nearest cycles in time

[26]:

master_1 = master_search[masters.sort_values(['aoi_intersec', 'days'], ascending=[False, False]).iloc[1].name]
master_2 = master_search[masters.sort_values(['aoi_intersec', 'days'], ascending=[False, False]).iloc[0].name]

[34]:

s1_identifiers = []
s1_references = []
locations = []

for product in [slave, master_1, master_2]:

    locations.append([t[::-1] for t in list(loads(product['wkt']).exterior.coords)])

    s1_identifiers.append(product['identifier'])
    s1_references.append(product['self'])

Plot the Sentinel-1 products (slave, master 1 and master 2), the earthquake point and its area of interest

[35]:

lat = loads(eq_search.earthquakes[0].wkt).y
lon = loads(eq_search.earthquakes[0].wkt).x

zoom_start = 7

m = folium.Map(location=[lat, lon], zoom_start=zoom_start)

radius = 4
folium.CircleMarker(
    location=[lat, lon],
    radius=radius,
    color='#FF0000',
    stroke=False,
    fill=True,
    fill_opacity=0.6,
    opacity=1,
    popup='{} pixels'.format(radius),
    tooltip='I am in pixels',
).add_to(m)


folium.PolyLine(
    locations=np.asarray([t[::-1] for t in list(loads(aoi_wkt).exterior.coords)]).tolist(),
    color='#FF0000',
    weight=2,
    tooltip='Japan flooding',
).add_to(m)

folium.PolyLine(
    locations=locations,
    color='orange',
    weight=1,
    opacity=1,
    smooth_factor=0,
).add_to(m)

m.save(os.path.join('results', '%s_final.html' % eq_search.earthquakes[0].id))

m

[35]:

[36]:

s1_references

[36]:

['https://catalog.terradue.com/sentinel1/search?format=atom&uid=S1A_IW_SLC__1SDV_20161030T163141_20161030T163208_013722_01603F_4094',
 'https://catalog.terradue.com/sentinel1/search?format=atom&uid=S1A_IW_SLC__1SDV_20161018T163141_20161018T163208_013547_015AEB_5994',
 'https://catalog.terradue.com/sentinel1/search?format=atom&uid=S1A_IW_SLC__1SDV_20161006T163141_20161006T163208_013372_01554F_F08D']

[37]:

pair_coseismic = [s1_references[0],
                 s1_references[1]]

pair_preseismic = [s1_references[1],
                 s1_references[2]]

[45]:

print ('aoi_wkt = \'%s\'' % aoi_wkt)

aoi_wkt = 'POLYGON ((22.6115 38.4244, 22.6115 39.4244, 21.6115 39.4244, 21.6115 38.4244, 22.6115 38.4244))'

[40]:

print ('pair_coseismic = %s' % pair_coseismic)

pair_coseismic = ['https://catalog.terradue.com/sentinel1/search?format=atom&uid=S1A_IW_SLC__1SDV_20161030T163141_20161030T163208_013722_01603F_4094', 'https://catalog.terradue.com/sentinel1/search?format=atom&uid=S1A_IW_SLC__1SDV_20161018T163141_20161018T163208_013547_015AEB_5994']

[41]:

print ('pair_preseismic = %s' % pair_preseismic)

pair_preseismic = ['https://catalog.terradue.com/sentinel1/search?format=atom&uid=S1A_IW_SLC__1SDV_20161018T163141_20161018T163208_013547_015AEB_5994', 'https://catalog.terradue.com/sentinel1/search?format=atom&uid=S1A_IW_SLC__1SDV_20161006T163141_20161006T163208_013372_01554F_F08D']

[ ]: