Understand Georelation

[1]:
import lxml.etree as etree
import requests
import cioppy

import numpy as np
import ipyleaflet
from shapely.geometry import box
from shapely.geometry import multipolygon
from shapely.wkt import loads
from ipyleaflet import Map, Polygon

ciop = cioppy.Cioppy()
[2]:
def get_params(osd):

    oss_ns = {'a':'http://www.w3.org/2001/XMLSchema',
          'b':'http://www.w3.org/2001/XMLSchema-instance',
          'c':'http://a9.com/-/opensearch/extensions/time/1.0/',
          'd':'http://www.opengis.net/eop/2.0',
          'e':'http://purl.org/dc/terms/',
          'f':'http://a9.com/-/spec/opensearch/extensions/parameters/1.0/',
          'g':'http://purl.org/dc/elements/1.1/',
          'h':'http://www.terradue.com/opensearch',
          'i':'http://a9.com/-/opensearch/extensions/geo/1.0/',
          'j':'http://a9.com/-/spec/opensearch/1.1/'}

    oss_content = etree.fromstring(requests.get(osd).content)

    url_template_element = oss_content.xpath('/j:OpenSearchDescription/j:Url[@type="application/atom+xml"]',
                                                 namespaces=oss_ns)[0]

    parameters = dict()

    for index, parameter in enumerate(url_template_element.xpath('.//f:Parameter', namespaces=oss_ns)):

        parameters[parameter.attrib['name']] = {'title' : parameter.attrib['title'],
                                                'value' : parameter.attrib['value']}

        options = []
        for option in parameter.xpath('.//f:Option', namespaces=oss_ns):

            options.append(option.attrib['value'])

        parameters[parameter.attrib['name']] = {'title' : parameter.attrib['title'],
                                                'value' : parameter.attrib['value'],
                                                'options' : options}
    return parameters
[3]:
def get_param_value(osd, os_parameter):

    params = get_params(osd)

    res = None

    for index, param in enumerate(params):

        if params[param]['value'] == os_parameter:

            res = params[param]
            res['name'] = param

    return res

Define the Sentinel-1 endpoint

[4]:
s1_osd_url = 'https://catalog.terradue.com/sentinel1/description'

Get the Sentinel-1 search parameters

[5]:
s1_parameters = get_params(s1_osd_url)

Print all the OpenSearch parameters

[6]:
for key, value in s1_parameters.iteritems():

    print '%24s' % s1_parameters[key]['value'],  s1_parameters[key]['title']
              {geo:uid?} The identifier of the resource within the search engine context (local reference)
    {t2:downloadOrigin?} a string that identifies the download origin (keyword, hostname...) to adapt the enclosure. If the parameter is enclosed between [] (e.g. [terradue]), enclosure will be returned only if there is a enclosure found for this source.
       {eop:cloudCover?} A number, set or interval requesting the cloud coverage
            {startIndex} index of the first search result desired
   {t2:vendorSpecifics?} A number, set or interval filtering vendor specific name:value
             {startPage} page number of the set of search results desired
       {eop:sensorType?} A string identifying the sensor type
     {eop:accessedFrom?} A string identifying the location from which the resource will be accessed. The catalogue shall return the download location in the enclosure atom link according to the parameter value.
         {t2:landCover?} A number, set or interval requesting the land coverage
      {eop:productType?} A string identifying the product type
           {time:start?} start of the temporal interval (RFC-3339)
 {eop:sensorResolution?} A string identifying the sensor spectral range
           {sct:source?} The URI of a source link (rel=via in ATOM).
  {eop:swathIdentifier?} Swath identifier that corresponds to precise incidence angles for the sensor
         {geo:geometry?} Geometry in WKT
         {geo:relation?} Spatial relation (possible values are “intersects”, “contains”, “disjoint”). The default is intersects.
 {eop:parentIdentifier?} A string identifying the collection of the entry in a hierarchy of dataset
  {eop:processingLevel?} A string identifying the processing level applied to the dataset
              {language} desired language of the results
{eop:platformSerialIdentifier?} A string with the Platform serial identifier
         {eop:platform?} A string with the platform short name
            {eop:track?} A number, set or interval requesting the range of orbit tracks
{t2:doubleCheckGeomtry?} Set to apply a finer geometry filtering
             {time:end?} stop of the temporal interval (RFC-3339)
         {dct:modified?} date after which dataset are updated (RFC-3339)
              {geo:box?} Rectangular bounding box
                 {count} number of search results per page desired
        {time:relation?} Temporal relation (possible values are “intersects”, “contains”, “during”, “disjoint”, “equals”)
   {eop:orbitDirection?} A string identifying the orbit direction
          {dct:subject?} The identifier of a category. Recommended best practice is to use a controlled vocabulary.
           {searchTerms} EO Free Text Search
         {t2:extension?} The value:name of a extension.
       {eop:instrument?} A string identifying the instrument
            {eop:title?} A name given to the resource
        {eop:orbitType?} A string identifying the orbit type

Find the parameter associated to ‘{eop:productType?}’ and its options

[7]:
geo_relation_parameter = get_param_value(s1_osd_url, '{geo:relation?}')
[8]:
geo_relation_parameter['name']
[8]:
'rel'
[9]:
options = s1_parameters[geo_relation_parameter['name']]['options']
[10]:
options
[10]:
['intersects', 'contains', 'disjoint']

Define an area of interest

[11]:
aoi_wkt = 'POLYGON ((20.9113 39.4866, 20.9113 40.0866, 20.3113 40.0866, 20.3113 39.4866, 20.9113 39.4866))'

Find the {geo:geometry?} parameter name

[12]:
wkt_parameter = get_param_value(s1_osd_url, '{geo:geometry?}')
[13]:
wkt_parameter
[13]:
{'name': 'geom',
 'options': [],
 'title': 'Geometry in WKT',
 'value': '{geo:geometry?}'}

Relation intersects

[14]:
search_params = dict([(geo_relation_parameter['name'],
                       options[0]),
                    (wkt_parameter['name'], aoi_wkt)])
[15]:
search_params
[15]:
{'geom': 'POLYGON ((20.9113 39.4866, 20.9113 40.0866, 20.3113 40.0866, 20.3113 39.4866, 20.9113 39.4866))',
 'rel': 'intersects'}
[16]:
search = ciop.search(end_point=s1_osd_url,
                     params=search_params,
                     output_fields='self,productType,track,enclosure,identifier,wkt,startdate',
                     model='EOP')

Show the first result returned

[17]:
m = Map(center=(39, 24), zoom=6)

aoi = Polygon(
    locations=np.asarray([t[::-1] for t in list(loads(aoi_wkt).exterior.coords)]).tolist(),
    color="green",
    fill_color="green",
    weight=2
)

for index, elem in enumerate(search):

    m += Polygon(locations=np.asarray([t[::-1] for t in list(loads(elem['wkt']).exterior.coords)]).tolist(), color="red", fill_color="red", weight=1, fill_opacity=0.1)

m.add_layer(aoi);

m

Relation contains

[18]:
buffer_size = 4

extended_aoi_wkt = box(*loads(aoi_wkt).buffer(buffer_size).bounds).wkt

[19]:
search_params = dict([(geo_relation_parameter['name'],
                       options[1]),
                     (wkt_parameter['name'], extended_aoi_wkt)])
[20]:
search_params
[20]:
{'geom': 'POLYGON ((24.9113 35.4866, 24.9113 44.0866, 16.3113 44.0866, 16.3113 35.4866, 24.9113 35.4866))',
 'rel': 'contains'}
[21]:
search = ciop.search(end_point=s1_osd_url,
                     params=search_params,
                     output_fields='self,productType,track,enclosure,identifier,wkt,startdate',
                     model='EOP')
[22]:
m = Map(center=(39, 24), zoom=5)

aoi = Polygon(
    locations=np.asarray([t[::-1] for t in list(loads(extended_aoi_wkt).exterior.coords)]).tolist(),
    color="green",
    fill_color="green",
    weight=2
)

for index, elem in enumerate(search):

    m += Polygon(locations=np.asarray([t[::-1] for t in list(loads(elem['wkt']).exterior.coords)]).tolist(), color="red", fill_color="red", weight=1, fill_opacity=0.1)

m.add_layer(aoi);

m

Relation disjoint

[23]:
buffer_size = 15

extended_aoi_wkt = box(*loads(aoi_wkt).buffer(buffer_size).bounds).wkt

[24]:
search_params = dict([(geo_relation_parameter['name'],
                       options[2]),
                     (wkt_parameter['name'], extended_aoi_wkt)])
[25]:
search_params
[25]:
{'geom': 'POLYGON ((35.9113 24.4866, 35.9113 55.0866, 5.311299999999999 55.0866, 5.311299999999999 24.4866, 35.9113 24.4866))',
 'rel': 'disjoint'}
[26]:
search = ciop.search(end_point=s1_osd_url,
                     params=search_params,
                     output_fields='self,productType,track,enclosure,identifier,wkt,startdate',
                     model='EOP')
[27]:
m = Map(center=(39, 24), zoom=1)

aoi = Polygon(
    locations=np.asarray([t[::-1] for t in list(loads(extended_aoi_wkt).exterior.coords)]).tolist(),
    color="green",
    fill_color="green",
    weight=2
)

for index, elem in enumerate(search):

    if not(isinstance(loads(elem['wkt']), multipolygon.MultiPolygon)):
        m += Polygon(locations=np.asarray([t[::-1] for t in list(loads(elem['wkt']).exterior.coords)]).tolist(), color="red", fill_color="red", weight=1, fill_opacity=0.1)

m.add_layer(aoi);

m