Přeskočit navigaci

Harmonized Landsat Sentinel-2

SatelliteImagery EarthObservation AIforEarth NASA ESA

Satelitní snímky Severní Ameriky z družic Landsat-8 a Sentinel-2.

Produkt Harmonizovaná data z družic Landsat a Sentinel-2 (HLS) zahrnuje data z družic Landsat-8 a Sentinel-2 přizpůsobená jednotnému systému dlaždic s rozlišením 30 m., a to od roku 2013 až do současnosti v případě družice Landsat a od roku 2015 až do současnosti v případě družice Sentinel-2. Produkt HLS spravuje Národní úřad pro letectví a kosmonautiku (NASA).

Tuto datovou sadu udržuje společnost Ag-Analytics®. Společnost Ag-Analytics® nabízí také rozhraní API, které přijímá mnohoúhelník oblasti zájmu, rozsah dat a další možnosti a vrací zpracované obrázky pro jednotlivá spektrální pásma MSI, normalizovaný diferenční vegetační index (NDVI) a další metriky a také mozaiku s odfiltrovanými mraky.

Tato datová sada se aktualizuje každý měsíc.

Prostředky úložiště

Data se uchovávají v objektech blob v datacentru v oblasti Východní USA 2, a to v následujícím kontejneru objektů blob:

https://hlssa.blob.core.windows.net/hls

V rámci tohoto kontejneru se data uspořádávají následovně:

<folder>/HLS.<product>.T<tileid>.<daynum>.<version>_<subdataset>.tif

  • folder je L309 v případě družice Landsat, S309 v případě družice Sentinel-2.
  • product je L30 v případě družice Landsat, S30 v případě družice Sentinel-2.
  • tileid je čtyřmístný kód dlaždice ze systému dlaždic Sentinel-2.
  • daynum se skládá ze čtyřmístného roku a třímístného čísla dne v roce (od 001 do 365), např. 2019001 představuje 1. leden 2019.
  • version je vždycky v1.4.
  • subdataset je dvoumístný řetězec s indexováním od 1 označující podřízenou datovou sadu (viz níže).

Mapování zeměpisné šířky a délky na ID dlaždic najdete tady. Poznámkový blok v části věnované “přístupu k datům” ukazuje použití této tabulky k vyhledání ID dlaždice podle zeměpisné šířky a délky. ID dlaždic je možné zjistit také pomocí rozhraní Ag-Analytics® API.

K dispozici jsou data pro následující primární dlaždice:

[‘10 U’,‘11 U’,‘12 U’,‘13 U’,‘14 U’,‘15 U’,‘16 U’,‘10 T’,‘11 T’,‘12 T’,‘13 T’,‘14 T’,‘15 T’,‘16 T’,‘17 T’,‘18 T’,‘19 T’,‘10 S’,‘11 S’,‘12 S’,‘13 S’,‘14 S’,‘15 S’,‘16 S’,‘17 S’,‘18 S’,‘12 R’,‘13 R’,‘14 R’,‘15 R’,‘16 R’,‘17 R’]

Spektrální pásma jsou následující:

Název spektrálního pásma Číslo spektrálního pásma OLI Číslo spektrálního pásma MSI Číslo podřízené datové sady L30 Číslo podřízené datové sady S30
Pobřežní aerosol 1 1 01 01
Blue 2 2 02 02
Green 3 3 03 03
Red 4 4 04 04
Okraj červeného pásma 1 5 05
Okraj červeného pásma 2 6 06
Okraj červeného pásma 3 7 07
Blízké infračervené záření (široké spektrum) 8 08
Blízké infračervené záření (úzké spektrum) 5 8A 05 09
Krátkovlnné infračervené záření 1 6 11 06 10
Krátkovlnné infračervené záření 2 7 12 07 11
Vodní pára 9 12
Cirrus 9 10 08 13
Tepelné infračervené záření 1 10 09
Tepelné infračervené záření 2 11 10
QA 11 14

Například soubor s názvem HLS.S30.T16TDL.2019206.v1.4_01.tif by se nacházel v umístění https://hlssa.blob.core.windows.net/hls/S309/HLS.S30.T16TDL.2019206.v1.4_03.tif a představoval by data HLS Sentinel-2 (S30) pro dlaždici 16TDL (primární dlaždici 16T a podřízenou dlaždici DL) pro datovou sadu spektrálního pásma 03 (spektrální pásmo MSI 3, zelená) pro 206. den roku 2019.

Nabízíme také token SAS (sdílený přístupový podpis) jen pro čtení, který umožňuje přístup k datům HLS například přes nástroj BlobFuse, který umožňuje připojit kontejnery objektů blob jako jednotky:

st=2019-08-07T14%3A54%3A43Z&se=2050-08-08T14%3A54%3A00Z&sp=rl&sv=2018-03-28&sr=c&sig=EYNJCexDl5yxb1TxNH%2FzILznc3TiAnJq%2FPvCumkuV5U%3D

Pokyny k připojení pro Linux najdete tady.

Data HLS můžou spotřebovávat stovky terabajtů, takže rozsáhlé zpracování je nejlepší provádět v datacentru Azure v oblasti Východní USA 2, kde jsou obrázky uložené. Pokud data HLS využíváte v oblasti věd o životním prostředí, zvažte možnost zažádat o grant AI for Earth, který vám pomůže s vašimi požadavky na výpočetní prostředky.

Kontakt

Pokud máte k této datové sadě nějaké dotazy, obraťte se na aiforearthdatasets@microsoft.com.

Sdělení

MICROSOFT POSKYTUJE SLUŽBU AZURE OPEN DATASETS TAK, JAK JE. MICROSOFT V SOUVISLOSTI S VAŠÍM POUŽÍVÁNÍM DATOVÝCH SAD NEPOSKYTUJE ŽÁDNÉ ZÁRUKY, AŤ UŽ VÝSLOVNÉ NEBO PŘEDPOKLÁDANÉ, ANI JEJ NIJAK NEPODMIŇUJE. V ROZSAHU POVOLENÉM MÍSTNÍM ZÁKONEM MICROSOFT ODMÍTÁ JAKOUKOLI ODPOVĚDNOST ZA ŠKODY A ZTRÁTY ZPŮSOBENÉ VAŠÍM POUŽÍVÁNÍM DATOVÝCH SAD, VČETNĚ PŘÍMÝCH, NÁSLEDNÝCH, ZVLÁŠTNÍCH, NEPŘÍMÝCH, NÁHODNÝCH NEBO TRESTNÍCH ŠKOD.

Na tuto datovou sadu se vztahují původní podmínky, které Microsoft přijal se zdrojovými daty. Datová sada může obsahovat data pocházející z Microsoftu.

Access

Available inWhen to use
Azure Notebooks

Quickly explore the dataset with Jupyter notebooks hosted on Azure or your local machine.

Select your preferred service:

Azure Notebooks

Azure Notebooks

Package: Language: Python

Demo notebook for accessing HLS data on Azure

This notebook provides an example of accessing HLS (Harmonized Landsat Sentinel-2) data from blob storage on Azure, extracting image metadata using GDAL, and displaying an image using GDAL and rasterio.

HLS data are stored in the East US 2 data center, so this notebook will run more efficiently on the Azure compute located in East US 2. You don't want to download hundreds of terabytes to your laptop! If you are using HLS data for environmental science applications, consider applying for an AI for Earth grant to support your compute requirements.

HLS data on Azure are managed by Ag-Analytics. Ag-Analytics also provides an API which allows the caller to query to perform spatial queries over the HLS archive, as well as querying for additional data such as cloud cover and Normalized Difference Vegetation Index (NDVI). Ag-Analytics also provides an API to retrieve tile IDs matching spatial queries.

Imports and environment

In [2]:
# Standard-ish packages
import requests
import re
import numpy as np
import urllib
import io
import matplotlib.pyplot as plt
import pandas as pd

# Less standard, but all of the following are pip- or conda-installable
import rasterio

# pip install azure-storage-blob
from azure.storage.blob import ContainerClient

from osgeo import gdal,osr

# Storage locations are documented at http://aka.ms/ai4edata-hls
hls_container_name = 'hls'
hls_account_name = 'hlssa'
hls_account_url = 'https://' + hls_account_name + '.blob.core.windows.net/'
hls_blob_root = hls_account_url + hls_container_name

# This file is provided by NASA; it indicates the lat/lon extents of each
# hls tile.
#
# The file originally comes from:
#
# https://hls.gsfc.nasa.gov/wp-content/uploads/2016/10/S2_TilingSystem2-1.txt
#
# ...but as of 8/2019, there is a bug with the column names in the original file, so we
# access a copy with corrected column names.
hls_tile_extents_url = 'https://ai4edatasetspublicassets.blob.core.windows.net/assets/S2_TilingSystem2-1.txt?st=2019-08-23T03%3A25%3A57Z&se=2028-08-24T03%3A25%3A00Z&sp=rl&sv=2018-03-28&sr=b&sig=KHNZHIJuVG2KqwpnlsJ8truIT5saih8KrVj3f45ABKY%3D'

# Load this file into a table, where each row is:
#
# Tile ID, Xstart, Ystart, UZ, EPSG, MinLon, MaxLon, MinLon, MaxLon
print('Reading tile extents...')
s = requests.get(hls_tile_extents_url).content
hls_tile_extents = pd.read_csv(io.StringIO(s.decode('utf-8')),delimiter=r'\s+')
print('Read tile extents for {} tiles'.format(len(hls_tile_extents)))

# Read-only shared access signature (SAS) URL for the hls container
hls_sas_token = 'st=2019-08-07T14%3A54%3A43Z&se=2050-08-08T14%3A54%3A00Z&sp=rl&sv=2018-03-28&sr=c&sig=EYNJCexDl5yxb1TxNH%2FzILznc3TiAnJq%2FPvCumkuV5U%3D'

hls_container_client = ContainerClient(account_url=hls_account_url, 
                                         container_name=hls_container_name,
                                         credential=None)
                                

%matplotlib inline
Reading tile extents...
Read tile extents for 56686 tiles

Functions

In [3]:
def get_hls_tile(blob_url):
    """
    Given a URL pointing to an HLS image in blob storage, load that image via GDAL
    and return both data and metadata.
    """    
    
    formatted_gdal_bloburl='/{}/{}'.format('vsicurl',blob_url)
    
    tile_open = gdal.Open(formatted_gdal_bloburl)
    data = tile_open.GetRasterBand(1)
    ndv,xsize,ysize = data.GetNoDataValue(),tile_open.RasterXSize,tile_open.RasterYSize
    
    projection = osr.SpatialReference()
    projection.ImportFromWkt(tile_open.GetProjectionRef())
    
    datatype = data.DataType
    datatype = gdal.GetDataTypeName(datatype)  
    data_array = data.ReadAsArray()

    return ndv,xsize,ysize,projection,data_array


def list_available_tiles(prefix):
    """
    List all blobs in an Azure blob container matching a prefix.  
    
    We'll use this to query tiles by location and year.
    """
    
    files = []
    generator = hls_container_client.list_blobs(name_starts_with=prefix)
    for blob in generator:
        files.append(blob.name)
    return files

    
def lat_lon_to_hls_tile_id(lat,lon):
    """
    Get the hls tile ID for a given lat/lon coordinate pair
    """  
    
    found_matching_tile = False

    for i_row,row in hls_tile_extents.iterrows():
        found_matching_tile = lat >= row.MinLat and lat <= row.MaxLat \
        and lon >= row.MinLon and lon <= row.MaxLon
        if found_matching_tile:
            break
    
    if not found_matching_tile:
        return None
    else:
        return row.TilID

Find a tile for a given location and date

In [4]:
# Specify a location and year of interest
lat = 47.6101; lon = -122.2015 # Bellevue, WA

year = '2019'
daynum = '109'    # 1-indexed day-of-year
folder = 'S309'   # 'S309' for Sentinel, 'L309' for Landsat
product = 'S30'   # 'S30' for Sentinel, 'L30' for Landsat
year = '2019'

tile_id = lat_lon_to_hls_tile_id(lat,lon)
assert tile_id is not None, 'Invalid lat/lon'
prefix = folder + '/HLS.' + product + '.T' + tile_id + '.' + year

print('Finding tiles with prefix {}'.format(prefix))
matches = list_available_tiles(prefix)
assert len(matches) > 0, 'No matching tiles'

blob_name = matches[0]
print('Found {} matching tiles, using file {}'.format(len(matches),blob_name))
Finding tiles with prefix S309/HLS.S30.T10TET.2019
Found 1918 matching tiles, using file S309/HLS.S30.T10TET.2019001.v1.4_01.tif

...or build a tile path from components

In [5]:
lat = 47.6101; lon = -122.2015 # Bellevue, WA

year    = '2019'
daynum  = '001'   # 1-indexed day-of-year
folder  = 'S309'  # 'S309' for Sentinel, 'L309' for Landsat
product = 'S30'   # 'S30' for Sentinel, 'L30' for Landsat
band    = '01'
tile_id = '10TET' # See hls.gsfc.nasa.gov/wp-content/uploads/2016/10/S2_TilingSystem2-1.txt
version = 'v1.4'  # Currently always v1.4

blob_name = folder + '/HLS.' + product + '.T' + tile_id + '.' + year + daynum + '.' + version \
    + '_' + band + '.tif'

print('Using file {}'.format(blob_name))
Using file S309/HLS.S30.T10TET.2019001.v1.4_01.tif

Access one band of the selected image using GDAL's virtual file system (vsicurl)

In [6]:
gdal.SetConfigOption('GDAL_HTTP_UNSAFESSL', 'YES')
blob_url = hls_blob_root + '/' + blob_name
print('Reading tile from {}'.format(blob_url))
ndv,xsize,ysize,projection,data_array = get_hls_tile(blob_url)

print('No-data value: {}'.format(ndv))
print('\nSize: {},{}'.format(xsize,ysize))
print('\nProjection:\n{}'.format(projection))
Reading tile from https://hlssa.blob.core.windows.net/hls/S309/HLS.S30.T10TET.2019001.v1.4_01.tif
No-data value: -1000.0

Size: 3660,3660

Projection:
PROJCS["UTM Zone 10, Northern Hemisphere",
    GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",
        DATUM["Not_specified_based_on_WGS_84_spheroid",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-123],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

Display a Sentinel-2 image using rasterio and vsicurl

In [7]:
# Bands 2, 3, and 4 are B, G, and R in Sentinel-2 HLS images

base_url = '/vsicurl/' + hls_blob_root + '/' + blob_name
band2_url = re.sub('_(\d+).tif','_02.tif',base_url)
band3_url = re.sub('_(\d+).tif','_03.tif',base_url)
band4_url = re.sub('_(\d+).tif','_04.tif',base_url)
print('Reading bands from:\n{}\n{}\n{}'.format(band2_url,band3_url,band4_url))

band2 = rasterio.open(band2_url)
band3 = rasterio.open(band3_url)
band4 = rasterio.open(band4_url)

norm_value = 2000
image_data = []
for band in [band4,band3,band2]:
    band_array = band.read(1)
    band_array = band_array / norm_value
    image_data.append(band_array)
    band.close()

rgb = np.dstack((image_data[0],image_data[1],image_data[2]))
np.clip(rgb,0,1,rgb)
plt.imshow(rgb)
Reading bands from:
/vsicurl/https://hlssa.blob.core.windows.net/hls/S309/HLS.S30.T10TET.2019001.v1.4_02.tif
/vsicurl/https://hlssa.blob.core.windows.net/hls/S309/HLS.S30.T10TET.2019001.v1.4_03.tif
/vsicurl/https://hlssa.blob.core.windows.net/hls/S309/HLS.S30.T10TET.2019001.v1.4_04.tif
Out[7]:
<matplotlib.image.AxesImage at 0x2249d8e6f48>