• »
  • NOWATORSKA TECHNOLOGIA »
  • Wykorzystanie obrazów Sentinel-2 do monitorowania ujścia rzeki Syr Daria i Jeziora Aralskiego oraz analizowanie ich za pomocą wskaźnika NDWI na NSIS

Wykorzystanie obrazów Sentinel-2 do monitorowania ujścia rzeki Syr Daria i Jeziora Aralskiego oraz analizowanie ich za pomocą wskaźnika NDWI na NSIS

W tym artykule przedstawimy prosty przykład, w jaki sposób można łatwo uzyskać dostęp do danych z NSIS za pomocą interfejsu API RESTO. Uzyskamy dostęp do obrazów Sentinel-2 zawierających dane dla rzeki Syr Daria i Jeziora Aralskiego w Kazachstanie. Naszym celem będzie monitorowanie dziennej ilości wody w pobliżu tego jeziora od maja do lipca 2023 roku.

Użyjemy pakietów Python zarówno do uzyskiwania obrazów z chmury, jak i analizowania obecności wody za pomocą znormalizowanego różnicowego wskaźnika wody (NDWI).

Czym jest NDWI?

NDWI to znormalizowany wskaźnik różnicowy obliczany na podstawie współczynnika odbicia światła w dwóch różnych pasmach spektralnych: zielonym (G) i bliskiej podczerwieni (NIR). Wzór na obliczenie NDWI jest następujący:

../_images/nwdi_formula.png

gdzie:

  • G - współczynnik odbicia w paśmie zielonym,

  • NIR - współczynnik odbicia w paśmie bliskiej podczerwieni.

Wartości NDWI wahają się od -1 do 1, gdzie wyższe wartości wskazują na obecność wody, podczas gdy niższe wartości mogą wskazywać na inne typy pokrycia terenu, takie jak roślinność lub obszary miejskie.

Zalety korzystania ze wskaźnika NDWI

Obrazy rastrowe zawierające wskaźnik NDWI mogą być wykorzystywane w różnych dziedzinach, na przykład takich jak:

Monitorowanie zasobów wodnych

Śledzenie zmian poziomu wody w jeziorach, rzekach i zbiornikach, szczególnie ważne na obszarach dotkniętych suszą lub nadmiernymi opadami deszczu.

Zarządzanie kryzysowe

Szybka identyfikacja zalanych obszarów podczas powodzi i ocena szkód.

Ochrona środowiska

Monitorowanie terenów podmokłych i innych ekosystemów wodnych w celu zachowania ich różnorodności biologicznej.

Planowanie urbanistyczne

Analiza wpływu rozwoju miast na pobliskie zasoby wodne.

Rolnictwo

Ocena nawadniania pól uprawnych i monitorowanie stanu roślin.

Co zostanie omówione?

  • Definiowanie zmiennych dla przechowywania danych

  • Praca z interfejsem API RESTO

  • Tworzenie szablonu dla adresu URL API

  • Tworzenie żądania z szablonu

  • Uzyskiwanie ścieżek produktów

  • Uzyskanie NDWI

  • Funkcje odczytu, obliczania i zapisywania danych rastrowych

  • Obliczanie rastrów NWDI

  • Analiza

Wymagania wstępne

Nr 1 Hosting

Jest wymagane konto hostingowe NSIS z dostępem do interfejsu Horizon: https://horizon.cloudferro.com.

Nr 2 Działający Jupyter Notebook

Kod przedstawiony w tym artykule działa w środowisku Jupyter Notebook i może zostać zainstalowany na wybranej platformie, postępując zgodnie z oficjalną stroną instalacji Jupyter Notebook.

Jednym ze szczególnych sposobów instalacji (i w żadnym wypadku nie jest to wymagane do śledzenia tego artykułu) jest zainstalowanie go na klastrze Kubernetes. W takim przypadku będzie przydatny artykuł:

Możesz pobrać cały kod tego artykułu i uruchomić go w swoim Jupyter Notebook, korzystając z tego linku:

sentinel2_aral_lake.ipynb

Nr 3 Wymagany działający pakiet Python

Aby połączyć się z RESTO, użyjemy tylko jednego pakietu Pythona - requests. Jest to standardowa biblioteka, która może wysyłać żądania GET/POST do serwerów HTTP i obsługiwać odpowiedzi.

Jeśli pakiet nie jest jeszcze zainstalowany, zainstaluj go, używając polecenia

pip install requests

a następnie zaimportuj go za pomocą następującego polecenia, gdzieś na początku pliku:

import requests

Nr 4 Zainstalowane biblioteki Python dla danych obrazowych

Pakiety Pythona zorientowane przestrzennie, których będziemy używać do danych obrazowych, są następujące:

rasterio

Umożliwia odczytywanie danych rastrowych i wizualizację utworzonych obrazów NDWI.

osgeo

API GDAL dla Pythona.

numpy

Popularny pakiet Pythona, umożliwiający tworzenie złożonych równań matematycznych i formatów, takich jak macierze.

os

Łączy Pythona z systemem operacyjnym.

Instalacja potrzebnych pakietów może być skomplikowana (dotyczy to zwłaszcza**osgeo**), dlatego warto skorzystać ze środowiska conda. Jest to open-source’owy menedżer pakietów i środowisk, współdzielony przez wiele pakietów geoprzestrzennych.

Definiowanie zmiennych dla przechowywania danych

Interesuje nas okres, w którym obserwujemy i analizujemy ilość wody w rzece Syr Darya i jeziorze Aralskim, a także odpowiadający im obszar. Współrzędne są definiowane za pomocą BBOX (Bounding Box) jako lista współrzędnych - Xmin, Ymin, Xmax, Ymax. Wszystkie współrzędne zostaną zdefiniowane w EPSG:4326.

W Pythonie definiujemy trzy zmienne:

# Timerange of data that we want to recieve
start_date = '2023-05-01'
end_date = '2023-07-31'
# Geometry in form of a BBOX
bbox = [60.2347,46.1731,61.1160,46.6541]

W wybranych notatnikach Jupyter będzie to zazwyczaj wyglądać tak:

../_images/the_look_of_jupyter_book.png

Praca z interfejsem API RESTO

API RESTO to sposób komunikacji z EODATA, które jest ogromnym repozytorium danych obserwacji Ziemi. RESTO jest interfejsem API opartym na RESTful API. Umożliwia ono użytkownikom obsługiwać żądania danych na podstawie parametrów przestrzennych i czasowych.

Tworzenie szablonu dla adresu URL API

Będziemy używać tego predefiniowanego szablonu do obsługi żądań danych. Niektóre parametry są zdefiniowane jako zmienne i będzie można je edytować.

API_URL = """http://catalogue.nsiscloud.polsa.gov.pl/resto/api/collections/Sentinel2/search.json?\
box={minx},{miny},{maxx},{maxy}\
&startDate={start_date}\
&completionDate={end_date}\
&sortParam=startDate\
&sortOrder=ascending\
&platform=S2A\
&cloudCover=[0,10]"""

Tworzenie żądania na podstawie szablonu

Teraz, gdy mamy już szablon żądania i wartości jego parametrów, możemy połączyć je w rzeczywiste polecenie:

specified_url = API_URL.format(minx=bbox[0], miny=bbox[1], maxx=bbox[2], maxy=bbox[3], start_date=start_date, end_date=end_date)

Wartość uzyskanego adresu URL to:

specified_url

  `https://catalogue.nsiscloud.polsa.gov.pl/resto/api/collections/Sentinel2/search.json?box=60.2347,46.1731,61.116,46.6541&startDate=2023-05-01&completionDate=2023-09-30&sortParam=startDate&sortOrder=ascending&platform=S2A&cloudCover=[0,10] <https://catalogue.nsiscloud.polsa.gov.pl/resto/api/collections/Sentinel2/search.json?box=60.2347,46.1731,61.116,46.6541&startDate=2023-05-01&completionDate=2023-09-30&sortParam=startDate&sortOrder=ascending&platform=S2A&cloudCover=[0,10]>`_

Po zdefiniowaniu żądania, wyślijmy je do API.

response = requests.get(specified_url).json()

Poniżej znajduje się wynik dla pierwszego znalezionego produktu.

{'type': 'Feature',
 'id': '259ff692-e36b-45d8-8baa-a3e49fabb5ce',
 'geometry': {'type': 'Polygon',
  'coordinates': [[[60.2782172365724, 45.9163137717613],
    [60.9930021665452, 45.8957352519379],
    [61.0656571326958, 46.8814596616377],
    [60.6234129381189, 46.8944007194611],
    [60.6119380882436, 46.8622024095855],
    [60.5598328017783, 46.7158712476585],
    [60.5079224466574, 46.5694873249003],
    [60.456176385744, 46.4230015400117],
    [60.4045426003916, 46.2764467160812],
    [60.3529921585166, 46.1298520112078],
    [60.3015821350397, 45.9831749949217],
    [60.2782172365724, 45.9163137717613]]]},
 'properties': {'collection': 'SENTINEL-2',
  'status': 'ONLINE',
  'license': {'licenseId': 'unlicensed',
   'hasToBeSigned': 'never',
   'grantedCountries': None,
   'grantedOrganizationCountries': None,
   'grantedFlags': None,
   'viewService': 'public',
   'signatureQuota': -1,
   'description': {'shortName': 'No license'}},
  'parentIdentifier': None,
  'title': 'S2A_MSIL2A_20230503T064621_N0509_R020_T40TGS_20230503T110902.SAFE',
  'description': 'The Copernicus Sentinel-2 mission consists of two polar-orbiting satellites that are positioned in the same sun-synchronous orbit, with a phase difference of 180°. It aims to monitor changes in land surface conditions. The satellites have a wide swath width (290 km) and a high revisit time. Sentinel-2 is equipped with an optical instrument payload that samples 13 spectral bands: four bands at 10 m, six bands at 20 m and three bands at 60 m spatial resolution [https://dataspace.copernicus.eu/explore-data/data-collections/sentinel-data/sentinel-2].',
  'organisationName': None,
  'startDate': '2023-05-03T06:46:21.024Z',
  'completionDate': '2023-05-03T06:46:21.024Z',
  'productType': 'S2MSI2A',
  'processingLevel': 'S2MSI2A',
  'platform': 'S2A',
  'instrument': 'MSI',
  'resolution': 0,
  'sensorMode': None,
  'orbitNumber': 41058,
  'quicklook': None,
  'thumbnail': 'https://catalogue.nsiscloud.polsa.gov.pl/get-object?path=/Sentinel-2/MSI/L2A/2023/05/03/S2A_MSIL2A_20230503T064621_N0509_R020_T40TGS_20230503T110902.SAFE/S2A_MSIL2A_20230503T064621_N0509_R020_T40TGS_20230503T110902-ql.jpg',
  'updated': '2024-06-18T12:16:41.475Z',
  'published': '2023-05-03T13:18:52.861Z',
  'snowCover': 0,
  'cloudCover': 9.513979,
  'gmlgeometry': '<gml:Polygon srsName="EPSG:4326"><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>60.2782172365724,45.9163137717613 60.9930021665452,45.8957352519379 61.0656571326958,46.8814596616377 60.6234129381189,46.8944007194611 60.6119380882436,46.8622024095855 60.5598328017783,46.7158712476585 60.5079224466574,46.5694873249003 60.456176385744,46.4230015400117 60.4045426003916,46.2764467160812 60.3529921585166,46.1298520112078 60.3015821350397,45.9831749949217 60.2782172365724,45.9163137717613</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs></gml:Polygon>',
  'centroid': {'type': 'Point',
   'coordinates': [60.7315835400736, 46.3584900636335]},
  'productIdentifier': '/eodata/Sentinel-2/MSI/L2A/2023/05/03/S2A_MSIL2A_20230503T064621_N0509_R020_T40TGS_20230503T110902.SAFE',
  'orbitDirection': None,
  'timeliness': None,
  'relativeOrbitNumber': 20,
  'processingBaseline': 5.09,
  'missionTakeId': 'GS2A_20230503T064621_041058_N05.09',
  'services': {'download': {'url': 'https://catalogue.nsiscloud.polsa.gov.pl/download/259ff692-e36b-45d8-8baa-a3e49fabb5ce',
    'mimeType': 'application/octet-stream',
    'size': 464127649}},
  'links': [{'rel': 'self',
    'type': 'application/json',
    'title': 'GeoJSON link for 259ff692-e36b-45d8-8baa-a3e49fabb5ce',
    'href': 'https://catalogue.nsiscloud.polsa.gov.pl/resto/collections/SENTINEL-2/259ff692-e36b-45d8-8baa-a3e49fabb5ce.json'}]}}

Uzyskiwanie ścieżek produktów

Po otrzymaniu odpowiedzi od RESTO możemy pozyskać ścieżki do katalogu EODATA dla tych produktów. Odczytamy je z parametru o nazwie productIdentifier:

dataset = [x['properties']['productIdentifier'] for x in response['features']]

Tutaj możesz zobaczyć ścieżki dla swoich produktów:

dataset

['/eodata/Sentinel-2/MSI/L2A/2023/05/03/S2A_MSIL2A_20230503T064621_N0509_R020_T40TGS_20230503T110902.SAFE',
 '/eodata/Sentinel-2/MSI/L2A/2023/05/26/S2A_MSIL2A_20230526T065621_N0509_R063_T41TLM_20230526T104855.SAFE',
 '/eodata/Sentinel-2/MSI/L1C/2023/05/26/S2A_MSIL1C_20230526T065621_N0509_R063_T41TLM_20230526T084740.SAFE',
 '/eodata/Sentinel-2/MSI/L1C/2023/05/26/S2A_MSIL1C_20230526T065621_N0509_R063_T40TGS_20230526T084740.SAFE',
 '/eodata/Sentinel-2/MSI/L2A/2023/05/26/S2A_MSIL2A_20230526T065621_N0509_R063_T40TGS_20230526T104855.SAFE',
 '/eodata/Sentinel-2/MSI/L2A/2023/06/05/S2A_MSIL2A_20230605T065621_N0509_R063_T40TGS_20230605T104755.SAFE',
 '/eodata/Sentinel-2/MSI/L1C/2023/06/05/S2A_MSIL1C_20230605T065621_N0509_R063_T40TGS_20230605T084725.SAFE',
 '/eodata/Sentinel-2/MSI/L1C/2023/06/05/S2A_MSIL1C_20230605T065621_N0509_R063_T41TLM_20230605T084725.SAFE',
 '/eodata/Sentinel-2/MSI/L2A/2023/06/05/S2A_MSIL2A_20230605T065621_N0509_R063_T41TLM_20230605T104755.SAFE',
 '/eodata/Sentinel-2/MSI/L1C/2023/06/15/S2A_MSIL1C_20230615T065621_N0509_R063_T41TLM_20230615T092752.SAFE',
 '/eodata/Sentinel-2/MSI/L2A/2023/06/15/S2A_MSIL2A_20230615T065621_N0509_R063_T41TLM_20230615T110856.SAFE',
 '/eodata/Sentinel-2/MSI/L1C/2023/06/15/S2A_MSIL1C_20230615T065621_N0509_R063_T40TGS_20230615T092752.SAFE',
 '/eodata/Sentinel-2/MSI/L1C/2023/07/02/S2A_MSIL1C_20230702T064631_N0509_R020_T40TGS_20230702T083747.SAFE',
 '/eodata/Sentinel-2/MSI/L2A/2023/07/02/S2A_MSIL2A_20230702T064631_N0509_R020_T40TGS_20230702T105455.SAFE',
 '/eodata/Sentinel-2/MSI/L2A/2023/07/05/S2A_MSIL2A_20230705T065631_N0509_R063_T41TLM_20230705T105352.SAFE',
 '/eodata/Sentinel-2/MSI/L2A/2023/07/05/S2A_MSIL2A_20230705T065631_N0509_R063_T40TGS_20230705T105352.SAFE',
 '/eodata/Sentinel-2/MSI/L1C/2023/07/05/S2A_MSIL1C_20230705T065631_N0509_R063_T40TGS_20230705T084704.SAFE',
 '/eodata/Sentinel-2/MSI/L1C/2023/07/05/S2A_MSIL1C_20230705T065631_N0509_R063_T41TLM_20230705T084704.SAFE',
 '/eodata/Sentinel-2/MSI/L1C/2023/07/15/S2A_MSIL1C_20230715T065631_N0509_R063_T40TGS_20230715T084624.SAFE',
 '/eodata/Sentinel-2/MSI/L2A/2023/07/15/S2A_MSIL2A_20230715T065631_N0509_R063_T40TGS_20230715T104959.SAFE']

Pozyskiwanie NDWI

Wykorzystajmy teraz dane dla rzeki Syr Daria i Jeziora Aralskiego. Najpierw należy

  • obliczyć wskaźnik NDWI dla każdego z pikseli, a następnie

  • utworzyć z niego raster.

Korzystając ze wszystkich znalezionych elementów, będziemy mogli monitorować stan jeziora przez cały okres.

Tutaj użyjemy uprzednio zainstalowanych pakietów:

import rasterio
from osgeo import gdal, gdal_array, osr
import numpy as np
import os

Funkcje odczytu, obliczania i zapisywania danych rastrowych

Tutaj przedstawiamy kilka funkcji do odczytu danych rastrowych w macierzy numpy, obliczania NDWI z macierzami NIR i GREEN oraz zapisywania wyniku jako nowego rastra. Przeprowadzimy takie obliczenia dla każdego pobranego elementu. Ostatecznie uzyskamy dane NDWI dla rzeki Syr Daria i Jeziora Aralskiego z okresu trzech miesięcy, od maja do lipca 2023 roku.

def getFullPath(dir: str, resolution: int, band: str):
    if not os.path.isdir(dir):
        raise ValueError(f"Provided path does not exist: {dir}")
    elif resolution not in [10,20,60]:
        raise ValueError(f"Provided resolution does not exist: R{resolution}m")
    else:
        full_path = dir
        while True:
            content = os.listdir(full_path)
            if len(content) == 0:
                raise ValueError(f"Directory empty: {full_path}")
            elif len(content) == 1:
                if full_path[-1] != '/':
                    full_path = full_path + '/' + content[0]
                else:
                    full_path = full_path + content[0]
            else:
                if 'GRANULE' in content:
                    full_path = full_path + '/' + 'GRANULE'
                    break
                else:
                    raise ValueError(f"Unsupported dir architecture: {full_path}")
        full_path = full_path + '/' + os.listdir(full_path)[0]
        full_path = full_path + '/' + "IMG_DATA"
        if len(os.listdir(full_path)) == 3:
            full_path = full_path + '/' + f'R{resolution}m'
            images = os.listdir(full_path)
            for img in images:
                if band in img:
                    return full_path + '/' + img
            raise ValueError(f'No such band {band} in directory: {full_path}')
        else:
            images = os.listdir(full_path)
            for img in images:
                if band in img:
                    return full_path + '/' + img
            raise ValueError(f'No such band {band} in directory: {full_path}')

# Get transformation matrix from raster
def getTransform(pathToRaster):
    dataset = gdal.Open(pathToRaster)
    transformation = dataset.GetGeoTransform()
    return transformation

# Read raster and return pixels' values matrix as int16, new transformation matrix, crs
def readRaster(path, resolution, band):
    path = getFullPath(path, resolution, band)
    trans = getTransform(path) # trzeba zdefiniować który kanał
    raster, crs = rasterToMatrix(path)
    return raster.astype(np.int16), crs, trans

def rasterToMatrix(pathToRaster):
    with rasterio.open(pathToRaster) as src:
        matrix = src.read(1)
    return matrix, src.crs.to_epsg()

# Transform numpy's matrix to geotiff; pass new raster's filepath, matrix with pixels' values, gdal file type, transformation matrix, projection, nodata value
def npMatrixToGeotiff(filepath, matrix, gdalType, projection, transformMatrix, nodata = None):
    driver = gdal.GetDriverByName('Gtiff')
    if len(matrix.shape) > 2:
        (bandNr, yRes, xRes) = matrix.shape
        image = driver.Create(filepath, xRes, yRes, bandNr, gdalType)
        for b in range(bandNr):
            b = b + 1
            band = image.GetRasterBand(b)
            if nodata is not None:
                band.SetNoDataValue(nodata)
            band.WriteArray(matrix[b-1,:,:])
            band.FlushCache
    else:
        bandNr = 1
        (yRes, xRes) = matrix.shape
        image = driver.Create(filepath, xRes, yRes, bandNr, gdalType)
        print(type(image))
        band = image.GetRasterBand(bandNr)
        if nodata is not None:
            band.SetNoDataValue(nodata)
        band.WriteArray(matrix)
        band.FlushCache
    image.SetGeoTransform(transformMatrix)
    image.SetProjection(projection)
    del driver, image, band

Obliczanie rastrów NWDI

Teraz użyjemy zdefiniowanych powyżej funkcji do wygenerowania rastrów NWDI. Jedynymi danymi potrzebnymi w tym kroku jest lista ścieżek do produktów (wyodrębnionych z archiwum zip). Funkcja readRaster wybierze określone pasmo z podanej ścieżki.

# Path to catalog for processed images with NDWI index
compution_output = '/dir/for/output'

# Iterating over single product
for item in dataset:
    # Reading name from path
    name = item.split('/')[-1]
    # Reading green band into matrix
    green = readRaster(item, 10, 'B03')
    # Reading NIR band into matrix
    nir = readRaster(item, 10, 'B08')
    # Calculating NDWI matrix
    ndwi = (green[0]-nir[0]) / (green[0]+nir[0])
    # Creating SpatialReference object and setting it to match original's raster CRS
    projection = osr.SpatialReference()
    projection.ImportFromEPSG(green[1])
    # Creating raster from matrix in GeoTiff format
    npMatrixToGeotiff(f'{compution_output}/{name}.tif', ndwi, gdal_array.NumericTypeCodeToGDALTypeCode(np.float32), projection.ExportToWkt(), green[2])

Po pomyślnym utworzeniu i zapisaniu nowych obrazów możemy zwizualizować je w Pythonie za pomocą pakietu rasterio.

img = rasterio.open('/path/to/one/of/the/new/rasters/with/NDWI')
from rasterio.plot import show
show(img)

Analiza

Teraz, gdy otrzymaliśmy obrazy, zobaczmy, czego możemy się z nich dowiedzieć. Oto trzy obrazy dla maja, czerwca i lipca 2023 roku:

../_images/the_end_of_the_lake.png

W tym stylu obrazów

  • kolor ciemnoniebieski oznacza wysokie prawdopodobieństwo obecności wody, podczas gdy

  • kolor żółty oznacza całkowity brak wody.

Jak widać, ujście rzeki zaczyna być coraz bardziej żółte z każdym miesiącem. Wnioskujemy, że w pobliżu ujścia rzeki Syr Daria prawie nie ma wody. Możemy założyć, że latem rzeki tracą dużo wody. Może się tak dziać z powodu:

  • wysokich temperatur, ale także dlatego, że

  • ludzie mieszkający wzdłuż tej rzeki wykorzystują więcej wpdy do nawadniania swoich pól.

Niestety oznacza to, że Jezioro Aralskie nadal będzie umierać i nie jest wykluczone, że będziemy świadkami jego końca – końca tego, co kiedyś było jednym z największych jezior na naszej planecie.

Co można zrobić dalej?

Chcesz zobaczyć dodatkowe przykłady badania przyrody za pomocą zdjęć satelitarnych?

Rozrost miast:

Monitorowanie rozrostu miast za pomocą danych Sentinel-1 SAR z wykorzystaniem wartości pikseli i progowania polaryzacji w NSIS

Powodzie:

Analiza i monitorowanie powodzi przy użyciu języka Python i zdjęć satelitarnych Sentinel-2 na NSIS