Python example for using sat-search

[vka012]

Hi

This is not an issue, but more of a question. Im fairly new to python, and looking for a neat and efficient way of searching for Landsat data. For now, I have been using the landsat-util, but if I understood correctly this is getting deprecated. I have therefore been directed to the sat-search library.

So, I was wondering if someone could take the time to provide me with an example or short tutorial on how to use this library for searching for satellite data in python (and not in command line). Im especially interested in the various input features in the search function (sensor, date, geometry, etc.)

Thanks, Vebjørn

[matthewhanson]

Hello @vka012 , thanks for checking out sat-search. It's true, we do need to provide documentation, but let me try to help in the meantime. The CLI code has some clues on how the Python library works: https://github.com/sat-utils/sat-search/blob/master/satsearch/main.py#L9

Notably, see the initial Search call:

search = Search(**kwargs)
scenes = Scenes(search.scenes(), properties=kwargs)

What this is doing is it's passing in kwargs to the Search, which performs the search. But what are kwargs you ask? kwargs is simply a dictionary of parameters that have been passed in from the command line, and the double asterisk is just converting the dictionary to a keywords. So these two calls are equivalent:

search = Search(**{"datetime": "2018-01-01"})
search = Search(datetime="2018-01-01")

Ok, but now you are probably wondering what the list of available parameters are that you can search on. The answer is that any field that appears in the properties of an Item or a Collection can be searched on:

Sentinel properties example:

properties: {
    id: "L1C_T41DMG_A005078_20180225T044720",
    c:id: "sentinel-2-l1c",
    datetime: "2018-02-25T04:47:20.460Z",
    eo:platform: "Sentinel-2B",
    eo:cloud_cover: 100,
    eo:epsg: "32741",
    sentinel:product_id: "S2B_MSIL1C_20180225T044719_N0206_R047_T41DMG_20180225T092707"
}

Sentinel collection example: https://sat-api.developmentseed.org/collections/sentinel-2-l1c/definition

properties: {
    c:id: "sentinel-2-l1c",
    c:name: "Sentinel 2 L1C",
    c:description: "Sentinel-2a and Sentinel-2b imagery",
    provider: "ESA",
    license: "https://sentinel.esa.int/documents/247904/690755/Sentinel_Data_Legal_Notice",
    eo:gsd: 10,
    eo:instrument: "MSI",
    eo:off_nadir: 0
}

Landsat properties example:

properties: {
    id: "LC08_L1TP_026009_20150729_20170406_01_T1",
    c:id: "landsat-8-l1",
    datetime: "2015-07-29T16:44:13.782Z",
    eo:cloud_cover: 89,
    eo:sun_azimuth: 177.64085481,
    eo:sun_elevation: 36.48367326,
    landsat:path: 26,
    landsat:row: 9
}

Landsat Collection example: https://sat-api.developmentseed.org/collections/landsat-8-l1/definition

properties: {
    c:id: "landsat-8-l1",
    c:name: "Landsat 8 L1",
    c:description: "Landat 8 imagery radiometrically calibrated and orthorectified using gound points 
        and Digital Elevation Model (DEM) data to correct relief displacement.",
    provider: "USGS",
    license: "PDDL-1.0",
    eo:gsd: 30,
    eo:platform: "landsat-8",
    eo:instrument: "OLI_TIRS",
    eo:off_nadir: 0
}

So any of these properties can be searched on. In the next version of STAC a more complex querying language will be allowed so we hope to get that into sat-api and sat-search soon.

Note that for properties with a colon in the name, such as eo:cloud_cover, you can't use those as keywords as it's not valid Python syntax:

# this will not work
search = Search(eo:cloud_cover="0/20")

# this will work
search = Search(**{"eo:cloud_cover": "0/20"})

Note that a range search, like 0/20 above, will search for cloud cover between 0 and 20%.

After getting the search results:

# create a scenes object, attaching the search query dictionary (optional) as properties
scenes = Scenes(search.scenes(), properties=kwargs)

# save the search results as a GeoJSON FeatureCollection
scenes.save(filename="mysearch.geojson")

Hope this helps, please don't hesitate to post more questions, it'll help us create better documentation.

Also note that things will change over the next month as the new STAC version comes out, watch out blog on Medium for announcements: https://medium.com/devseed

[vka012]

Hi again

Thank you so much for this detailed reply. I now understand more how the sat-search library can be used in Python. However, I still dont manage to use it exactly how I want. I would therefore kindly like to ask if you could write an example code from scratch filling my requirements (preferably including the imports you do), such as:

from satsearch import Scene, Search
...

I want to search for Landsat-8 scenes in the arctic region covering Svalbard. Is it possible to provide some sort of bounding box to the search function, say a geojson file defining this area of interest? Or is it only possible to define the search geometry by paths and rows? Further, the cloud cover percentage should preferably be below 30%, i.e., 0/30. The search should also include a starttime and endtime defining the period of time the search should cover.

From the possible scenes found in this search I would like to retrieve information such as the productid, sensingtime, and footprint (and possibly some url showing a preview of the satellite image).

Hope to hear from you, Vebjørn

[matthewhanson]

@vka012 Yes, you can search by area of interest.

from satsearch import Search, Scene
import json

filename = 'myaoi.geojson'

# load geojson from file
with open(filename) as f:
    geoj = json.dumps(json.loads(f.read()))

# search for matching scenes
results = Search(intersects=geoj, datetime="2018-01-01/2018-06-30", **{"eo:cloud_cover": "0/30"})

scenes = Scenes(results.scenes)

# save scenes as GeoJSON for viewing footprints in QGIS or similar
scenes.save("myscenes.geojson")

The myscenes.geojson file is a FeatureCollection with each scenes as a Feature. so if you want to print out the timestamp or scene ID you can just iterate through the Features. The variables are in the properties field in each Feature.

To take this one step further, if you want to download the preview images:

scenes.download('thumbnail')

This will download all of the thumbnail images. The directory which it saves them in depends on the SATUTILS_DATADIR and SATUTILS_FILENAME environment variables. Read more about using those in the sat-search README under Download options and Metadata Patterns.

[vka012]

@matthewhanson Thank you. This is exactly what I am looking for.

However, there is something wrong with the example. My guess is that the error occurs when defining the scenes variable.

scenes = Scenes(results.scenes)

The following error appears: TypeError: argument of type 'method' is not iterable.

Hope to hear from you, Vebjørn

[matthewhanson]

@vka012 ah sorry it should be

myscenes = Scenes(results.scenes())

since it's a function. I also renamed variable to myscenes so it wouldn't conflict with the Scenes class.

[vka012]

Thanks for being so helpful. It now works perfectly. Here is my example code:

from satsearch import Search, Scenes
import json

import logging
logging.basicConfig(level=logging.DEBUG)

filename = 'myaoi.geojson'

# load geojson from file
with open(filename) as f:
    geoj = json.dumps(json.loads(f.read()))

    # search for matching scenes
    results = Search(intersects=geoj, datetime="2018-06-01/2018-06-30", **{"eo:cloud_cover": "0/30"})

    scenes = Scenes(results.scenes())

    # save scenes as GeoJSON for viewing footprints in QGIS or similar
    scenes.save("myscenes.geojson")

The geoj is a string with value:

{"type": "FeatureCollection", "features": [{"type": "Feature", "properties": {}, "geometry": {"type": "Polygon", "coordinates": [[[-5.712890625, 75.95223506623554], [65.7421875, 75.95223506623554], [65.7421875, 82.04574006217713], [-5.712890625, 82.04574006217713], [-5.712890625, 75.95223506623554]]]}}]}

Best, Vebjørn

[vka012]

Hi @matthewhanson

Again, thanks for being so helpfull previously. Once again I need to ask you for some guidance. I am frequently using sat-search for finding Landsat-8 data. With that in mind, I wonder if there exists a way of downloading the full data set of one landsat scene? As for now, I have found that a Scene object have a download and download_file utility, but I only manage to download one band seperately with these.

Also, I am not able to search for landsat scenes by identifier.

results = Search(id="LC08_L1TP_026009_20150729_20170406_01_T1")
scn = results.scenes()

The output is empty. I am using the id of the landsat properties example you listed above.

Best, Vebjørn

[matthewhanson]

@vka012

I just published sat-search 0.2.0 which is using the newly deployed version of sat-api so you will want to update as the old version will no longer work with the deployed version of sat-api.developmentseed.org

I added a Jupyter notebook tutorial in the top level of this repo which should help with examples of how to use as a library. Note that the scenes() function is now items().

To search by ids you can now:

search = Search.search(ids=['id1', 'id2'], collection='landsat-8-l1')
items = search.items()

Note you need the collection name.

To download the full dataset (all asset keys) you can do this with the CLI now by specifying --download ALL but for the library you'll just need to loop through all the assets:

# get all keys from all items
keys = set([k for i in items for k in i.assets])

for key in keys:
    items.download(key=key)

I'm going to close this issue but please don't hesitate to start a new one if you have any more questions or issues!

[dempsey-ryan]

@vka012

I just published sat-search 0.2.0 which is using the newly deployed version of sat-api so you will want to update as the old version will no longer work with the deployed version of sat-api.developmentseed.org

I added a Jupyter notebook tutorial in the top level of this repo which should help with examples of how to use as a library. Note that the scenes() function is now items().

To search by ids you can now:

search = Search.search(ids=['id1', 'id2'], collection='landsat-8-l1')
items = search.items()

Note you need the collection name.

To download the full dataset (all asset keys) you can do this with the CLI now by specifying --download ALL but for the library you'll just need to loop through all the assets:

# get all keys from all items
keys = set([k for i in items for k in i.assets])

for key in keys:
    items.download(key=key)

I'm going to close this issue but please don't hesitate to start a new one if you have any more questions or issues!

Following this procedure, I get the warning Unable to download s3://sentinel-s2-l1c/tiles/18/T/VQ/2018/6/30/0/B04.jp2: No connection adapters were found for 's3://sentinel-s2-l1c/tiles/18/T/VQ/2018/6/30/0/B04.jp2' for every single jp2 file.

While I am here, I will ask: my goal is to obtain a satellite image (as an array of pixels) given a polygon. Is this the right package for that task? I will have a large number of polygons, probably close to 1M, each a square of about 1 metre x 1 metre. If I am only wanting to grab the satellite image data, what extra argument(s) would I pass to items.download? Thanks!

Upon inspection, it seems that items.download(key='thumbnail') is what I'm looking for. Does that seem right? I tried key='visual', key='visual_20m', and key='visual_60m' since those were in keys and the names sound like maybe they're what I want but it says no such asset.

And my last question, when using Search, what EPSG should the parameters be? For example, if I want a certain bbox, should that be in WGS84? Or another?