78 read raster

Interactive Extraction and Visualization of AWS Open Geospatial Data

Uncomment the following line to install leafmap if needed.

In [1]:

# %pip install -U leafmap

# %pip install -U leafmap

In [2]:

from leafmap import leafmap

from leafmap import leafmap

Set custom STAC endpoints.

In [3]:

catalogs = {
    "Element84 Earth Search": "https://earth-search.aws.element84.com/v1",
    "Microsoft Planetary Computer": "https://planetarycomputer.microsoft.com/api/stac/v1",
}

catalogs = { "Element84 Earth Search": "https://earth-search.aws.element84.com/v1", "Microsoft Planetary Computer": "https://planetarycomputer.microsoft.com/api/stac/v1", }

Display the STAC search GUI. Pan and zoom to the area of interest and use the drawing tools to draw a bounding box or polygon. Select a STAC catalog and click the Collections button to retrieve the collections, then click on the Items button to retrieve the items within the bounding box or polygon. Select an item from the dropdown list and click the Display button to display the item on the map.

In [4]:

m = leafmap.Map(center=[37.7452, -122.4108], zoom=12, catalog_source=catalogs)
m.add("stac")
m

m = leafmap.Map(center=[37.7452, -122.4108], zoom=12, catalog_source=catalogs) m.add("stac") m

Display the information of the selected item.

In [5]:

# m.stac_item

# m.stac_item

Alternatively, you can search the STAC catalog programmatically by providing a bounding box, time range, and other filters. The example below use the Earth Search STAC endpoint by Element 84 for searching for AWS Open Data.

In [6]:

url = "https://earth-search.aws.element84.com/v1/"
collection = "sentinel-2-l2a"
time_range = "2023-04-01/2023-07-31"
bbox = [-122.491, 37.7208, -122.411, 37.7786]

url = "https://earth-search.aws.element84.com/v1/" collection = "sentinel-2-l2a" time_range = "2023-04-01/2023-07-31" bbox = [-122.491, 37.7208, -122.411, 37.7786]

Search the STAC catalog and return the results as an ItemCollection.

In [7]:

search = leafmap.stac_search(
    url=url,
    max_items=5,
    collections=[collection],
    bbox=bbox,
    datetime=time_range,
    query={"eo:cloud_cover": {"lt": 10}},
    get_collection=True,
)
# search

search = leafmap.stac_search( url=url, max_items=5, collections=[collection], bbox=bbox, datetime=time_range, query={"eo:cloud_cover": {"lt": 10}}, get_collection=True, ) # search

Search the STAC catalog and return the results as a dictionary of assets.

In [8]:

search = leafmap.stac_search(
    url=url,
    max_items=5,
    collections=[collection],
    bbox=bbox,
    datetime=time_range,
    query={"eo:cloud_cover": {"lt": 10}},
    get_assets=True,
)
# search

search = leafmap.stac_search( url=url, max_items=5, collections=[collection], bbox=bbox, datetime=time_range, query={"eo:cloud_cover": {"lt": 10}}, get_assets=True, ) # search

Get the first item in the collection.

In [9]:

name, item = next(iter(search.items()))
name

name, item = next(iter(search.items())) name

Out[9]:

'S2A_10SEG_20230618_0_L2A'

Retrieve the item's assets, which are links to the actual data files.

In [10]:

item

item

Out[10]:

{'aot': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/AOT.tif',
 'blue': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B02.tif',
 'coastal': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B01.tif',
 'granule_metadata': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/granule_metadata.xml',
 'green': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B03.tif',
 'nir': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B08.tif',
 'nir08': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B8A.tif',
 'nir09': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B09.tif',
 'red': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B04.tif',
 'rededge1': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B05.tif',
 'rededge2': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B06.tif',
 'rededge3': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B07.tif',
 'scl': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/SCL.tif',
 'swir16': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B11.tif',
 'swir22': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B12.tif',
 'thumbnail': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/thumbnail.jpg',
 'tileinfo_metadata': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/tileinfo_metadata.json',
 'visual': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/TCI.tif',
 'wvp': 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/WVP.tif',
 'aot-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/AOT.jp2',
 'blue-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/B02.jp2',
 'coastal-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/B01.jp2',
 'green-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/B03.jp2',
 'nir-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/B08.jp2',
 'nir08-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/B8A.jp2',
 'nir09-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/B09.jp2',
 'red-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/B04.jp2',
 'rededge1-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/B05.jp2',
 'rededge2-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/B06.jp2',
 'rededge3-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/B07.jp2',
 'scl-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/SCL.jp2',
 'swir16-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/B11.jp2',
 'swir22-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/B12.jp2',
 'visual-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/TCI.jp2',
 'wvp-jp2': 's3://sentinel-s2-l2a/tiles/10/S/EG/2023/6/18/0/WVP.jp2'}

Retrieve the STAT item's URLs.

In [11]:

search = leafmap.stac_search(
    url=url,
    max_items=5,
    collections=[collection],
    bbox=bbox,
    datetime=time_range,
    query={"eo:cloud_cover": {"lt": 10}},
    get_links=True,
)
search

search = leafmap.stac_search( url=url, max_items=5, collections=[collection], bbox=bbox, datetime=time_range, query={"eo:cloud_cover": {"lt": 10}}, get_links=True, ) search

Out[11]:

['https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/S2A_10SEG_20230618_0_L2A',
 'https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/S2B_10SEG_20230424_0_L2A',
 'https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a/items/S2B_10SEG_20230414_0_L2A']

Check the band names of the selected item.

In [12]:

url = search[0]
bands = leafmap.stac_bands(url)
bands[:10]

url = search[0] bands = leafmap.stac_bands(url) bands[:10]

Out[12]:

['aot',
 'blue',
 'coastal',
 'green',
 'nir',
 'nir08',
 'nir09',
 'red',
 'rededge1',
 'rededge2']

Display the selected item on the map.

In [13]:

m = leafmap.Map()
m.add_stac_layer(url, bands=["nir", "red", "green"], name="Sentinel-2")
m

m = leafmap.Map() m.add_stac_layer(url, bands=["nir", "red", "green"], name="Sentinel-2") m

Use the drawing tools to draw a small bounding box on the image.

In [14]:

if m.user_roi is not None:
    roi = m.user_roi_bounds()
else:
    roi = [-122.5315, 37.6882, -122.3523, 37.8166]

if m.user_roi is not None: roi = m.user_roi_bounds() else: roi = [-122.5315, 37.6882, -122.3523, 37.8166]

Specify the bands to use.

In [15]:

bands = ["nir", "red", "green"]

bands = ["nir", "red", "green"]

Display the COG URL.

In [16]:

item["nir"]

item["nir"]

Out[16]:

'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/10/S/EG/2023/6/S2A_10SEG_20230618_0_L2A/B08.tif'

Extract one single band within the bounding box as an a numpy array.

In [17]:

array = leafmap.read_raster(item["nir"], window=roi, coord_crs="epsg:4326")

array = leafmap.read_raster(item["nir"], window=roi, coord_crs="epsg:4326")

Check the shape of the array.

In [18]:

array.shape

array.shape

Out[18]:

(1, 1434, 1570)

Extract multiple bands within the bounding box as an a numpy array.

In [19]:

sources = [item["nir"], item["red"], item["green"]]
array = leafmap.read_rasters(sources, window=roi, coord_crs="epsg:4326")

sources = [item["nir"], item["red"], item["green"]] array = leafmap.read_rasters(sources, window=roi, coord_crs="epsg:4326")

Check the shape of the array.

In [20]:

array.shape

array.shape

Out[20]:

(3, 1434, 1570)

Convert the numpy array to a Cloud Optimized GeoTIFF (COG).

In [21]:

leafmap.numpy_to_cog(
    array, "s2.tif", bounds=roi, profile=item["nir"], coord_crs="epsg:4326"
)

leafmap.numpy_to_cog( array, "s2.tif", bounds=roi, profile=item["nir"], coord_crs="epsg:4326" )

Display the image on the map.

In [22]:

m = leafmap.Map()
m.add_raster("s2.tif", band=[1, 2, 3], vmin=0, vmax=4000, layer_name="Subset")
m

m = leafmap.Map() m.add_raster("s2.tif", band=[1, 2, 3], vmin=0, vmax=4000, layer_name="Subset") m