5.2 Zpracování obrazových dat pomocí Rasterio a EODAG¶
Sentinel-2¶
Sentinel-2 je družicová mise v rámci programu Evropské unie Copernicus, jejímž cílem je monitorovat zemský povrch. Skládá se ze dvou identických družic Sentinel-2A a Sentinel-2B, které byly vypuštěny v roce 2015, resp. 2017. Společně poskytují multispektrální snímky zemského povrchu s vysokým rozlišením.
Data Sentinel-2 můžeme stáhnout např. pomocí Sentinel Hub (https://apps.sentinel-hub.com/eo-browser/) anebo Python balíčku eodag (https://eodag.readthedocs.io/).
Zájmové území (AOI)¶
Můžeme použít vrstvu obcí dostupnou z datasetu Data250 (dokumentace):
In [2]:
Copied!
import fiona
nuts_name = "Katovice"
with fiona.open('/mnt/repository/155UZPR/05/PolbndMunDA.shp') as nuts:
nuts_crs = nuts.crs
for feature in nuts:
if feature["properties"]["NAMN"] == nuts_name:
aoi = feature["geometry"]
break
nuts_crs, aoi.type
import fiona
nuts_name = "Katovice"
with fiona.open('/mnt/repository/155UZPR/05/PolbndMunDA.shp') as nuts:
nuts_crs = nuts.crs
for feature in nuts:
if feature["properties"]["NAMN"] == nuts_name:
aoi = feature["geometry"]
break
nuts_crs, aoi.type
Out[2]:
(CRS.from_epsg(5514), 'Polygon')
Zobrazení AOI¶
Na tomto místě můžeme použít knihovnu GeoPandas.
In [3]:
Copied!
import geopandas as gpd
aoi_g = gpd.GeoSeries([aoi])
aoi_g.boundary.plot()
import geopandas as gpd
aoi_g = gpd.GeoSeries([aoi])
aoi_g.boundary.plot()
Out[3]:
<Axes: >
In [4]:
Copied!
aoi_g.set_crs(5514).explore()
aoi_g.set_crs(5514).explore()
Out[4]:
Make this Notebook Trusted to load map: File -> Trust Notebook
Výběr dat z Copernicus Data Space na základě AOI¶
Propojení s Copernicus Data Space poskytuje balíček eodag.
V našem případě se zaměříme na produkt L2A (snímky, které byly korigovány na atmosférické vlivy).
In [5]:
Copied!
from eodag import EODataAccessGateway
dag = EODataAccessGateway()
product = "S2_MSI_L2A"
dag.available_providers(product)
from eodag import EODataAccessGateway
dag = EODataAccessGateway()
product = "S2_MSI_L2A"
dag.available_providers(product)
Out[5]:
['cop_dataspace', 'creodias', 'planetary_computer', 'sara']
Prostrorový rozsah definující AOI musí být poskytnut ve WGS-84 (EPSG:4326).
In [6]:
Copied!
extent = tuple(aoi_g.set_crs(5514).to_crs(4326).total_bounds)
extent
extent = tuple(aoi_g.set_crs(5514).to_crs(4326).total_bounds)
extent
Out[6]:
(np.float64(13.80108831369076), np.float64(49.25553081328358), np.float64(13.851104863967173), np.float64(49.29378731553107))
V našem případě se budeme dotazovat služby Copernicus Data Space.
In [8]:
Copied!
dag.set_preferred_provider("cop_dataspace")
results = dag.search(
productType="S2_MSI_L2A",
start="2025-10-01",
end="2025-10-31",
geom=tuple(extent),
count=True
)
results.number_matched
dag.set_preferred_provider("cop_dataspace")
results = dag.search(
productType="S2_MSI_L2A",
start="2025-10-01",
end="2025-10-31",
geom=tuple(extent),
count=True
)
results.number_matched
Out[8]:
18
Hledáme data Sentinel-2 Level 2A s nízkým podílem mraků (méně než 20%).
In [12]:
Copied!
results_cc = results.filter_property(operator="le", cloudCover=20)
results_cc
results_cc = results.filter_property(operator="le", cloudCover=20)
results_cc
Out[12]:
| SearchResult (1) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
0 EOProduct(id=S2C_MSIL2A_20251029T101151_N0511_R022_T33UVQ_20251029T141217, provider=cop_dataspace)
|
Pro další kroky již budeme potřebovat nastavit přístupové údaje.
- Vytvoříme si učet: zde
- Vytvoříme konfigurační soubor
In [13]:
Copied!
import os
username = "???"
password = "???"
dag.update_providers_config(f"""
cop_dataspace:
download:
extract: False
output_dir: {os.environ['HOME']}/Downloads/sentinel2
delete_archive: False
auth:
credentials:
username: {username}
password: {password}
""")
import os
username = "???"
password = "???"
dag.update_providers_config(f"""
cop_dataspace:
download:
extract: False
output_dir: {os.environ['HOME']}/Downloads/sentinel2
delete_archive: False
auth:
credentials:
username: {username}
password: {password}
""")
Scény před stažením prohlédneme (zdroj):
In [14]:
Copied!
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
fig = plt.figure(figsize=(10, 8))
for i, product in enumerate(results_cc, start=1):
# This line takes care of downloading the quicklook
quicklook_path = product.get_quicklook()
# Plot the quicklook
img = mpimg.imread(quicklook_path)
ax = fig.add_subplot(3, 4, i)
ax.set_title(i - 1)
plt.imshow(img)
plt.tight_layout()
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
fig = plt.figure(figsize=(10, 8))
for i, product in enumerate(results_cc, start=1):
# This line takes care of downloading the quicklook
quicklook_path = product.get_quicklook()
# Plot the quicklook
img = mpimg.imread(quicklook_path)
ax = fig.add_subplot(3, 4, i)
ax.set_title(i - 1)
plt.imshow(img)
plt.tight_layout()
0.00B [00:00, ?B/s]
Stažení vybraného produktu Sentinel-2¶
In [15]:
Copied!
data_path = dag.download(results_cc[0])
data_path = dag.download(results_cc[0])
0.00B [00:00, ?B/s]
In [16]:
Copied!
from pathlib import Path
data_path, Path(data_path).exists()
from pathlib import Path
data_path, Path(data_path).exists()
Out[16]:
('/mnt/home/landamar/Downloads/sentinel2/S2C_MSIL2A_20251029T101151_N0511_R022_T33UVQ_20251029T141217.zip',
True)
Výpočet NDVI¶
Pro výpočet normalizavaného vegetačního indexu je potřeba červený a blízký infračervený kanál. V případě Sentinel-2 jde o 4 (červený) a 8 (blízký infračervený) ḱanál (viz wikipedia).
Příprava vstupních dat¶
Rastrová data otevřeme pomocí balíčku rasterio.
In [18]:
Copied!
import rasterio as rio
with rio.open(data_path) as ds:
sd = ds.subdatasets
meta = ds.meta
print(meta)
for idx in range(len(sd)):
print(idx, sd[idx])
import rasterio as rio
with rio.open(data_path) as ds:
sd = ds.subdatasets
meta = ds.meta
print(meta)
for idx in range(len(sd)):
print(idx, sd[idx])
{'driver': 'SENTINEL2', 'dtype': 'float_', 'nodata': None, 'width': 512, 'height': 512, 'count': 0, 'crs': None, 'transform': Affine(1.0, 0.0, 0.0,
0.0, 1.0, 0.0)}
0 SENTINEL2_L2A:/vsizip//mnt/home/landamar/Downloads/sentinel2/S2C_MSIL2A_20251029T101151_N0511_R022_T33UVQ_20251029T141217.zip/S2C_MSIL2A_20251029T101151_N0511_R022_T33UVQ_20251029T141217.SAFE/MTD_MSIL2A.xml:10m:EPSG_32633
1 SENTINEL2_L2A:/vsizip//mnt/home/landamar/Downloads/sentinel2/S2C_MSIL2A_20251029T101151_N0511_R022_T33UVQ_20251029T141217.zip/S2C_MSIL2A_20251029T101151_N0511_R022_T33UVQ_20251029T141217.SAFE/MTD_MSIL2A.xml:20m:EPSG_32633
2 SENTINEL2_L2A:/vsizip//mnt/home/landamar/Downloads/sentinel2/S2C_MSIL2A_20251029T101151_N0511_R022_T33UVQ_20251029T141217.zip/S2C_MSIL2A_20251029T101151_N0511_R022_T33UVQ_20251029T141217.SAFE/MTD_MSIL2A.xml:60m:EPSG_32633
3 SENTINEL2_L2A:/vsizip//mnt/home/landamar/Downloads/sentinel2/S2C_MSIL2A_20251029T101151_N0511_R022_T33UVQ_20251029T141217.zip/S2C_MSIL2A_20251029T101151_N0511_R022_T33UVQ_20251029T141217.SAFE/MTD_MSIL2A.xml:TCI:EPSG_32633
In [19]:
Copied!
with rio.open(sd[0]) as b10m:
print("CRS:", b10m.crs)
print("Bounds:", b10m.bounds)
print("Bands:", b10m.count)
with rio.open(sd[0]) as b10m:
print("CRS:", b10m.crs)
print("Bounds:", b10m.bounds)
print("Bands:", b10m.count)
CRS: EPSG:32633 Bounds: BoundingBox(left=399960.0, bottom=5390220.0, right=509760.0, top=5500020.0) Bands: 6
In [20]:
Copied!
with rio.open(sd[0]) as b10m:
meta = b10m.meta
for b in b10m.indexes:
print(b, b10m.tags(b)["BANDNAME"])
with rio.open(sd[0]) as b10m:
meta = b10m.meta
for b in b10m.indexes:
print(b, b10m.tags(b)["BANDNAME"])
1 B4 2 B3 3 B2 4 B8 5 AOT 6 WVP
Data lze vizualizovat pomocí modulu rasterio.plot.
In [23]:
Copied!
with rio.open(sd[0]) as b10m:
b4 = b10m.read(1)
from rasterio.plot import show
show(b4, cmap='pink')
with rio.open(sd[0]) as b10m:
b4 = b10m.read(1)
from rasterio.plot import show
show(b4, cmap='pink')
Out[23]:
<Axes: >
Úkol: Vyzkoušejte odstíny šedi (plt.cm.Greys_r) a navrhněte proceduru pro vyrovnání odstínů šedi dle histogramu.
Ořez vstupních dat pomocí AOI¶
Pro tento účel je třeba transformovat AOI do UTM 33N.
In [24]:
Copied!
with rio.open(sd[0]) as b10m:
print(b10m.crs)
with rio.open(sd[0]) as b10m:
print(b10m.crs)
EPSG:32633
In [25]:
Copied!
from shapely.geometry import box
aoi_utm = aoi_g.set_crs(5514).to_crs(epsg=b10m.crs.to_epsg())[0]
print(aoi_utm)
from shapely.geometry import box
aoi_utm = aoi_g.set_crs(5514).to_crs(epsg=b10m.crs.to_epsg())[0]
print(aoi_utm)
POLYGON ((414415.23495359445 5460782.02198694, 414665.4920849454 5460463.936878342, 414885.5097333678 5460224.785940796, 415023.2335143471 5460156.144412951, 415432.2788881825 5460031.120844939, 415609.3183309596 5459847.612494443, 415845.0967847256 5459746.653657803, 416170.5322425036 5459710.069922876, 416212.334955329 5459671.361150695, 416281.25354680244 5459650.531563609, 416329.1333147645 5459277.740589341, 416291.07707708055 5459057.6603579465, 416221.63009067485 5458888.716304019, 416031.1153676759 5458669.780376454, 415806.3081740029 5458260.0879329955, 415411.4216765355 5457985.713588415, 415423.1415911109 5457949.7642296525, 415615.2365594797 5457849.527558946, 415619.99107942317 5457740.016088234, 415691.44267991977 5457664.212432121, 415806.8168281668 5457577.933969953, 416188.20587965345 5457389.175878327, 416411.0493734572 5457326.424576229, 416263.583050246 5457221.166570061, 416164.0352813146 5457155.637038321, 415973.23765619425 5457240.222565783, 415614.59981618926 5457249.660854141, 415604.4127139483 5456978.945014593, 415573.7874442504 5456865.610744734, 415068.46528963407 5456599.92200434, 414951.6334712077 5456587.438216983, 414867.9808618081 5456534.830252225, 414715.2318538756 5456523.468669757, 414592.4164581152 5456580.796869376, 414487.05292553833 5456591.119101015, 414385.435096011 5456588.12975445, 414268.76314589765 5456546.76710763, 414259.123015705 5456641.536383862, 414162.85455595644 5456602.3615340665, 414135.22717111965 5456553.38001437, 413971.98748490604 5456654.925042947, 413931.9951672666 5456716.413096102, 413958.0402121966 5456869.367244104, 414042.2650516567 5456914.77437116, 414041.6177374156 5456981.371856687, 413856.08569724334 5457109.013166103, 413800.6970915466 5457256.7045930885, 413629.40400839155 5457385.5939385155, 413445.14661577495 5457421.787378978, 413283.92184041464 5457421.787586755, 413182.7618268375 5457443.992728189, 413228.0531508613 5457644.49841364, 413294.4814617025 5457780.790688897, 413288.6775453149 5457995.396959149, 413300.2560447471 5458084.638929602, 413105.8625515222 5458154.081693937, 412953.69286875275 5458175.662588881, 412846.9490903621 5458165.800941253, 412861.2236009652 5458507.731866156, 412794.73008523043 5458808.334991757, 412917.8413336353 5459032.088638158, 412917.3481903559 5459305.391025772, 413025.5740182328 5459314.111251961, 413148.26036700106 5459283.366625501, 413242.1985648548 5459437.491191088, 413185.8996477725 5459543.948561694, 413211.28207764484 5459595.527583643, 413265.26083709934 5459705.21602713, 413281.05852440244 5459837.1970285475, 413342.5726058419 5459882.065009499, 413484.4718724403 5459847.1471826285, 413675.2893589437 5460056.544068752, 413903.3941448833 5460266.966507941, 414038.08098894905 5460668.975665191, 414098.89451123273 5460784.9055157695, 414415.23495359445 5460782.02198694))
Rastrová data lze jednoduše ořezat pomocí modulu rasterio.mask.
In [26]:
Copied!
import rasterio.mask
with rio.open(sd[0]) as b10m:
out_image, out_transform = rasterio.mask.mask(b10m, [aoi_utm], crop=True)
out_image.shape
import rasterio.mask
with rio.open(sd[0]) as b10m:
out_image, out_transform = rasterio.mask.mask(b10m, [aoi_utm], crop=True)
out_image.shape
Out[26]:
(6, 427, 363)
Pro učel ořezu rastrových dat vytvoříme funkci, kromě toho z důvodu výpočtu NDVI:
- rastrové hodnoty převedeme na čísla s plovoucí desetinou čárkou,
- hodnotu nula nahradíme za NaN (Not a Number)
In [27]:
Copied!
import numpy as np
def clip_raster(path, aoi):
with rio.open(data_path) as ds:
with rasterio.open(ds.subdatasets[0]) as b10m:
image, transform = rasterio.mask.mask(b10m, [aoi], crop=True)
data = image.astype(float)
data[data==0] = np.nan
return data, transform
data_aoi, clip_transform = clip_raster(data_path, aoi_utm)
data_aoi.shape, clip_transform
import numpy as np
def clip_raster(path, aoi):
with rio.open(data_path) as ds:
with rasterio.open(ds.subdatasets[0]) as b10m:
image, transform = rasterio.mask.mask(b10m, [aoi], crop=True)
data = image.astype(float)
data[data==0] = np.nan
return data, transform
data_aoi, clip_transform = clip_raster(data_path, aoi_utm)
data_aoi.shape, clip_transform
/opt/jupyterhub-venv/lib/python3.12/site-packages/rasterio/__init__.py:356: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned. dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
Out[27]:
((6, 427, 363),
Affine(10.0, 0.0, 412790.0,
0.0, -10.0, 5460790.0))
Výpočet NDVI¶
Výpočet bude proveden pomocí balíčku NumPy.
In [28]:
Copied!
import numpy as np
b8_aoi = data_aoi[3]
b4_aoi = data_aoi[0]
ndvi = (b8_aoi - b4_aoi) / (b8_aoi + b4_aoi)
np.nanmin(ndvi), np.nanmax(ndvi)
import numpy as np
b8_aoi = data_aoi[3]
b4_aoi = data_aoi[0]
ndvi = (b8_aoi - b4_aoi) / (b8_aoi + b4_aoi)
np.nanmin(ndvi), np.nanmax(ndvi)
Out[28]:
(np.float64(-0.06702127659574468), np.float64(0.7226027397260274))
Vizualizace NDVI¶
In [29]:
Copied!
show(ndvi, cmap='Greens')
show(ndvi, cmap='Greens')
Out[29]:
<Axes: >
Reklasifikace NDVI¶
In [31]:
Copied!
ndvi[ndvi < 0.1] = 1
ndvi[ndvi < 0.5] = 2
ndvi[ndvi < 1.0] = 3
np.unique(ndvi, return_counts=True)
ndvi[ndvi < 0.1] = 1
ndvi[ndvi < 0.5] = 2
ndvi[ndvi < 1.0] = 3
np.unique(ndvi, return_counts=True)
Out[31]:
(array([ 1., 2., 3., nan]), array([15938, 64139, 15421, 59503]))
In [32]:
Copied!
show(ndvi, cmap='Greens')
show(ndvi, cmap='Greens')
Out[32]:
<Axes: >
Uložení reklasifikovaného NDVI do souboru GeoTIFF¶
In [35]:
Copied!
with rio.open(sd[0]) as b10m:
meta = b10m.meta
meta.update({
"driver": 'GTiff',
"dtype": "int8",
"height": ndvi.shape[0],
"width": ndvi.shape[1],
"transform": clip_transform,
"nodata": 0
})
np.nan_to_num(ndvi, copy=False, nan=0)
output_path = '/tmp/ndvi_cat.tif'
with rasterio.open(output_path, 'w', **meta) as dst:
dst.write_band(1, ndvi.astype(rasterio.int8))
with rio.open(sd[0]) as b10m:
meta = b10m.meta
meta.update({
"driver": 'GTiff',
"dtype": "int8",
"height": ndvi.shape[0],
"width": ndvi.shape[1],
"transform": clip_transform,
"nodata": 0
})
np.nan_to_num(ndvi, copy=False, nan=0)
output_path = '/tmp/ndvi_cat.tif'
with rasterio.open(output_path, 'w', **meta) as dst:
dst.write_band(1, ndvi.astype(rasterio.int8))
In [36]:
Copied!
Path(output_path).exists()
Path(output_path).exists()
Out[36]:
True