Bài 29: Truy cập dữ liệu Copernicus Data Space Ecosystem (CDSE)

conda env config vars set GDAL_HTTP_TCP_KEEPALIVE=YES \
  AWS_S3_ENDPOINT=eodata.dataspace.copernicus.eu \
  AWS_ACCESS_KEY_ID=your_access_key \
  AWS_SECRET_ACCESS_KEY=your_secret_key \
  AWS_HTTPS=YES \
  AWS_VIRTUAL_HOSTING=FALSE \
  GDAL_HTTP_UNSAFESSL=YES

# Thay thế bằng keys thực tế của bạn
import os 
access_key = 'your_actual_access_key'  # Thay bằng Access Key từ CDSE
secret_key = 'your_actual_secret_key'  # Thay bằng Secret Key từ CDSE

# Cấu hình environment variables cho GDAL và AWS S3
os.environ["GDAL_HTTP_TCP_KEEPALIVE"] = "YES"
os.environ["AWS_S3_ENDPOINT"] = "eodata.dataspace.copernicus.eu"
os.environ["AWS_HTTPS"] = "YES"
os.environ["AWS_VIRTUAL_HOSTING"] = "FALSE"
os.environ["GDAL_HTTP_UNSAFESSL"] = "YES"

# Uncomment và sử dụng keys của bạn
# os.environ["AWS_ACCESS_KEY_ID"] = access_key
# os.environ["AWS_SECRET_ACCESS_KEY"] = secret_key

from pystac_client import Client

# Connect to the Copernicus Data Space STAC catalog
catalog = Client.open("https://catalogue.dataspace.copernicus.eu/stac")
# you can see all collections in the catalog
collections = catalog.get_all_collections()
collections_ids = [col.id for col in collections]
print(collections_ids[10:15])  # print some collection ids

['cop-dem-glo-90-dged-cog', 'sentinel-3-olci-2-wrr-nrt', 'sentinel-2-l2a', 'sentinel-1-slc', 'sentinel-2-gri-l1c-gcp']

# define a bounding box or area of interest
bbox = [11.439263980173, 47.81384831137465, 11.475186504136818, 47.83147903246192] # this location is in Germany
# search for items in sentinel-2-l2a collection
items = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=bbox,
    datetime="2023-05-01/2023-09-30",
    query={"eo:cloud_cover": {"lt": 20}},
).item_collection()
print(f"Found {len(items)} items")

Found 31 items

# get the first image id as example 
# print the properties of the first item
list(items)[:1][0].properties # you can use this info to filter queries

{'gsd': 10,
 'created': '2025-03-21T17:41:10.043000Z',
 'updated': '2025-04-21T00:16:33.282435Z',
 'datetime': '2023-09-28T10:17:19.024Z',
 'platform': 'sentinel-2b',
 'grid:code': 'MGRS-32UPU',
 'published': '2025-04-21T00:16:33.282435Z',
 'statistics': {'water': 1.483443,
  'nodata': 1.7e-05,
  'dark_area': 0.003301,
  'vegetation': 80.213308,
  'thin_cirrus': 0.805717,
  'cloud_shadow': 0.0,
  'unclassified': 0.000269,
  'not_vegetated': 17.480843,
  'high_proba_clouds': 0.002183,
  'medium_proba_clouds': 0.010677,
  'saturated_defective': 0.0},
 'instruments': ['msi'],
 'auth:schemes': {'s3': {'type': 's3'},
  'oidc': {'type': 'openIdConnect',
   'openIdConnectUrl': 'https://identity.dataspace.copernicus.eu/auth/realms/CDSE/.well-known/openid-configuration'}},
 'end_datetime': '2023-09-28T10:17:19.024Z',
 'product:type': 'S2MSI2A',
 'view:azimuth': 269.6725837191976,
 'constellation': 'sentinel-2',
 'eo:snow_cover': 0.000255,
 'eo:cloud_cover': 0.818577,
 'start_datetime': '2023-09-28T10:17:19.024Z',
 'sat:orbit_state': 'descending',
 'storage:schemes': {'cdse-s3': {'type': 'custom-s3',
   'title': 'Copernicus Data Space Ecosystem S3',
   'platform': 'https://eodata.dataspace.copernicus.eu',
   'description': 'This endpoint provides access to EO data which is stored on the object storage of both CloudFerro Cloud and OpenTelekom Cloud (OTC). See the [documentation](https://documentation.dataspace.copernicus.eu/APIs/S3.html) for more information, including how to get credentials.',
   'requester_pays': False},
  'creodias-s3': {'type': 'custom-s3',
   'title': 'CREODIAS S3',
   'platform': 'https://eodata.cloudferro.com',
   'description': 'Comprehensive Earth Observation Data (EODATA) archive offered by CREODIAS as a commercial part of CDSE, designed to provide users with access to a vast repository of satellite data without predefined quota limits.',
   'requester_pays': True}},
 'eopf:datatake_id': 'GS2B_20230928T101719_034268_N05.10',
 'processing:level': 'L2',
 'view:sun_azimuth': 167.450831069454,
 'eopf:datastrip_id': 'S2B_OPER_MSI_L2A_DS_S2RP_20241026T142623_S20230928T102123_N05.10',
 'processing:version': '05.10',
 'product:timeliness': 'PT24H',
 'sat:absolute_orbit': 34268,
 'sat:relative_orbit': 65,
 'view:sun_elevation': 39.0424772062148,
 'processing:datetime': '2024-10-26T14:26:23.000000Z',
 'processing:facility': 'ESA',
 'eopf:instrument_mode': 'INS-NOBS',
 'eopf:origin_datetime': '2025-03-21T17:41:10.043000Z',
 'view:incidence_angle': 4.842546212382268,
 'product:timeliness_category': 'NRT',
 'sat:platform_international_designator': '2017-013A'}

from odc.stac import stac_load
# read the image data from the STAC items
data = stac_load(
    list(items),
    chunks={"x": 1024, "y": 1024},
    resolution=10,
    bbox=bbox, # you can also use selected band names here. here we load all bands
) # always return xarray.Dataset
data

img_id = list(items)[:1][0].id
search = catalog.search(collections=["sentinel-2-l2a"], ids=[img_id])
items = search.item_collection()
data = stac_load(
    list(items),
    chunks={"x": 1024, "y": 1024},
    resolution=10
) # always return xarray.Dataset
data

# Trích xuất band Blue (B02) với độ phân giải 10m
blue = data['B02_10m'].squeeze()  # Loại bỏ dimensions đơn lẻ
blue = blue.compute()  # Tải dữ liệu vào memory

print(f"🔵 Band Blue (B02):")
print(f"   Kích thước: {blue.shape}")
print(f"   Kiểu dữ liệu: {blue.dtype}")
print(f"   Giá trị min: {blue.min().values}")
print(f"   Giá trị max: {blue.max().values}")
print(f"   Giá trị trung bình: {blue.mean().values:.2f}")

# Tạo figure với subplots
fig, axes = plt.subplots(1, 2, figsize=(15, 6))

# Plot 1: Band Blue
im1 = axes[0].imshow(blue, cmap='Blues', vmin=0, vmax=4000)
axes[0].set_title('Sentinel-2 Band Blue (B02)', fontsize=12, fontweight='bold')
axes[0].set_xlabel('Pixel X')
axes[0].set_ylabel('Pixel Y')
plt.colorbar(im1, ax=axes[0], label='Reflectance')

# Plot 2: Histogram
blue_flat = blue.values.flatten()
blue_flat = blue_flat[~np.isnan(blue_flat)]  # Loại bỏ NaN values
axes[1].hist(blue_flat, bins=50, alpha=0.7, color='blue', edgecolor='black')
axes[1].set_title('Histogram - Band Blue Distribution', fontsize=12, fontweight='bold')
axes[1].set_xlabel('Reflectance Value')
axes[1].set_ylabel('Frequency')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("📊 Đã hoàn thành trực quan hóa band Blue!")

# Trích xuất các band cần thiết cho tính NDVI
red = data['B04_10m'].squeeze().compute()    # Band Red
nir = data['B08_10m'].squeeze().compute()    # Band NIR

# Tính NDVI (Normalized Difference Vegetation Index)
# NDVI = (NIR - Red) / (NIR + Red)
ndvi = (nir - red) / (nir + red)

print(f"🌱 NDVI Statistics:")
print(f"   Min: {ndvi.min().values:.3f}")
print(f"   Max: {ndvi.max().values:.3f}")
print(f"   Mean: {ndvi.mean().values:.3f}")
print(f"   Std: {ndvi.std().values:.3f}")

# Trực quan hóa NDVI
fig, ax = plt.subplots(figsize=(10, 8))
im = ax.imshow(ndvi, cmap='RdYlGn', vmin=-0.5, vmax=1.0)
ax.set_title('NDVI - Chỉ số thực vật chuẩn hóa', fontsize=14, fontweight='bold')
ax.set_xlabel('Pixel X')
ax.set_ylabel('Pixel Y')
cbar = plt.colorbar(im, ax=ax, label='NDVI Value')

# Thêm legend cho NDVI
ax.text(0.02, 0.98, 'NDVI Interpretation:\n• < 0: Water/Snow\n• 0-0.2: Bare soil\n• 0.2-0.5: Low vegetation\n• > 0.5: Dense vegetation', 
        transform=ax.transAxes, fontsize=10, verticalalignment='top',
        bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

plt.tight_layout()
plt.show()

print("🎯 Đã tính và trực quan hóa NDVI thành công!")

# Trích xuất các band RGB
green = data['B03_10m'].squeeze().compute()  # Band Green
# red và blue đã có từ trước

# Chuẩn hóa giá trị về 0-1 cho hiển thị RGB
def normalize_band(band, percentile_range=(2, 98)):
    """Chuẩn hóa band sử dụng percentile stretch"""
    p_low, p_high = np.percentile(band.values[~np.isnan(band.values)], percentile_range)
    band_norm = (band - p_low) / (p_high - p_low)
    return np.clip(band_norm, 0, 1)

# Chuẩn hóa các band
red_norm = normalize_band(red)
green_norm = normalize_band(green)
blue_norm = normalize_band(blue)

# Tạo RGB composite
rgb = np.stack([red_norm, green_norm, blue_norm], axis=-1)

# Hiển thị RGB composite
fig, axes = plt.subplots(1, 2, figsize=(16, 8))

# Plot RGB
axes[0].imshow(rgb)
axes[0].set_title('Sentinel-2 True Color (RGB)', fontsize=12, fontweight='bold')
axes[0].set_xlabel('Pixel X')
axes[0].set_ylabel('Pixel Y')

# Plot NDVI để so sánh
im2 = axes[1].imshow(ndvi, cmap='RdYlGn', vmin=-0.5, vmax=1.0)
axes[1].set_title('NDVI Comparison', fontsize=12, fontweight='bold')
axes[1].set_xlabel('Pixel X')
axes[1].set_ylabel('Pixel Y')
plt.colorbar(im2, ax=axes[1], label='NDVI Value')

plt.tight_layout()
plt.show()

print("🌈 Đã tạo và hiển thị RGB composite thành công!")

# Tìm kiếm ảnh trong collection sentinel-2-l2a
items = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=bbox,
    datetime="2023-05-01/2023-09-30",  # Thời gian từ tháng 5 đến tháng 9/2023
    query={"eo:cloud_cover": {"lt": 20}},  # Độ phủ mây < 20%
).item_collection()

print(f"🔍 Đã tìm kiếm xong với các tiêu chí:")
print(f"   - Collection: sentinel-2-l2a")
print(f"   - Thời gian: 2023-05-01 đến 2023-09-30")
print(f"   - Độ phủ mây: < 20%")

access_key = 'example_access_key' # Please replace with your actual access key
secret_key = 'example_secret_key' # Please replace with your actual secret key
os.environ["GDAL_HTTP_TCP_KEEPALIVE"] = "YES"
os.environ["AWS_S3_ENDPOINT"] = "eodata.dataspace.copernicus.eu"
os.environ["AWS_HTTPS"] = "YES"
os.environ["AWS_VIRTUAL_HOSTING"] = "FALSE"
os.environ["GDAL_HTTP_UNSAFESSL"] = "YES"
# os.environ["AWS_ACCESS_KEY_ID"] = access_key
# os.environ["AWS_SECRET_ACCESS_KEY"] = secret_key

# Connect to the Copernicus Data Space STAC catalog
catalog = Client.open("https://catalogue.dataspace.copernicus.eu/stac")
# you can see all collections in the catalog
collections = catalog.get_all_collections()
collections_ids = [col.id for col in collections]
print(collections_ids)  # print all collection ids

['sentinel-3-sl-2-aod-nrt', 'sentinel-1-global-mosaics', 'sentinel-3-olci-2-wfr-nrt', 'cop-dem-glo-30-dged-cog', 'sentinel-1-slc-wv', 'sentinel-2-gri-l1c', 'sentinel-2-l1c', 'sentinel-2-global-mosaics', 'sentinel-3-sl-2-wst-nrt', 'sentinel-3-olci-1-efr-nrt', 'cop-dem-glo-90-dged-cog', 'sentinel-3-olci-2-wrr-nrt', 'sentinel-2-l2a', 'sentinel-1-slc', 'sentinel-2-gri-l1c-gcp', 'sentinel-3-sl-2-lst-ntc', 'sentinel-3-olci-1-efr-ntc', 'sentinel-3-olci-2-wfr-ntc', 'sentinel-3-olci-2-lfr-nrt', 'sentinel-3-sl-2-wst-ntc', 'sentinel-3-olci-1-err-nrt', 'sentinel-3-olci-2-lfr-ntc', 'sentinel-3-sl-2-frp-ntc', 'sentinel-3-olci-1-err-ntc', 'sentinel-3-sl-1-rbt-ntc', 'sentinel-3-olci-2-lrr-ntc', 'sentinel-3-sl-2-lst-nrt', 'sentinel-3-sl-1-rbt-nrt', 'sentinel-3-olci-2-lrr-nrt', 'sentinel-3-sl-2-frp-nrt', 'sentinel-6-p4-1b-nrt', 'sentinel-5p-l1-ra-bd2-offl', 'sentinel-3-sr-1-sra-a-nrt', 'sentinel-6-p4-1b-ntc', 'sentinel-5p-l1-ra-bd2-nrti', 'sentinel-3-sr-1-sra-a-stc', 'sentinel-6-p4-1b-stc', 'sentinel-5p-l1-ra-bd1-rpro', 'sentinel-3-sr-1-sra-a-ntc', 'sentinel-6-p4-2-nrt', 'sentinel-5p-l1-ra-bd5-nrti', 'sentinel-1-grd', 'sentinel-6-p4-2-ntc', 'sentinel-5p-l1-ra-bd6-offl', 'sentinel-6-p4-2-stc', 'ccm-dem', 'sentinel-5p-l1-ra-bd4-nrti', 'sentinel-3-olci-2-wrr-ntc', 'sentinel-3-sr-1-sra-stc', 'ccm-sar', 'sentinel-5p-l1-ra-bd1-nrti', 'sentinel-3-sr-1-sra-ntc', 'ccm-optical', 'sentinel-5p-l1-ra-bd3-nrti', 'sentinel-3-sr-1-sra-nrt', 'sentinel-5p-l1-ra-bd6-rpro', 'sentinel-3-sr-2-lan-hy-nrt', 'opengeohub-landsat-bimonthly-mosaic-v1.0.1', 'sentinel-5p-l1-ra-bd5-rpro', 'sentinel-3-sr-2-lan-hy-stc', 'sentinel-6-amr-c-nrt', 'sentinel-5p-l1-ra-bd8-rpro', 'sentinel-3-sr-2-lan-hy-ntc', 'sentinel-6-amr-c-ntc', 'sentinel-5p-l1-ra-bd3-rpro', 'sentinel-3-sr-2-lan-li-ntc', 'sentinel-6-amr-c-stc', 'sentinel-5p-l1-ra-bd7-offl', 'sentinel-3-sr-2-lan-li-stc', 'sentinel-5p-l1-ra-bd6-nrti', 'sentinel-3-sr-2-lan-li-nrt', 'sentinel-5p-l1-ra-bd5-offl', 'sentinel-3-sr-2-lan-si-stc', 'sentinel-5p-l1-ra-bd8-offl', 'sentinel-3-sr-2-lan-si-nrt', 'sentinel-5p-l1-ra-bd2-rpro', 'sentinel-3-sr-2-lan-si-ntc', 'sentinel-5p-l1-ra-bd4-offl', 'sentinel-3-sr-2-lan-nrt', 'sentinel-5p-l1-ra-bd7-rpro', 'sentinel-3-sr-2-lan-stc', 'sentinel-5p-l1-ra-bd8-nrti', 'sentinel-3-sr-2-lan-ntc', 'sentinel-5p-l1-ra-bd4-rpro', 'sentinel-3-sr-2-wat-stc', 'sentinel-5p-l1-ra-bd7-nrti', 'sentinel-3-sr-2-wat-nrt', 'sentinel-5p-l1-ra-bd3-offl', 'sentinel-3-sr-2-wat-ntc', 'sentinel-5p-l1-ra-bd1-offl', 'sentinel-3-syn-2-aod-ntc', 'sentinel-5p-l2-aer-ai-nrti', 'sentinel-3-syn-2-syn-stc', 'sentinel-5p-l2-aer-ai-offl', 'sentinel-3-syn-2-syn-ntc', 'sentinel-5p-l2-aer-ai-rpro', 'sentinel-3-syn-2-v10-ntc', 'sentinel-5p-l2-aer-lh-nrti', 'sentinel-3-syn-2-v10-stc', 'sentinel-5p-l2-aer-lh-offl', 'sentinel-3-syn-2-vg1-ntc', 'sentinel-5p-l2-aer-lh-rpro', 'sentinel-3-syn-2-vg1-stc', 'sentinel-5p-l2-ch4-offl', 'sentinel-3-syn-2-vgp-stc', 'sentinel-5p-l2-ch4-rpro', 'sentinel-3-syn-2-vgp-ntc', 'sentinel-5p-l2-cloud-nrti', 'sentinel-5p-l2-cloud-offl', 'sentinel-5p-l2-cloud-rpro', 'sentinel-5p-l2-co-rpro', 'sentinel-5p-l2-co-nrti', 'sentinel-5p-l2-co-offl', 'sentinel-5p-l2-hcho-nrti', 'sentinel-5p-l2-hcho-rpro', 'sentinel-5p-l2-hcho-offl', 'sentinel-5p-l2-no2-offl', 'sentinel-5p-l2-no2-rpro', 'sentinel-5p-l2-no2-nrti', 'sentinel-5p-l2-np-bd3-rpro', 'sentinel-5p-l2-np-bd3-offl', 'sentinel-5p-l2-np-bd6-rpro', 'sentinel-5p-l2-np-bd6-offl', 'sentinel-5p-l2-np-bd7-rpro', 'sentinel-5p-l2-np-bd7-offl', 'sentinel-5p-l2-o3-pr-nrti', 'sentinel-5p-l2-o3-pr-rpro', 'sentinel-5p-l2-o3-pr-offl', 'sentinel-5p-l2-o3-tcl-offl', 'sentinel-5p-l2-o3-tcl-nrti', 'sentinel-5p-l2-o3-tcl-rpro', 'sentinel-5p-l2-o3-rpro', 'sentinel-5p-l2-o3-nrti', 'sentinel-5p-l2-o3-offl', 'sentinel-5p-l2-so2-offl', 'sentinel-5p-l2-so2-nrti', 'sentinel-5p-l2-so2-rpro']

# define a bounding box or area of interest
bbox = [11.439263980173, 47.81384831137465, 11.475186504136818, 47.83147903246192] # this location is in Germany
# search for items in sentinel-2-l2a collection
items = catalog.search(
    collections=["sentinel-2-l2a"],
    bbox=bbox,
    datetime="2023-05-01/2023-09-30",
    query={"eo:cloud_cover": {"lt": 20}},
).item_collection()

# check how many items found
print(f'Number of items found: {len(items)}')
# get the first image id as example
items = list(items)[:1]
# print the properties of the first item
items[0].properties # you can use this info to filter queries

# read the image data from the STAC items
data = stac_load(
    items,
    chunks={"x": 1024, "y": 1024},
    resolution=10,
    bbox=bbox, # you can also use selected band names here. here we load all bands
) # always return xarray.Dataset
data

# load B02_10m band only
blue = data['B02_10m'].squeeze()
blue = blue.compute() # load data into memory