Navigatie overslaan

Harmonized Landsat Sentinel-2

SatelliteImagery EarthObservation AIforEarth NASA ESA

Satellietbeelden van de satellieten Landsat-8 en Sentinel-2 voor Noord-Amerika.

Het HLS-product (Harmonized Landsat Sentinel-2) bevat gegevens die afkomstig zijn van de satellieten Landsat-8 en Sentinel-2, die zijn uitgelijnd op een algemeen tegelsysteem bij een resolutie van 30 m, van 2013 tot heden voor Landsat en van 2015 tot heden voor Sentinel-2. HLS wordt beheerd door de National Aeronautics and Space Administration (NASA).

Deze gegevensset wordt beheerd door Ag-Analytics®. Ag-Analytics® biedt ook een API waarmee een AOI-veelhoek (Area of Interest), gegevensbereik en andere opties worden geaccepteerd en verwerkte afbeeldingen worden geretourneerd voor afzonderlijke MSI-banden, alsmede een Normalized Difference Vegetation Index en andere metrische gegevens en mozaïeken waarin de wolken zijn gefilterd.

Deze gegevensset wordt wekelijks bijgewerkt.

Opslagresources

Gegevens worden opgeslagen in blobs in het datacenter in de regio US - oost 2, in de volgende blobcontainer:

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

In die container worden gegevens als volgt geordend:

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

  • folder is L309 voor Landsat, S309 voor Sentinel-2
  • product is L30 voor Landsat, S30 voor Sentinel-2
  • tileid is een viercijferige tegelcode die afkomstig is uit het tegelsysteem van Sentinel-2
  • daynum is een viercijferig jaar met een driecijferige dag van het jaar (van 001 tot en met 365); 2019001 staat bijvoorbeeld voor 1 januari 2019
  • version is altijd v1.4
  • subdataset is een op 1 geïndexeerde tekenreeks die uit twee tekens bestaat en die aangeeft dat er een subgegevensset wordt gebruikt (zie onder)

Kijk hier voor toewijzing van breedte-/lengtegraden aan tegel-id’s; in het notitieblok onder “Gegevenstoegang” wordt het gebruik van deze tabel voor het opzoeken van een tegel-id op basis van breedte-/lengtegraden getoond. Tegel-id’s vindt u ook met behulp van de Ag-Analytics® API.

Gegevens worden opgegeven voor de volgende primaire tegels:

[‘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’]

De volgende banden zijn beschikbaar:

Bandnaam OLI-bandnummer MSI-bandnummer L30-subgegevenssetnummer S30-subgegevenssetnummer
Aerosolgehalte aan de kust 1 1 01 01
Blue 2 2 02 02
Green 3 3 03 03
Red 4 4 04 04
Rood-rand 1 5 05
Rood-rand 2 6 06
Rood-rand 3 7 07
NIR breed 8 08
NIR smal 5 8A 05 09
SWIR 1 6 11 06 10
SWIR 2 7 12 07 11
Waterdamp 9 12
Cirrus 9 10 08 13
Thermisch infrarood 1 10 09
Thermisch infrarood 2 11 10
QA 11 14

De bestandsnaam HLS.S30.T16TDL.2019206.v1.4_01.tif is bijvoorbeeld te vinden via https://hlssa.blob.core.windows.net/hls/S309/HLS.S30.T16TDL.2019206.v1.4_03.tif en dit bestand vertegenwoordigt Sentinel-2 (S30) HLS-gegevens voor tegel 16TDL (primaire tegel 16T, subtegel DL) voor gegevenssetband 03 (MSI-band 3, groen) voor de 206e dag van 2019.

Ook beiden een alleen-lezen SAS-token (Shared Access Signature) om toegang tot HLS-gegevens toe te staan via bijvoorbeeld BlobFuse, waardoor u blobcontainers als schijven kunt koppelen:

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

Hier vindt u koppelingsinstructies voor Linux.

Voor HLS-gegevens kunnen honderden terabytes worden verbruikt, dus u kunt grootschalige verwerking het beste uitvoeren in het Azure-datacenter in de regio US - oost 2, waarin de afbeeldingen worden opgeslagen. Als u HLS-gegevens gebruikt voor milieuwetenschappelijke toepassingen, is het raadzaam om u aan te melden voor een AI for Earth-toelage ter ondersteuning van uw computevereisten.

Contactpersoon

Voor vragen over deze gegevensset kunt contact opnemen met aiforearthdatasets@microsoft.com.

Mededelingen

AZURE OPEN GEGEVENSSETS WORDEN DOOR MICROSOFT ONGEWIJZIGD GELEVERD. MICROSOFT GEEFT GEEN GARANTIES, EXPLICIET OF IMPLICIET, ZEKERHEDEN OF VOORWAARDEN MET BETREKKING TOT HET GEBRUIK VAN DE GEGEVENSSETS. VOOR ZOVER IS TOEGESTAAN ONDER HET TOEPASSELIJKE RECHT, WIJST MICROSOFT ALLE AANSPRAKELIJKHEID AF VOOR SCHADE OF VERLIEZEN, WAARONDER GEVOLGSCHADE OF DIRECTE, SPECIALE, INDIRECTE, INCIDENTELE OF PUNITIEVE SCHADE DIE VOORTVLOEIT UIT HET GEBRUIK VAN DE GEGEVENSSETS.

Deze gegevensset wordt geleverd onder de oorspronkelijke voorwaarden dat Microsoft de brongegevens heeft ontvangen. De gegevensset kan gegevens bevatten die afkomstig zijn van Microsoft.

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 [1]:
# 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'

# 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)))

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 [2]:
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 [ ]:
# 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))

...or build a tile path from components

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

year    = '2019'
daynum  = '001'   # 1-indexed day-of-year
folder  = 'S30'   # 'S30' for Sentinel, 'L30' 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 S30/HLS.S30.T10TET.2019001.v1.4_01.tif

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

In [4]:
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/S30/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,
            AUTHORITY["EPSG","9122"]]],
    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"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]

Display a Sentinel-2 image using rasterio and vsicurl

In [5]:
# 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/S30/HLS.S30.T10TET.2019001.v1.4_02.tif
/vsicurl/https://hlssa.blob.core.windows.net/hls/S30/HLS.S30.T10TET.2019001.v1.4_03.tif
/vsicurl/https://hlssa.blob.core.windows.net/hls/S30/HLS.S30.T10TET.2019001.v1.4_04.tif
Out[5]:
<matplotlib.image.AxesImage at 0x159f5af1a60>

Read the jpeg thumbnail from a COG file without reading the whole file

In [13]:
rgb_urls = [band4_url, band3_url, band2_url]

thumbnail_data = []

# url = rgb_urls[0]
for url in rgb_urls:
    
    # From:
    #
    # https://automating-gis-processes.github.io/CSC/notebooks/L5/read-cogs.html
    with rasterio.open(url) as raster:
        
        # List of overviews from largest to smallest
        oviews = raster.overviews(1)
    
        # Retrieve the second-largest thumbnail
        decimation_level = oviews[1]
        h = int(raster.height/decimation_level)
        w = int(raster.width/decimation_level)
        
        thumbnail_channel = raster.read(1, out_shape=(1, h, w)) / norm_value
        thumbnail_data.append(thumbnail_channel)

rgb = np.dstack((thumbnail_data[0],thumbnail_data[1],thumbnail_data[2]))
np.clip(rgb,0,1,rgb)
plt.imshow(rgb)
Out[13]:
<matplotlib.image.AxesImage at 0x15993868a30>