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

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

W tym artykule wykorzystasz techniki progowania i porównasz wartości polaryzacji VV i VH do analizy obszarów miejskich. Jako ilustrację, spróbujesz pobrać obrazy Sentinel-1 zawierające dane obszaru miejskiego Warszawy.

Analiza obrazów radarowych z syntetyczną aperturą (SAR) umożliwia wykrywanie zmian w odbiciach fal radarowych, co ma kluczowe znaczenie dla identyfikacji obszarów zurbanizowanych. Obrazy SAR wykazują wyraźne wzorce odbić, które pozwalają na wyznaczenie obszarów miejskich, w tym aglomeracji i mniejszych osiedli miejskich. Obszary zurbanizowane wyróżniają się na tych obrazach przede wszystkim ze względu na intensywne odbicie fal radarowych przez konstrukcje stworzone przez człowieka, drogi i inne sztuczne elementy infrastruktury, co ułatwia ich wyraźną identyfikację w przeciwieństwie do naturalnego terenu. Taka analiza obszarów takich jak aglomeracja warszawska może posłużyć do

  • monitorowania wzrostu i rozwoju miasta,

  • planowania urbanistycznego,

  • oceniania zagrożeń naturalnych i antropogenicznych oraz

  • zarządzania środowiskiem.

Oto obraz, który otrzymasz:

../_images/exampleimage2.png

Wyraźnie przedstawia on obszar aglomeracji warszawskiej przylegający do regionu między aglomeracjami warszawską i łódzką. Obecnie planuje się zagospodarowanie tej przestrzeni w celu połączenia obu aglomeracji. W centralnej części ramki planowane jest utworzenie dużego kompleksu lotniskowego.

Co zostanie omówione?

  • Dostęp do danych z NSIS Cloud

  • Wyodrębnianie obszarów miejskich przez progowanie

  • Zapewnienie analizy ilościowej poprzez analizę poszczególnych pikseli

Wymagania wstępne

  1. Konto NSIS Cloud

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

  1. Dostęp do S3

Upewnij się, że masz poświadczenia dostępu do S3. Do interakcji z magazynem NSIS Cloud S3 potrzebny będzie klucz dostępu i klucz tajny. Zobacz artykuł

Jak wygenerować poświadczenia EC2 i zarządzać nimi na NSIS Cloud

  1. Dostępność Jupyter Notebook

Kod w tym przewodniku działa w środowisku Jupyter Notebook. Można go zainstalować na wybranej platformie, klikając łącze: https://jupyter.org/install

Można także skonfigurować Jupyter Notebook na klastrze NSIS Kubernetes:

Instalacja JupyterHub na klastrze Magnum Kubernetes w NSIS Cloud

  1. Biblioteki Python zainstalowane w Jupyter Notebook

Kod w tym artykule opiera się na następujących bibliotekach w języku Python zainstalowanych w Jupyter Notebook:

json  requests os boto3 datetime
numpy rasterio matplotlib zipfile io

Możesz pobrać kod z tego artykułu i uruchomić go bezpośrednio w swojej instalacji JupyterLab:

urban_sprawl.ipynb

Krok 1 Uzyskanie dostępu do danych z NSIS Cloud

1. Import i deklaracja zmiennych

import json
import requests
import os
import boto3
import datetime
import numpy as np
import rasterio
from rasterio.transform import Affine
import matplotlib.pyplot as plt
import zipfile
import io

# Settings for the API query
collection = 'SENTINEL-1'
product_type = 'GRD'
start_date = '2024-05-15T00:00:00.000Z'
end_date = '2024-05-17T00:05:00.000Z'
polygon_coords = '20.851215 52.261859, 21.271151 52.261859, 21.271151
52.017798, 20.851215 52.017798, 20.851215 52.261859'

# Creating URL with query parameters
url = f"https://catalogue.nsiscloud.polsa.gov.pl/odata/v1/Products?$filter=Collection/Name eq '{collection}' and Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'productType' and att
OData.CSC.StringAttribute/Value eq '{product_type}') and ContentDate/Start gt {start_date} and ContentDate/End lt {end_date} and
OData.CSC.Intersects(area=geography'SRID=4326;POLYGON(({polygon_coords}))')&$expand=Attributes&$top=1"

2. Konfiguracja S3 i funkcja tworzenia listy plików TIFF

# S3 access keys
access_key = 'your_access_key'
secret_key = 'your_secret_key'
host = 'http://eodata.cloudferro.com'

# Initialize the S3 client with the custom endpoint
s3 = boto3.client('s3',
                  aws_access_key_id=access_key,
                  aws_secret_access_key=secret_key,
                  endpoint_url=host)

# Function to list TIFF files in the S3 bucket
def list_tiff_files_in_s3(prefix):
    tiff_files = []
    response = s3.list_objects_v2(Bucket='DIAS', Prefix=prefix)

    for obj in response.get('Contents', []):
        if obj['Key'].lower().endswith('.tiff') or obj['Key'].lower().endswith('.tif'):
            tiff_files.append(obj['Key'])

    return tiff_files

3. Pobieranie ścieżki produktu z odpowiedzi API

# Fetch product path from the API response
response = requests.get(url)
products = json.loads(response.text)

# Extract the S3 path from the first product in the response
product_path = ''
for item in products.get('value', []):
   product_path = item.get('S3Path')
   break

# Replace '/eodata/' in the product_path if needed
product_path = product_path.replace('/eodata/', '') + '/'

4. Pobieranie plików TIFF z S3 i zapisywanie ich lokalnie

# List TIFF files in the S3 bucket
output_dir = os.path.join(os.getcwd(), 'sen_1')
if not os.path.exists(output_dir):
   os.makedirs(output_dir)

# Create directory if it doesn't exist
tiff_files = list_tiff_files_in_s3(product_path)

# Download TIFF files
for file in tiff_files:
file_name = os.path.basename(file)
   local_file_path = os.path.join(output_dir, file_name)
   s3.download_file('DIAS', file, local_file_path)
   print(f'Downloaded: {file_name}')

Downloaded: s1a-iw-grd-vh-20240516t162754-20240516t162819-053899-068d12-002.tiff
Downloaded: s1a-iw-grd-vv-20240516t162754-20240516t162819-053899-068d12-001.tiff

Krok 2 Wyodrębnianie obszarów miejskich przez progowanie

1. Odczytywanie i wstępne przetwarzanie danych

Pierwszym krokiem w analizie będzie odczyt danych. Dane zostaną odczytane przy użyciu biblioteki rasterio ze ścieżki folderu wyjściowego podanej w poprzedniej sekcji.

# Function to find files containing specific substrings ("vv" and "vh") in their names
def find_files_with_names(folder_path):
   vv_path = None
   vh_path = None

   for root, dirs, files in os.walk(folder_path):
      print(files)
      for file_name in files:
         if "vv" in file_name:
            vv_path = os.path.join(root, file_name)
         if "vh" in file_name:
            vh_path = os.path.join(root, file_name)
   return vv_path, vh_path
output_dir = 'sen_1'

# Find the "vv" and "vh" files in the specified directory
input_file_vv, input_file_vh = find_files_with_names(output_dir)
print(input_file_vv)

# Open the "vv" file and read its data and profile
with rasterio.open(input_file_vv) as src_vv:
   vv_band = src_vv.read(1)
   profile_vv = src_vv.profile

# Open the "vh" file and read its data and profile
with rasterio.open(input_file_vh) as src_vh:
   vh_band = src_vh.read(1)
   profile_vh = src_vh.profile

 # Set CRS and transform manually if they are missing
 crs = 'EPSG:4326'
 transform = Affine.translation(0.0, 0.0) * Affine.scale(30.0, -30.0)
 profile_vv.update(crs=crs, transform=transform)
 profile_vh.update(crs=crs, transform=transform)

['s1a-iw-grd-vv-20240516t162754-20240516t162819-053899-068d12-001.tiff',
's1a-iw-grd-vh-20240516t162754-20240516t162819-053899-068d12-002.tiff']
sen_1/s1a-iw-grd-vv-20240516t162754-20240516t162819-053899-068d12-001.tiff

2. Konwersja wartości pikseli na decybele

Po odczytaniu danych możemy przekonwertować wartości pikseli na decybele. Dane z Sentinel-1 są konwertowane na skalę decybelową (dB), aby lepiej wizualizować różnice w intensywności odbicia fal radarowych, co pomaga w analizie obszarów miejskich ze względu na ich specyficzne właściwości odbicia.

Najwyższy szczyt histogramu polaryzacji VV znajduje się w zakresie 15–20 dB, co sugeruje, że większość powierzchni na badanym obszarze ma umiarkowane rozproszenie radarowe. Próg 23,5 dB został wybrany w celu wyróżnienia pikseli o najwyższych wartościach rozproszenia radarowego, które mogą odpowiadać obszarom z intensywną infrastrukturą miejską lub innymi obiektami o wysokim współczynniku odbicia.

Najwyższy szczyt histogramu polaryzacji VH znajduje się w zakresie 13–18  B, co sugeruje, że kanał ten jest bardziej wrażliwy na powierzchnie o niższym rozproszeniu radarowym. Próg 21,5 dB dla kanału VH jest odpowiedni do identyfikacji pikseli o wyższym rozproszeniu wstecznym, co pomaga lepiej rozróżnić struktury terenu o nieco niższym, ale wciąż znaczącym rozproszeniu.

def convert_to_db(image):
    # Create a copy of the image to avoid modifying the original array
    image = np.copy(image)

    # Replace non-positive values with a small positive value to avoid taking log of zero or negative
    image[image <= 0] = np.finfo(float).eps

    # Calculate dB values
    with np.errstate(divide='ignore', invalid='ignore'):
        image_db = 10 * np.log10(image)

    # Replace infinite values with the maximum finite value
    max_value = np.nanmax(image_db[np.isfinite(image_db)])
    image_db[np.isinf(image_db)] = max_value

    return image_db

vv_band_db = convert_to_db(vv_band)
vh_band_db = convert_to_db(vh_band)

3. Analiza danych i progowanie

Następnym krokiem jest wykreślenie histogramów decybelowych wartości pikseli dla pasm VV i VH. Te histogramy zapewniają wgląd w rozkład intensywności odbić fal radarowych i pomagają w zrozumieniu charakterystyki obserwowanego obszaru. Wybrane progi (23,5 dB dla VV i 21,5 dB dla VH) zostały dobrane na podstawie wcześniejszej wiedzy lub wymagań analizy specyficznych dla obszarów miejskich w celu segmentacji obszarów zainteresowania na obrazach. Zastosowana tutaj metoda progowania konwertuje wartości pikseli powyżej progu na 1, a te poniżej na 0, co pomaga w wyznaczeniu cech lub obszarów o specyficznych charakterystykach odbicia fal radarowych typowych dla środowisk miejskich.

fig, axs = plt.subplots(1, 2, figsize=(18, 5))

# Histogram of VV Image in dB
axs[0].hist(vv_band_db.flatten(), bins=50, color='blue', alpha=0.7, label='VV dB')
axs[0].set_title('Histogram of VV Image in dB')
axs[0].set_xlabel('Radar Reflection Value (dB)')
axs[0].set_ylabel('Number of Pixels')
axs[0].legend()
axs[0].grid(True, linestyle='--', alpha=0.5)

# Histogram of VH Image in dB
axs[1].hist(vh_band_db.flatten(), bins=50, color='red', alpha=0.7, label='VH dB')
axs[1].set_title('Histogram of VH Image in dB')
axs[1].set_xlabel('Radar Reflection Value (dB)')
axs[1].set_ylabel('Number of Pixels')
axs[1].legend()
axs[1].grid(True, linestyle='--', alpha=0.5)

plt.show()

# Function to apply a threshold to an image
# Pixels greater than the threshold are set to 1, others to 0
def thresholding(image, threshold):
   return (image > threshold).astype(np.uint8)

# Define threshold values for VV and VH bands
threshold_value_vv = 23.5
threshold_value_vh = 21.5

# Apply thresholding to the VV and VH band
vv_band_threshold = thresholding(vv_band_db, threshold_value_vv)
vh_band_threshold = thresholding(vh_band_db, threshold_value_vh)

# Combine the thresholded images using a logical AND operation
# The result is 1 where both VV and VH are above their respective thresholds

combined_threshold = np.logical_and(vv_band_threshold, vh_band_threshold).astype(np.uint8)
../_images/histogram.png

4. Wykreślanie obrazów

Następnie na wykresach cząstkowych można wizualizować różnice w odbiciu sygnału radiowego dla różnych typów polaryzacji. Wyświetlając oryginalne i przekształcone decybelowo obrazy dla polaryzacji VV i VH, można zaobserwować wyraźne cechy i intensywność na obszarach miejskich. Zmiana odbicia sygnału na oryginalnych obrazach jest ledwo widoczna, co dowodzi, że stosowanie skali decybelowej jest przydatne. Co więcej, łącząc informacje z obu polaryzacji poprzez progowanie, zwiększ sięi zdolność do maskowania regionów zurbanizowanych, zapewniając bardziej kompleksowe zrozumienie cech miejskich obserwowanego obszaru.

# Create a figure with a 3x2 grid of subplots, setting the overall figure size
fig, ax = plt.subplots(3, 2, figsize=(18, 18))

# Display the original VV image
ax[0, 1].imshow(vv_band, cmap='viridis')
ax[0, 1].set_title('Original VV Image')

# Display the VV image in decibels
ax[1, 1].imshow(vv_band_db, cmap='viridis')
ax[1, 1].set_title('VV Image in dB')

# Display the VV image after thresholding
ax[2, 1].imshow(vv_band_threshold, cmap='inferno')
ax[2, 1].set_title('VV Image after Thresholding')

# Display the original VH image
ax[0, 0].imshow(vh_band, cmap='viridis')
ax[0, 0].set_title('Original VH Image')

# Display the VH image in decibels
ax[1, 0].imshow(vh_band_db, cmap='viridis')
ax[1, 0].set_title('VH Image in dB')

# Display the VH image after thresholding
ax[2, 0].imshow(vh_band_threshold, cmap='inferno')
ax[2, 0].set_title('VH Image after Thresholding')
# Turn off the axes for all subplots to improve the visual presentation
   for i in range(2):
      for j in range(3):
         ax[j, i].axis('off')
# Show the complete figure
plt.show()
../_images/exampleimages.png 
# Display the combined threshold image using a specific colormap
plt.figure(figsize=(18, 9))
im = plt.imshow(combined_threshold, cmap='twilight_shifted')
plt.title('Combined Threshold Image')
cbar = plt.colorbar(im)
cbar.set_label('Signal Strength')
plt.axis('off')
plt.show()
../_images/exampleimage2.png

Ta metoda analizy danych ułatwia proste maskowanie obszarów zurbanizowanych. Na przetworzonym obrazie aglomeracja warszawska jest wyraźnie widoczna po prawej stronie. Gwiaździste skupiska ludności są widoczne w Warszawie wzdłuż głównych korytarzy transportowych, takich jak linie kolejowe i drogi ekspresowe. Na obrazie są także widoczne mniejsze ośrodki miejskie. Ponadto można zaobserwować, że liczba maskowanych pikseli rośnie w miarę zbliżania się do centrum miasta, co wskazuje na gęstszą zabudowę miejską w śródmieściu w porównaniu do przedmieść

5. Analiza ilościowa

Po przetworzeniu zdjęć satelitarnych możemy przeprowadzić analizę ilościową. Charakterystyka danych satelitarnych pozwala analizować wartości poszczególnych pikseli, co jest kluczowe dla badań szczegółowych i statystycznych. W naszym przypadku, mając dane z progowaniem przy rozmiarze każdego piksela wynoszącym 10x10 metrów, możemy obliczyć procent i odpowiedni rozmiar obszaru zurbanizowanego.

# Calculate the number of urban pixels by summing up the binary values in the combined threshold image
urban_pixels = np.sum(combined_threshold)

# Define the area covered by a single pixel (in square meters)
pixel_area = 10 * 10 # Assuming each pixel represents a 10m x 10m area

# Calculate the total urban area in square meters
urban_area = urban_pixels * pixel_area

# Calculate the total number of pixels in the combined threshold image
total_pixels = combined_threshold.size

# Calculate the percentage of urbanized area in the image
urban_percentage = (urban_pixels / total_pixels) * 100

# Convert urban area from square meters to hectares
urban_area_hectares = urban_area / 10_000

# Convert total area from square meters to hectares
total_area_hectares = (total_pixels * pixel_area) / 10_000

print(f"Percentage of Urbanized Area: {urban_percentage:.2f}%")
print(f"Urbanized Area: {urban_area_hectares:.2f} ha")
print(f"Total Image Area: {total_area_hectares:.2f} ha")

Percentage of Urbanized Area: 9.08%
Urbanized Area: 401728.18 ha
Total Image Area: 4426150.50 ha

6. Pobieranie wyników

# Define output file names
output_file_db_vv = 'vv_db_urban_warsaw.tif'
output_file_db_vh = 'vh_db_urban_warsaw.tif'
output_file_threshold_combined = 'warsaw_urban_threshold_combined.tif'

# Update the data type in the profile to float32 for dB conversion
profile_vv.update(dtype=rasterio.float32)
profile_vh.update(dtype=rasterio.float32)

# Write the result data in dB to the output file
with rasterio.open(output_file_db_vv, 'w', **profile_vv) as dst:
   dst.write(vv_band_db.astype(np.float32), 1)

with rasterio.open(output_file_db_vh, 'w', **profile_vh) as dst:
   dst.write(vh_band_db.astype(np.float32), 1)

# Update the data type in the VV profile to uint8 for the combined
thresholded image
profile_vv.update(dtype=rasterio.uint8)

# Write the combined thresholded image data to the output file
with rasterio.open(output_file_threshold_combined, 'w', **profile_vv) as dst:
   dst.write(combined_threshold, 1)

Co można zrobić dalej?

Również interesujący jest artykuł:

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