There are so many ways to download Sentinel 2 data. I am not going to list them all here, but most involve an interaction with a website that contains a map. This is fine for many users. Sometimes though it is preferable to access via the command line or a script and again there are tools that allow you to do this. One of these was recently updated – sentinelsat. I haven’t used this before and having talked about it on #scenefromabove podcast a couple of weeks ago, I thought I would give it a go.
In this post and accompanying notebook I’ll show you how with a scihub account you can build a script to search, download and unzip Sentinel2 level 2A data using Python. The only thing to note is that the boundary file you will search on needs to be projected in WGS 84 or which is EPSG code 4326. I am using a geojson file for my boundary as that follows the example in the documentation. I assume you can use other vector formats for a boundary, but if you cannot then use ogr2ogr to convert to geojson (or QGIS – its pretty good!).
Installation
If you run your Python environment via Anaconda (as I do), then install sentinelsat using pip install. Open the anaconda prompt and type
pip install sentinelsat
That is it. For any installation issues information have a look at the docs here.
Searching
In this example I have created a geojson file that covers the Isle of Wight – I’ve called it bounds_IOW.geojson. We can use GeoPandas to plot this file but in this example I am going to use folium to plot the bounding box over OpenStreetMap. Using the code below in a Jupyter Notebook:
m = folium.Map([50.6, -1.3], zoom_start=10)
boundsdata = r'D:\sentinelsat\bounds_IOW.geojson'
folium.GeoJson(boundsdata).add_to(m)
m
I am using the folium example from the documentation, so the variable ‘m’ is my map. I’ve centered it roughly on the Isle of Wight and zoomed in (zoom_start=10). Then I have added my bounds_IOW.geojson to the map (‘m’) and finally called m in the last line to plot it.
I really like folium, I have written a little about in the past, if you are interested it is here.
Now that we know our geojson file is the location we are expecting, lets use it with sentinelsat to discover if there is any imagery available. You’ll need a username and password to the scihub, but once you have these just paste them into the code below:
user = 'user' ## change this!
password = 'password' ## change this!
api = SentinelAPI(user, password, 'https://scihub.copernicus.eu/dhus')
footprint = geojson_to_wkt(read_geojson(boundsdata))
print (footprint)
Sentinelsat is looking for a boundary in wkt (well known text) format, so the penultimate line here coverts our boundary in geojson to wkt. Now we are ready to search.
products = api.query(footprint,
date = ('20190225', '20190227'),
platformname = 'Sentinel-2',
processinglevel = 'Level-2A',
cloudcoverpercentage = (0, 20))
Basically we parse a query to the api to search for products between a date, platform, the processing level and a cloud cover %. This will return an ordered dictionary. Based on the results you may wish to adjust the values in your query, or even add other parameters to the query. Use print(products) to inspect the ordered dictionary if unsure.
The particular query will return 4 records. Which one to download? We can use the plot method with sentinelsat to see where the image boundaries lie.
areas = api.to_geodataframe(products)
areas.plot(column='uuid', cmap=None)
We could also label the polygons using this code (acknowledgement https://stackoverflow.com/questions/38899190/geopandas-label-polygons)
ax = areas.plot(column='uuid', cmap=None, figsize=(20, 20))
areas.apply(lambda x: ax.annotate(s=x.uuid, xy=x.geometry.centroid.coords[0], ha='center'),axis=1)
This is only partially useful though, I have shaded the products by their id (uuid), but we need the boundary of the Isle of Wight plotted to give us the context.
If we use Geopandas we can plot two geodataframes on top of each other (like a GIS)
gdf2 = gpd.read_file(boundsdata)
f, ax = plt.subplots(1)
areas.plot(ax=ax,column='uuid',cmap=None,)
gdf2.plot(ax=ax)
plt.show()
I’ve written about this code before so for more explanation take a look here.
Based on the graphical results of the code from above, its clear that the product id we want is ’30b9e88b-94a7-45ec-9ceb-4fa8bf86be04′. You might wish to further explore intersection methods in geopandas to do this without plotting the data, but for now this method works nicely in a notebook.
Now that we have the product id we need we can just run this line to download the product to the location the notebook is (be careful on file size, its about 1gb)
api.download('30b9e88b-94a7-45ec-9ceb-4fa8bf86be04')
After download you can adapt the code found here to unzip and extract the data and you are ready to use the imagery. Below works on my example – you will have to change paths and potentially the Out[12] to match your notebook.
path_tozip = Out[12]['path'] ##notebook way - will not work in standalone script
import zipfile # https://stackoverflow.com/questions/3451111/unzipping-files-in-python
zip_ref = zipfile.ZipFile(path_tozip, 'r')
zip_ref.extractall(r'D:\sentinelsat')
zip_ref.close()
Recap
In this blog we have looked at using the sentinelsat Python API to download a Sentinel 2 level 2A image containing the Isle of Wight. We have used Folium and sentinelsat and geopandas to explore and download data.
All the code is available here.
Did I miss something? Let me know and I will add it to the list.
I am a freelancer able to help you with your projects. I offer consultancy, training and writing. I’d be delighted to hear from you. Please check out the books I have written on QGIS 3.4
https://www.packtpub.com/application-development/learn-qgis-fourth-edition
https://www.packtpub.com/application-development/qgis-quick-start-guide