0

I'm trying to obtain satellite data for given coordinates by AWS, following tutorial here: https://sentinelhub-py.readthedocs.io/en/latest/examples/aws_request.html

What I'm getting is multiple tiles, not all of which cover the area of my interest, and I have the two following questions:

  1. can I somewhere specify (using sentinelhub&WebFeatureService that I only need tiles that cover the whole area?
  2. The search returns me coordinates as well, but I wonder what those coordinates really tell me, eg:
S2B_OPER_MSI_L1C_TL_MPS__20170812T112511_A002264_T34UGU_N02.05
coordinates info: [[[[23.71882836395446, 48.72091839211982], [25.208270245719866, 48.676085970164756], [25.128248164103578, 47.69086051716279], [23.66702320446364, 47.73417667514742], [23.71882836395446, 48.72091839211982]]]] 

Here we have 5 pairs of coordinates - what are those? If there were four pairs, I would understant it as coordinates of the corners, but there is more, sometimes even 11:

S2B_OPER_MSI_L1C_TL_MTI__20170815T113408_A002307_T35ULP_N02.05
coordinates info: [[[[24.280655565325187, 48.72090620487599], [25.18477006495381, 48.73870515546696], [25.074163489640565, 48.500333285001794], [25.007053685251048, 48.35557424560157], [24.939654082919624, 48.210945750598185], [24.873506894376266, 48.066054171483124], [24.8065216723558, 47.92142438867608], [24.742086017207008, 47.776261805908064], [24.727020474010263, 47.74266252879528], [24.332470529432726, 47.734164900431786], [24.280655565325187, 48.72090620487599]]]] 

My code looks like this:

def get_aws_data(bbox, config):
    search_bbox = BBox(bbox=bbox, crs=CRS.WGS84)
    search_time_interval = ('2017-07-01T00:00:00', '2017-09-30T23:59:59')

    wfs_iterator = WebFeatureService(
        search_bbox,
        search_time_interval,
        data_collection=DataCollection.SENTINEL2_L1C,
        maxcc=0.1,
        config=config
    )

    ids = []
    coords = []
    for tile_info in wfs_iterator:
        print("\n New tile: \n")
        print("coordinates info: {} \n".format(tile_info['geometry']['coordinates']))
        print("id info: {} \n".format(tile_info['properties']['id']))
        print("cloud info: {} \n".format(tile_info['properties']['cloudCoverPercentage']))
        ids.append(tile_info['properties']['id'])
        coords.append("coordinates info: {} \n".format(tile_info['geometry']['coordinates']))
        
    return ids, coords

1 Answer 1

1

Disclaimer: I work at Sinergise (company running the SentinelHub services and developing Python tools you mention) You might also have some more luck getting answers on a dedicated forum.

The WFS returns the Sentinel-2 tiles that intersect (have non empty intersection) with your bounding box. If you want to get just the tiles that completely cover your bounding box you will have to do another sweep over the results, compute intersections yourself and filter out the ones where the intersecting area is equal to your bounding box.

You can do that by converting the results from WFS to geodataframe with something like

import geopandas as gpd

geo_df = gpd.GeoDataFrame.from_features(list(wfs_iterator), crs='epsg:4326')

and then look for full intersections with covers or similar:

gdf.geometry.apply(lambda g: g.covers(search_bbox.geometry))

Once you have your responses in geodataframe, you can also visualise the results (tiles) rather quickly:

gdf.plot(alpha=0.2)

As a side note: 5 pairs of coordinates are needed to construct (closed) polygon (first and last pair are equal). More than that just means the tile geometry is more complex (have more vertices), or can even be multipolygon.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.