Zpracování obrazových dat pomocí Rasterio a Sentinelsat¶
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 buď 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)/Default.aspx?mode=TextMeta&side=mapy_data250&metadataID=CZ-CUZK-DATA250-V&head_tab=sekce-02-gp&menu=22910) (dokumentace):
In [1]:
Copied!
import fiona
nuts_name = "Katovice"
with fiona.open('/home/martin/geodata/vyuka/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, len(aoi.coordinates[0])
import fiona
nuts_name = "Katovice"
with fiona.open('/home/martin/geodata/vyuka/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, len(aoi.coordinates[0])
Out[1]:
(CRS.from_epsg(5514), 76)
Zobrazení AOI¶
Na tomto místě můžeme použít knihovnu GeoPandas.
In [2]:
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[2]:
<Axes: >
In [25]:
Copied!
aoi_map = aoi_g.set_crs(epsg=5514)
aoi_map.explore()
aoi_map = aoi_g.set_crs(epsg=5514)
aoi_map.explore()
Out[25]:
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 [4]:
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)
/home/martin/git/k155cvut/uzpr/venv/lib/python3.10/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html from .autonotebook import tqdm as notebook_tqdm wekeo skipped: could not be loaded from user configuration
Out[4]:
['astraea_eod', 'cop_dataspace', 'creodias', 'onda', 'planetary_computer', 'sara']
Prostrorový rozsah definující AOI musí být poskytnut ve WGS-84 (EPSG:4326).
In [5]:
Copied!
extent = tuple(aoi_map.to_crs(4326).total_bounds)
extent
extent = tuple(aoi_map.to_crs(4326).total_bounds)
extent
Out[5]:
(np.float64(13.80108831369076), np.float64(49.25553081328358), np.float64(13.851104863967173), np.float64(49.29378731553108))
V našem případě se budeme dotazovat služby Copernicus Data Space.
In [6]:
Copied!
dag.set_preferred_provider("cop_dataspace")
results = dag.search(
productType="S2_MSI_L2A",
start="2024-10-01",
end="2024-10-23",
geom=tuple(extent),
count=True
)
print(results.number_matched)
dag.set_preferred_provider("cop_dataspace")
results = dag.search(
productType="S2_MSI_L2A",
start="2024-10-01",
end="2024-10-23",
geom=tuple(extent),
count=True
)
print(results.number_matched)
5
Hledáme data Sentinel-2 Level 2A s nízkým podílem mraků (méně než 25%).
In [7]:
Copied!
results_cc = results.filter_property(operator="le", cloudCover=5)
results_cc
results_cc = results.filter_property(operator="le", cloudCover=5)
results_cc
Out[7]:
| SearchResult (1) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
0 EOProduct(id=S2B_MSIL2A_20241009T100829_N0511_R022_T33UVQ_20241009T132806, 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 [8]:
Copied!
username = "???"
password = "???"
dag.update_providers_config(f"""
cop_dataspace:
download:
extract: False
outputs_prefix: /tmp/sentinel
delete_archive: False
auth:
credentials:
username: {username}
password: {password}
""")
username = "???"
password = "???"
dag.update_providers_config(f"""
cop_dataspace:
download:
extract: False
outputs_prefix: /tmp/sentinel
delete_archive: False
auth:
credentials:
username: {username}
password: {password}
""")
Scény před stažením prohlédneme (zdroj):
In [9]:
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()
quicklooks/S2B_MSIL2A_20241009T100829_N0511_R022_T33UVQ_20241009T132806: 0.00B [00:00, ?B/s]
Stažení vybraného produktu Sentinel-2¶
In [10]:
Copied!
data_path = dag.download(results_cc[0])
data_path = dag.download(results_cc[0])
0.00B [00:00, ?B/s] S2B_MSIL2A_20241009T100829_N0511_R022_T33UVQ_20241009T132806: 0.00B [00:00, ?B/s] S2B_MSIL2A_20241009T100829_N0511_R022_T33UVQ_20241009T132806: 0file [00:00, ?file/s] S2B_MSIL2A_20241009T100829_N0511_R022_T33UVQ_20241009T132806: 100%|██████████████████████| 1/1 [00:00<00:00, 658.65file/s]
In [11]:
Copied!
from pathlib import Path
print(data_path, Path(data_path).exists())
from pathlib import Path
print(data_path, Path(data_path).exists())
/tmp/S2B_MSIL2A_20241009T100829_N0511_R022_T33UVQ_20241009T132806.SAFE.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 [12]:
Copied!
import rasterio as rio
with rio.open(data_path) as ds:
sd = ds.subdatasets
meta = ds.meta
print(meta)
print(len(sd), sd[0])
import rasterio as rio
with rio.open(data_path) as ds:
sd = ds.subdatasets
meta = ds.meta
print(meta)
print(len(sd), sd[0])
{'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)}
4 SENTINEL2_L2A:/vsizip//tmp/S2B_MSIL2A_20241009T100829_N0511_R022_T33UVQ_20241009T132806.SAFE.zip/S2B_MSIL2A_20241009T100829_N0511_R022_T33UVQ_20241009T132806.SAFE/MTD_MSIL2A.xml:10m:EPSG_32633
/home/martin/git/k155cvut/uzpr/venv/lib/python3.10/site-packages/rasterio/__init__.py:355: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned. dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
In [13]:
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 [14]:
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 vizualiyovat pomocí raster.plot.
In [15]:
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[15]:
<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 [16]:
Copied!
with rio.open(sd[0]) as b10m:
print(b10m.crs)
with rio.open(sd[0]) as b10m:
print(b10m.crs)
EPSG:32633
In [17]:
Copied!
from shapely.geometry import box
aoi_utm = aoi_map.to_crs(epsg=32633)[0]
print(aoi_utm)
from shapely.geometry import box
aoi_utm = aoi_map.to_crs(epsg=32633)[0]
print(aoi_utm)
POLYGON ((414415.2349535946 5460782.021986941, 414665.4920849456 5460463.936878342, 414885.5097333678 5460224.785940798, 415023.23351434694 5460156.144412949, 415432.2788881825 5460031.120844939, 415609.3183309595 5459847.612494443, 415845.0967847256 5459746.653657802, 416170.5322425035 5459710.069922877, 416212.334955329 5459671.361150695, 416281.25354680244 5459650.531563609, 416329.1333147645 5459277.740589341, 416291.07707708055 5459057.6603579465, 416221.63009067485 5458888.716304019, 416031.1153676757 5458669.780376454, 415806.30817400274 5458260.0879329955, 415411.4216765354 5457985.713588415, 415423.1415911108 5457949.7642296525, 415615.2365594797 5457849.527558946, 415619.9910794232 5457740.016088234, 415691.44267991977 5457664.212432121, 415806.81682816683 5457577.933969952, 416188.2058796536 5457389.175878327, 416411.04937345727 5457326.42457623, 416263.583050246 5457221.166570061, 416164.0352813146 5457155.637038321, 415973.2376561943 5457240.222565784, 415614.5998161894 5457249.660854142, 415604.4127139484 5456978.945014593, 415573.7874442503 5456865.610744733, 415068.4652896342 5456599.92200434, 414951.63347120787 5456587.438216983, 414867.9808618081 5456534.830252225, 414715.2318538756 5456523.468669757, 414592.4164581152 5456580.796869377, 414487.05292553833 5456591.119101014, 414385.43509601103 5456588.12975445, 414268.76314589765 5456546.76710763, 414259.123015705 5456641.536383862, 414162.85455595655 5456602.3615340665, 414135.22717111954 5456553.38001437, 413971.98748490604 5456654.925042947, 413931.9951672668 5456716.413096102, 413958.0402121966 5456869.367244104, 414042.2650516567 5456914.77437116, 414041.6177374156 5456981.371856687, 413856.08569724346 5457109.0131661035, 413800.6970915465 5457256.704593089, 413629.40400839155 5457385.5939385155, 413445.14661577507 5457421.787378978, 413283.9218404145 5457421.787586757, 413182.7618268375 5457443.9927281905, 413228.0531508612 5457644.49841364, 413294.4814617025 5457780.790688897, 413288.6775453147 5457995.396959147, 413300.2560447472 5458084.638929601, 413105.8625515222 5458154.081693937, 412953.69286875287 5458175.662588881, 412846.9490903622 5458165.800941253, 412861.2236009653 5458507.731866157, 412794.73008523043 5458808.334991757, 412917.84133363544 5459032.088638159, 412917.3481903559 5459305.391025773, 413025.57401823293 5459314.111251961, 413148.26036700094 5459283.366625501, 413242.1985648548 5459437.491191088, 413185.8996477725 5459543.948561694, 413211.2820776447 5459595.527583643, 413265.2608370992 5459705.21602713, 413281.05852440244 5459837.197028547, 413342.572605842 5459882.065009499, 413484.4718724403 5459847.1471826285, 413675.2893589437 5460056.544068752, 413903.3941448833 5460266.966507941, 414038.08098894905 5460668.975665191, 414098.89451123273 5460784.9055157695, 414415.2349535946 5460782.021986941))
Rastrová data lze jednoduše ořezat pomocí modulu rasterio.mask.
In [18]:
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[18]:
(6, 427, 363)
Pro učel ořezu rastrových dat vytvoříme funkci, kromě toho:
- rastrové hodnoty převedeme na čísla s plovoucí desetinou čárkou,
- hodnotu nula nahradíme za NaN (Not a Number)
In [19]:
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
/home/martin/git/k155cvut/uzpr/venv/lib/python3.10/site-packages/rasterio/__init__.py:355: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned. dataset = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
Out[19]:
((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 [20]:
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[20]:
(np.float64(-0.07288401253918496), np.float64(0.6990387802452768))
Vizualizace NDVI¶
In [21]:
Copied!
show(ndvi, cmap='Greens')
show(ndvi, cmap='Greens')
Out[21]:
<Axes: >
Reklasifikace NDVI¶
In [22]:
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[22]:
(array([ 1., 2., 3., nan]), array([ 6605, 67179, 21714, 59503]))
In [23]:
Copied!
show(ndvi, cmap='Greens')
show(ndvi, cmap='Greens')
Out[23]:
<Axes: >
Uložení reklasifikovaného NDVI do souboru GeoTIFF¶
In [24]:
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)
with rasterio.open('/tmp/ndvi_cat.tif', '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)
with rasterio.open('/tmp/ndvi_cat.tif', 'w', **meta) as dst:
dst.write_band(1, ndvi.astype(rasterio.int8))