Zpracování obrazových dat pomocí Rasterio a Sentinelsat¶
Sentinel-2¶
Data Sentinel-2 můžeme stáhnout buď pomocí Sentinel Hub (https://apps.sentinel-hub.com/eo-browser/) anebo Python balíčku sentinelsat (https://sentinelsat.readthedocs.io/).
Zájmové území (AOI)¶
Můžeme použít vrstvu obcí dostupnou z datasetu Data250: PolbndMunDA
import fiona
nuts_name = "Katovice"
with fiona.open('/home/martin/geodata/data250/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])
(CRS.from_epsg(5514), 76)
Zobrazení AOI¶
Na tomto místě můžeme použít knihovnu GeoPandas.
aoi_map = aoi_g.set_crs(epsg=5514)
aoi_map.explore()
Výběr dat z Copernicus Open Access Hub na základě AOI¶
Propojení s Copernicus Open Access Hub poskytuje balíček sentinelsat.
from sentinelsat import SentinelAPI
from datetime import date
user = '???'
password = '???'
api = SentinelAPI(user, password, 'https://scihub.copernicus.eu/dhus/')
Polygon definující AOI musí být poskytnut ve WGS-84 (EPSG:4326) a v reprezentaci Well Known Text (WKT).
from shapely.geometry import shape
shape(aoi)
Transformace z S-JTSK (EPSG:5514) do WGS-84 může být provedena pomocí balíčku pyproj e. Pořadí os v S-JTSK (yx) vyžaduje zadat parametr always_xy
.
import pyproj
from shapely.ops import transform
wgs84 = pyproj.CRS('EPSG:4326')
project = pyproj.Transformer.from_crs(nuts_crs, wgs84, always_xy=True).transform
aoi_wgs84 = transform(project, shape(aoi))
aoi_wgs84.wkt
'POLYGON ((13.82293979970452 49.29378731553108, 13.826448782496241 49.29096169392838, 13.829524806581228 49.28884162799265, 13.831432930640515 49.288243481518464, 13.837083190261211 49.28717584777203, 13.839555926014192 49.285549971082425, 13.842818636329449 49.28467452564179, 13.847300290005018 49.28439023825107, 13.847883089119554 49.28404784879536, 13.848834910125367 49.28386996490329, 13.849571139469106 49.280523878183764, 13.849094061927076 49.278539406747605, 13.8481748314762 49.27701051832525, 13.845602059544255 49.27501541209587, 13.842598305742673 49.27129998774989, 13.837229015713076 49.26877794700398, 13.837397682922232 49.26845626315535, 13.840058807437863 49.26758133662277, 13.840147237565011 49.26659711734293, 13.841145138743796 49.26592524208866, 13.842748811605027 49.265165202130724, 13.848029484206469 49.263520011347346, 13.851104863967173 49.26298615936998, 13.849100447729283 49.26201936256567, 13.847746245382487 49.26141639429876, 13.845106699748191 49.2621509239092, 13.840176524288749 49.26218643372698, 13.840093607659867 49.25975038260884, 13.83969669210826 49.25872690033793, 13.832809709372693 49.25626752232297, 13.831207124300258 49.25613902126496, 13.830068948104621 49.2556542654965, 13.827972667568863 49.25553081328358, 13.82627302009805 49.25602925210631, 13.82482316084171 49.256107365446695, 13.823427607921115 49.256066269498184, 13.821833434944253 49.25567794432083, 13.821680693454097 49.256528882975715, 13.820366387203523 49.25616307224135, 13.819997299443186 49.2557186896414, 13.817732653109989 49.25660898263055, 13.817169953041978 49.25715633794693, 13.817494942406096 49.258535563268055, 13.818642467582018 49.258955763393104, 13.818619275509771 49.25955460334651, 13.816042498186166 49.26067641960492, 13.815249610731733 49.26199684758112, 13.812868025939906 49.26313182141011, 13.81032822281889 49.263431272012234, 13.808112737292614 49.26340843746765, 13.806717822896866 49.263593783545666, 13.807296733361227 49.26540340856068, 13.808180074698482 49.266638540280795, 13.808053815937088 49.268567724980855, 13.808193604614553 49.269371941705415, 13.805506930228514 49.26996886535398, 13.803410908853504 49.270141305244294, 13.801946031489013 49.27003741310387, 13.802067744317371 49.273114508418786, 13.80108831369076 49.275808417114746, 13.802731765273963 49.277838229983, 13.80266547919254 49.280296026369776, 13.80415128263931 49.28038985737549, 13.805844441675236 49.28013080789998, 13.80710230083473 49.28153022431773, 13.806305269650135 49.28247962563239, 13.806643003849969 49.282947091692314, 13.80736125697369 49.28394120629608, 13.80754980674181 49.28513038216314, 13.808385754035024 49.28554261487919, 13.810344104922914 49.28524869420771, 13.812922187580014 49.287158831226826, 13.81601290222796 49.289083389168965, 13.817778218299225 49.29271771669948, 13.81858947183848 49.29376885463424, 13.82293979970452 49.29378731553108))'
Hledáme data Sentinel-2 Level 2A s nízkým podílem mraků (méně než 25%).
products = api.query(aoi_wgs84.wkt,
date=(date(2023, 10, 1), date(2023, 10, 29)),
platformname='Sentinel-2',
producttype='S2MSI2A',
cloudcoverpercentage=(0, 25))
df = api.to_dataframe(products)
df
title | link | link_alternative | link_icon | summary | ondemand | generationdate | beginposition | endposition | ingestiondate | ... | s2datatakeid | producttype | platformidentifier | orbitdirection | platformserialidentifier | processinglevel | datastripidentifier | granuleidentifier | identifier | uuid | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
43e5a483-cc10-4f6b-a8b9-fa8a60a04c6b | S2B_MSIL2A_20231025T101029_N0509_R022_T33UVQ_2... | https://scihub.copernicus.eu/dhus/odata/v1/Pro... | https://scihub.copernicus.eu/dhus/odata/v1/Pro... | https://scihub.copernicus.eu/dhus/odata/v1/Pro... | Date: 2023-10-25T10:10:29.024Z, Instrument: MS... | false | 2023-10-25 11:46:51 | 2023-10-25 10:10:29.024 | 2023-10-25 10:10:29.024 | 2023-10-25 15:49:30.431 | ... | GS2B_20231025T101029_034654_N05.09 | S2MSI2A | 2017-013A | DESCENDING | Sentinel-2B | Level-2A | S2B_OPER_MSI_L2A_DS_2BPS_20231025T114651_S2023... | S2B_OPER_MSI_L2A_TL_2BPS_20231025T114651_A0346... | S2B_MSIL2A_20231025T101029_N0509_R022_T33UVQ_2... | 43e5a483-cc10-4f6b-a8b9-fa8a60a04c6b |
e086b61e-6362-4bd4-9f4c-b1b78ce9ebb6 | S2B_MSIL2A_20231002T095739_N0509_R122_T33UVQ_2... | https://scihub.copernicus.eu/dhus/odata/v1/Pro... | https://scihub.copernicus.eu/dhus/odata/v1/Pro... | https://scihub.copernicus.eu/dhus/odata/v1/Pro... | Date: 2023-10-02T09:57:39.024Z, Instrument: MS... | false | 2023-10-02 13:00:29 | 2023-10-02 09:57:39.024 | 2023-10-02 09:57:39.024 | 2023-10-02 17:07:21.232 | ... | GS2B_20231002T095739_034325_N05.09 | S2MSI2A | 2017-013A | DESCENDING | Sentinel-2B | Level-2A | S2B_OPER_MSI_L2A_DS_2BPS_20231002T130029_S2023... | S2B_OPER_MSI_L2A_TL_2BPS_20231002T130029_A0343... | S2B_MSIL2A_20231002T095739_N0509_R122_T33UVQ_2... | e086b61e-6362-4bd4-9f4c-b1b78ce9ebb6 |
2 rows × 42 columns
Vybere první nalezený produkt.
df.iloc[0]
title S2B_MSIL2A_20231025T101029_N0509_R022_T33UVQ_2... link https://scihub.copernicus.eu/dhus/odata/v1/Pro... link_alternative https://scihub.copernicus.eu/dhus/odata/v1/Pro... link_icon https://scihub.copernicus.eu/dhus/odata/v1/Pro... summary Date: 2023-10-25T10:10:29.024Z, Instrument: MS... ondemand false generationdate 2023-10-25 11:46:51 beginposition 2023-10-25 10:10:29.024000 endposition 2023-10-25 10:10:29.024000 ingestiondate 2023-10-25 15:49:30.431000 orbitnumber 34654 relativeorbitnumber 22 illuminationazimuthangle 171.809081 illuminationzenithangle 61.583316 vegetationpercentage 51.092267 notvegetatedpercentage 11.333832 waterpercentage 1.115555 unclassifiedpercentage 1.21178 mediumprobacloudspercentage 8.707906 highprobacloudspercentage 9.928188 snowicepercentage 0.000788 cloudcoverpercentage 21.617395 level1cpdiidentifier S2B_OPER_MSI_L1C_TL_2BPS_20231025T110431_A0346... gmlfootprint <gml:Polygon srsName="http://www.opengis.net/g... footprint MULTIPOLYGON (((13.641546225080209 48.65702077... format SAFE processingbaseline 05.09 platformname Sentinel-2 filename S2B_MSIL2A_20231025T101029_N0509_R022_T33UVQ_2... instrumentname Multi-Spectral Instrument instrumentshortname MSI size 1.08 GB s2datatakeid GS2B_20231025T101029_034654_N05.09 producttype S2MSI2A platformidentifier 2017-013A orbitdirection DESCENDING platformserialidentifier Sentinel-2B processinglevel Level-2A datastripidentifier S2B_OPER_MSI_L2A_DS_2BPS_20231025T114651_S2023... granuleidentifier S2B_OPER_MSI_L2A_TL_2BPS_20231025T114651_A0346... identifier S2B_MSIL2A_20231025T101029_N0509_R022_T33UVQ_2... uuid 43e5a483-cc10-4f6b-a8b9-fa8a60a04c6b Name: 43e5a483-cc10-4f6b-a8b9-fa8a60a04c6b, dtype: object
Stažení vybraného produktu Sentinel¶
data = api.download(df.iloc[0].uuid)
data
Widget Javascript not detected. It may not be installed or enabled properly. Reconnecting the current kernel may help.
Widget Javascript not detected. It may not be installed or enabled properly. Reconnecting the current kernel may help.
{'id': '43e5a483-cc10-4f6b-a8b9-fa8a60a04c6b', 'title': 'S2B_MSIL2A_20231025T101029_N0509_R022_T33UVQ_20231025T114651', 'size': 1162042944, 'md5': '2f27ba8de74c1ed73f9d6f0b047ccf4a', 'date': datetime.datetime(2023, 10, 25, 10, 10, 29, 24000), 'footprint': 'POLYGON((15.132795357138537 48.75548909987771,15.130892980718372 48.75143072455255,15.09032559893374 48.664730447818165,13.641546225080209 48.6570207750582,13.614272368185702 49.64443007346099,15.13521359578094 49.652643868751746,15.132795357138537 48.75548909987771))', 'url': "https://scihub.copernicus.eu/dhus/odata/v1/Products('43e5a483-cc10-4f6b-a8b9-fa8a60a04c6b')/$value", 'Online': True, 'Creation Date': datetime.datetime(2023, 10, 25, 15, 50, 10, 60000), 'Ingestion Date': datetime.datetime(2023, 10, 25, 15, 49, 30, 431000), 'manifest_name': 'manifest.safe', 'product_root_dir': 'S2B_MSIL2A_20231025T101029_N0509_R022_T33UVQ_20231025T114651.SAFE', 'quicklook_url': "https://scihub.copernicus.eu/dhus/odata/v1/Products('43e5a483-cc10-4f6b-a8b9-fa8a60a04c6b')/Products('Quicklook')/$value", 'path': 'S2B_MSIL2A_20231025T101029_N0509_R022_T33UVQ_20231025T114651.zip', 'downloaded_bytes': 1162042944}
import zipfile
with zipfile.ZipFile(data['path']) as zip:
zip.extractall()
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¶
from pathlib import Path
b4_path = list(Path(data['product_root_dir']).rglob('*B04_10m.jp2'))[0]
b8_path = list(Path(data['product_root_dir']).rglob('*B08_10m.jp2'))[0]
b4_path, b8_path
(PosixPath('S2B_MSIL2A_20231025T101029_N0509_R022_T33UVQ_20231025T114651.SAFE/GRANULE/L2A_T33UVQ_A034654_20231025T101025/IMG_DATA/R10m/T33UVQ_20231025T101029_B04_10m.jp2'), PosixPath('S2B_MSIL2A_20231025T101029_N0509_R022_T33UVQ_20231025T114651.SAFE/GRANULE/L2A_T33UVQ_A034654_20231025T101025/IMG_DATA/R10m/T33UVQ_20231025T101029_B08_10m.jp2'))
Rastrová data otevřeme pomocí balíčku rasterio.
import rasterio
b4 = rasterio.open(b4_path)
b4
<open DatasetReader name='S2B_MSIL2A_20231025T101029_N0509_R022_T33UVQ_20231025T114651.SAFE/GRANULE/L2A_T33UVQ_A034654_20231025T101025/IMG_DATA/R10m/T33UVQ_20231025T101029_B04_10m.jp2' mode='r'>
print("CRS:", b4.crs)
print("Bounds:", b4.bounds)
print("Bands:", b4.count)
CRS: EPSG:32633 Bounds: BoundingBox(left=399960.0, bottom=5390220.0, right=509760.0, top=5500020.0) Bands: 1
Data lze vizualiyovat pomocí raster.plot.
Ořez vstupních dat pomocí AOI¶
Pro tento účel je třeba transformovat AOI do UTM 33N.
utm = pyproj.CRS(b4.crs)
project = pyproj.Transformer.from_crs(nuts_crs, utm, always_xy=True).transform
aoi_utm = transform(project, shape(aoi))
show(b4, cmap='pink')
import matplotlib as mpl
ax = mpl.pyplot.gca()
# TODO: https://gis.stackexchange.com/questions/193653/plot-shapefile-on-top-of-raster-using-plot-and-imshow-from-matplotlib
#from descartes import PolygonPatch
#patches = [PolygonPatch(aoi_utm, edgecolor="red", facecolor="none", linewidth=2)]
#ax.add_collection(mpl.collections.PatchCollection(patches, match_original=True))
Rastrová data lze jednoduše ořezat pomocí modulu rasterio.mask.
import rasterio.mask
with rasterio.open(b4_path) as src:
out_image, out_transform = rasterio.mask.mask(src, [aoi_utm], crop=True)
out_image.shape
(1, 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)
import numpy as np
def clip_raster(path):
with rasterio.open(path) as src:
image, transform = rasterio.mask.mask(src, [aoi_utm], crop=True)
data = image.astype(float)[0]
data[data==0] = np.nan
return data, transform
b4_aoi, clip_transform = clip_raster(b4_path)
b8_aoi, _ = clip_raster(b8_path)
b4_aoi.shape, clip_transform
((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.
import numpy as np
ndvi = (b8_aoi - b4_aoi) / (b8_aoi + b4_aoi)
np.nanmin(ndvi), np.nanmax(ndvi)
(-0.0654988575780655, 0.6845070422535211)
Vizualizace NDVI¶
Reklasifikace NDVI¶
ndvi[ndvi < 0.1] = 1
ndvi[ndvi < 0.5] = 2
ndvi[ndvi < 1.0] = 3
np.unique(ndvi, return_counts=True)
(array([ 1., 2., 3., nan]), array([13147, 61714, 20637, 59503]))
Uložení reklasifikovaného NDVI do souboru GeoTIFF¶
meta = b4.meta
meta["driver"] = 'GTiff'
meta["dtype"] = "int8"
meta["height"] = ndvi.shape[0]
meta["width"] = ndvi.shape[1]
meta["transform"] = clip_transform
np.nan_to_num(ndvi, copy=False, nan=0)
meta["nodata"] = 0
with rasterio.open('ndvi_cat.tif', 'w', **meta) as dst:
dst.write_band(1, ndvi.astype(rasterio.int8))